Context Navigation

trends.pro @ 108

Last change on this file since 108 was 108, checked in by ericg, 16 years ago
Bug in computation of 1D stddev in trends.pro corrected
File size: 15.2 KB

Line
1	FUNCTION trends, fld, trend_int, type, silent = silent
2
3	@common
4	@com_eg
5
6	; fld = field to work on
7	; trend_int = 1,2,3,4,...9 (see cases after)
8	; type = 't', 'xt', 'yt', 'zt' or 'xyt'
9
10	IF NOT keyword_set(silent) THEN print, ' Enter trends.pro'
11
12	trend_type = strtrim(string(trend_int), 2)
13
14	; select subset of time serie according to [sel1,sel2,freq]
15	; =========================================================
16
17	month1 = 11
18	month2 = 12
19	sample = 0
20	field_intl = field_int
21
22	CASE trend_type OF
23	'7': timeunit = 30/6
24	'8': timeunit = 30
25	'71': BEGIN
26	trend_type = '4'+strtrim(string(6*(month2-month1+1)), 2)
27	timeunit = 30/6
28	END
29	'81': BEGIN
30	trend_type = '4'+strtrim(string(month2-month1+1), 2)
31	timeunit = 30
32	END
33	ELSE: timeunit = 1
34	ENDCASE
35
36	delt = 30/timeunit
37	sel1 = (month1-1)*delt
38	sel2 = (month2)*delt-1
39	freq = 360/timeunit
40
41	; build indexes to keep in [0,N-1] = serie of [sel1,sel2] mod(freq)
42
43	IF (size(fld))[0] EQ 1 AND sample EQ 1 THEN BEGIN ; 1D time serie
44	IF NOT keyword_set(silent) THEN print, ' Extract months : ', month1, month2
45	Ndim = (size(fld))[1]
46	idx0 = floor(findgen(sel2-sel1+1)+sel1)
47	ndelt = n_elements(idx0)
48	idx = idx0
49	ninter = Ndim/freq
50	FOR t = 1, ninter-1 DO BEGIN
51	idx0 = idx0+freq
52	idx = [idx, idx0]
53	ENDFOR
54
55	; new array
56
57	fld = fld[idx]
58	jpt = n_elements(fld)
59
60	; change time accordingly
61
62	time = time[idx(0)]+(findgen(jpt)-1)*360/ndelt
63
64	ENDIF
65
66	; normalise 1d time series if c_normal=1
67
68	IF c_normal EQ 1 AND (size(fld))[0] EQ 1 THEN BEGIN
69	IF NOT keyword_set(silent) THEN print, ' Normalizing 1D time series with stddev = ',sqrt((moment(fld))[1])
70	fld = fld/sqrt((moment(fld))[1])
71	ENDIF
72
73	; modify time serie according to trend_typ
74	; ========================================
75
76	CASE strmid(trend_type, 0, 1) OF
77	'1': BEGIN ; trend = remove first value
78	; ---------------------------
79	IF (size(fld))[0] EQ 1 THEN BEGIN ; 1D time serie
80	fld = fld-fld[0]
81	ENDIF ELSE BEGIN ; hovmoeller
82	fldrem = fld[*, 0]#replicate(1, (size(fld))[2])
83	fld = fld-fldrem
84	ENDELSE
85	END
86	'2': BEGIN ; drift = remove previous value
87	; ------------------------------
88	IF (size(fld))[0] EQ 1 THEN BEGIN ; 1D time serie
89	fldrem = shift(fld, 1)
90	fldrem[0]= fld[0]
91	fld = fld-fldrem
92	ENDIF ELSE BEGIN ; hovmoeller
93	fldrem = shift(fld, 0, 1)
94	fldrem[, 0]= fld[, 0]
95	fld = fld-fldrem
96	ENDELSE
97	END
98	'3': BEGIN ; inverse trend = remove mean of <n> last values
99	; ----------------------------------------------
100	IF strlen(trend_type) GT 1 THEN $
101	mean_n = long(strmid(trend_type, 1, strlen(trend_type)-1)) ELSE mean_n = 1
102	IF (size(fld))[0] EQ 1 THEN BEGIN ; 1D time serie
103	fld = mean(fld[(size(fld))[1]-mean_n:(size(fld))[1]-1])-fld
104	ENDIF ELSE BEGIN ; hovmoeller
105	fldrem = reform(fld[*, (size(fld))[2]-mean_n:(size(fld))[2]-1],(size(fld))[1],mean_n)
106	fldrem = total(fldrem, 2)/float(mean_n)#replicate(1, (size(fld))[2])
107	idx = where(fld EQ valmask)
108	fld = fldrem-fld
109	IF idx(0) NE -1 THEN fld(idx) = valmask
110	ENDELSE
111	END
112	'4': BEGIN ; anomaly = remove average running mean
113	; -------------------------------------
114	IF strlen(trend_type) GT 1 THEN $
115	running = long(strmid(trend_type, 1, strlen(trend_type)-1)) ELSE running = 1
116
117	IF (size(fld))[0] EQ 1 THEN BEGIN ; 1D time serie
118	; build array to remove from time serie
119	IF running GE 2 THEN BEGIN
120	lenght = (size(fld))[1]
121	fldrem = fltarr(running)
122	FOR t = 0, running-1 DO BEGIN
123	fldrem(t) = mean(fld(long(findgen(lenght/running)*running+t)))
124	ENDFOR
125	IF fld_flag EQ 1 THEN BEGIN
126	fldrem_t1 = fldrem
127	ENDIF
128	IF fld_flag EQ 2 THEN fldrem_t2 = fldrem
129	sc_ampl = sqrt((moment(fldrem))[1])
130	sc_ampl2 = max(fldrem)-min(fldrem)
131	fldrem = reform(fldrem#replicate(1, long(lenght/running)), lenght)
132
133	; Compute mean error-bar using running bootstrap window (width = 10 months)
134	temperr = stat_error(fld, 10L, /mean, niter = 999)
135	err_std = sqrt( total( (temperr-mean(fld))^2 )/n_elements(temperr) )
136
137	IF NOT keyword_set(silent) THEN BEGIN
138	print, ' 1D times-serie info : mean +- err ', (moment(fld))[0], err_std
139	print, ' Mean seasonal cycle amplitude stdev / max-min', sc_ampl, sc_ampl2
140	ENDIF
141
142	; Compute std error-bar using running bootstrap window (width in plt_def)
143	IF jpt GT boot_win(0) THEN BEGIN
144	temperr = stat_error(fld, boot_win(0))
145	err_std_1 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
146	ENDIF ELSE err_std_1 = 0
147	IF jpt GT boot_win(1) THEN BEGIN
148	temperr = stat_error(fld, boot_win(1))
149	err_std_2 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
150	ENDIF ELSE err_std_2 = 0
151	IF jpt GT boot_win(2) THEN BEGIN
152	temperr = stat_error(fld, boot_win(2))
153	err_std_3 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
154	ENDIF ELSE err_std_3 = 0
155
156	IF NOT keyword_set(silent) THEN print, ' 1D times-serie anomaly info : var/stddev +- error (24,36,48)', (moment(fld - fldrem))[1], sqrt((moment(fld - fldrem))[1]), err_std_1, err_std_2, err_std_3
157	IF NOT keyword_set(silent) THEN print, ' '
158	ENDIF ELSE fldrem = mean(fld)
159
160	idx = where(fld EQ valmask)
161	fld = fld - fldrem
162	; keep mean seasonal cycle in common
163	mean_sc = fldrem
164
165	IF idx(0) NE -1 THEN BEGIN
166	fld(idx) = valmask
167	ENDIF
168
169	; compute time lag correlation
170	IF ioverchk EQ 2 AND c_correl EQ 1 THEN BEGIN
171	IF (size(fld))(1) EQ (size(fld_prev_t))(1) THEN BEGIN
172	lag_array = findgen(2*lag_correl+1)-lag_correl
173	lag_correlation = C_CORRELATE(fld, fld_prev_t, lag_array)
174	print, ' 1D time series lag correlation = ', lag_correlation
175	indx = where (lag_correlation EQ max(lag_correlation))
176	IF indx[0] NE -1 THEN BEGIN
177	print, ' Lag of max correlation =', lag_array(indx),lag_correlation(indx)
178	print, ' Lag correlation limit [-N,N] =', lag_correl
179	ENDIF
180	ENDIF
181	ENDIF
182
183	;
184	; apply running mean stdev
185	;
186	IF run_stddev GT 0 THEN BEGIN
187	print, ' Compute running std dev, window = ',run_stddev
188	print, ' '
189	full = 0
190	nyears = lenght/running
191	bootfoo = run_stddev/running
192	stat = spectra(fld, nyears, tukey, run_stddev/running, full, bootfoo)
193	stdarr = fld
194	stdarr[*] = !VALUES.F_NAN
195	size_data = (size(stat.run_std))[1]
196	stdarr[run_stddev/2:(run_stddev/2)+size_data-1]= stat.run_std
197	fld = stdarr
198	ENDIF
199
200	fld_prev_t = fld
201
202	ENDIF ELSE BEGIN ; hovmoeller
203	; build array to remove from time serie
204	IF debug_w THEN print, ' anomaly = remove average running mean (hovmoeller)'
205	CASE type OF
206	'xt': BEGIN
207	nx = (size(fld))[1]
208	lenght = (size(fld))[2]
209	fldrem = reform(fld, nx, running, lenght/running)
210	fldrem = total(fldrem, 3)/(lenght/running)
211	fldrem = fldrem[*]
212	fldrem = reform(fldrem#replicate(1, long(lenght/running)),nx, lenght)
213	END
214	'yt': BEGIN
215	ny = (size(fld))[1]
216	lenght = (size(fld))[2]
217	fldrem = reform(fld, ny, running, lenght/running)
218	fldrem = total(fldrem, 3)/(lenght/running)
219	fldrem = fldrem[*]
220	fldrem = reform(fldrem#replicate(1, long(lenght/running)),ny, lenght)
221	END
222	'zt': BEGIN
223	nz = (size(fld))[1]
224	lenght = (size(fld))[2]
225	fldrem = reform(fld, nz, running, lenght/running)
226	fldrem = total(fldrem, 3)/(lenght/running)
227	fldrem = fldrem[*]
228	fldrem = reform(fldrem#replicate(1, long(lenght/running)),nz, lenght)
229	END
230	'xyt': BEGIN
231	nx = (size(fld))[1]
232	ny = (size(fld))[2]
233	lenght = (size(fld))[3]
234	fldrem = reform(fld, nx*ny, running, lenght/running)
235	fldrem = total(fldrem, 3)/(lenght/running)
236	fldrem = fldrem[*]
237	fldrem = reform(fldrem#replicate(1, long(lenght/running)),nx*ny, lenght)
238	fldrem = reform(fldrem, nx, ny, lenght)
239	END
240	ELSE: stop, 'The type '+type+' is not handled in trends to build the seasonal cycle fldrem. Please change your cmd.plt'
241	ENDCASE
242
243	idx = where(fld EQ valmask)
244	fld = fld - fldrem
245	; keep mean seasonal cycle in common
246	mean_sc = fldrem
247	IF idx(0) NE -1 THEN fld(idx) = valmask
248	ENDELSE
249
250	END
251	'5': BEGIN ; Apply digital filter
252	; --------------------
253	IF strlen(trend_type) EQ 1 THEN BEGIN
254	print, ' *** specify filter width as <time_int>_<lower>_<upper>'
255	return, -1
256	ENDIF
257	IF (size(fld))[0] EQ 1 THEN BEGIN ; 1D time serie
258	print, ' 1D times-serie info : mean ', (moment(fld))[0]
259	print, ' 1D times-serie anomaly info : var/stddev ', (moment(fld))[1], sqrt((moment(fld))[1])
260	print, ' '
261	print, ' *** filter on 1D data not ready'
262
263	ENDIF ELSE BEGIN
264	print, ' *** filter on 2D data not ready'
265	return, -1
266	ENDELSE
267
268	; filter = digital_filter(0.15, 0.66666, 50, 15)
269	; fraction of niquisk freq (upper, lower
270	; limit indays))
271	; output_1d = convol(input_1d, filter)
272	END
273	'6': field_intl = 1; time integral done after
274	ELSE:
275	ENDCASE
276
277	; compute field time integral if required
278	IF n_elements(field_int) EQ 0 THEN field_int = 0
279	IF field_intl EQ 1 THEN BEGIN
280	; running integral ?
281	IF strlen(trend_type) GT 1 THEN $
282	int_win = long(strmid(trend_type, 1, strlen(trend_type)-1)) ELSE int_win = 0
283	CASE (size(fld))[0] OF
284	1: BEGIN
285	print, ' -> compute 1D data field time integral / lenght = ', strtrim(string(int_win), 2)
286	fldint = fld
287	lenght = (size(fld))[1]
288	IF int_win EQ 0 THEN BEGIN ; integral from start
289	FOR t = 0, lenght-1 DO BEGIN
290	fldint(t) = total(fld[0:t])
291	ENDFOR
292	ENDIF ELSE BEGIN ; running integral
293	FOR t = 0, int_win-1 DO BEGIN
294	fldint(t) = total(fld[0:t])
295	ENDFOR
296	FOR t = int_win,lenght-1 DO BEGIN
297	fldint(t) = total(fld[t-int_win+1:t])
298	ENDFOR
299	ENDELSE
300	fld = fldint
301	END
302	2: BEGIN
303	print, ' -> compute 2D data field time integral / lenght = ', strtrim(string(int_win), 2)
304	CASE type OF
305	'yt': BEGIN
306	ny = (size(fld))[1]
307	lenght = (size(fld))[2]
308	FOR i = 0, ny-1 DO BEGIN
309	zwrk = fld(i, *)
310	zint = zwrk
311	IF int_win EQ 0 THEN BEGIN ; integral from start
312	FOR t = 0, lenght-1 DO BEGIN
313	zint(t) = total(zwrk[0:t])
314	ENDFOR
315	ENDIF ELSE BEGIN ; running integral
316	FOR t = 0, int_win-1 DO BEGIN
317	zint(t) = total(zwrk[0:t])
318	ENDFOR
319	FOR t = int_win,lenght-1 DO BEGIN
320	zint(t) = total(zwrk[t-int_win+1:t])
321	ENDFOR
322	ENDELSE
323	fld(i, *) = zint
324	ENDFOR
325	END
326	'xt': BEGIN
327	nx = (size(fld))[1]
328	lenght = (size(fld))[2]
329	FOR i = 0, nx-1 DO BEGIN
330	zwrk = fld(i, *)
331	zint = zwrk
332	IF int_win EQ 0 THEN BEGIN ; integral from start
333	FOR t = 0, lenght-1 DO BEGIN
334	zint(t) = total(zwrk[0:t])
335	ENDFOR
336	ENDIF ELSE BEGIN ; running integral
337	FOR t = 0, int_win-1 DO BEGIN
338	zint(t) = total(zwrk[0:t])
339	ENDFOR
340	FOR t = int_win,lenght-1 DO BEGIN
341	zint(t) = total(zwrk[t-int_win+1:t])
342	ENDFOR
343	ENDELSE
344	fld(i, *) = zint
345	ENDFOR
346	END
347	ELSE: BEGIN
348	print, ' -> 2D data field type ', type, ' not ready'
349	END
350	ENDCASE
351	END
352	ENDCASE
353	ENDIF
354
355	; if write_data then write to ascii/cdf file
356
357	IF n_elements(write_data) EQ 0 THEN write_data = 0
358	IF write_data NE 0 THEN BEGIN
359	IF write_data GE 1 THEN BEGIN ; ascii
360	scale = 1.
361	get_lun, nuldat
362	filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'.asc'
363	openw, nuldat, asciidir+filename
364	print, ' -> writing 1D data to ', asciidir+filename
365	print, ' '
366	printf, nuldat, fld*scale, format = '(f8.4)'
367	free_lun, nuldat
368	close, nuldat
369	ENDIF ELSE BEGIN ; netCDF (tcdf case)
370	flddesc = field
371	flddesc.data = fld
372	cdf_description = nc_build(cmd_wrk, flddesc, type, vargrid)
373	fldcdf = {data:fld, units:field.units, short_name:field.legend+'_'+legbox, long_name:field.legend+' ('+legbox+')', missing_value:valmask, direc:type}
374	file_out = cmd_wrk.var+'_'+strtrim(string(FORMAT = '(I2.2)', ABS(write_data)), 2)+'.nc'
375	pltcmd ='nc_put, fldcdf, file_out, ncdf_db = hom_idl'+cdf_description
376	printf, nulhis, strcompress(pltcmd)
377	res = execute(pltcmd)
378	ENDELSE
379	ENDIF
380
381	IF debug_w THEN print, ' Vargrid : ', vargrid
382	IF debug_w THEN print, ' Leaving trends.pro'
383
384	return, fld
385
386	END
387

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format