Context Navigation

trends.pro @ 169

Last change on this file since 169 was 169, checked in by pinsard, 15 years ago
add compile_opt idl2, strictarrsubs and subsequent modifications
Property svn:keywords set to `Id`
File size: 15.9 KB

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

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

Download in other formats:

Original Format