Context Navigation

trends.pro @ 140

Last change on this file since 140 was 140, checked in by ericg, 15 years ago
Added error message in trends if jpt MOD running <> 0
File size: 15.8 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	; check that jpt is a multiple of running
118	IF NOT keyword_set(silent) THEN print, ' jpt MOD running = ', jpt, running, jpt MOD running
119
120	IF (jpt MOD running) NE 0 THEN BEGIN
121	print, ' *** ERROR: jpt (number of dates) is not a multiple of <n> in @t4<n>', jpt, running, jpt MOD running
122	stop
123	ENDIF
124
125	IF (size(fld))[0] EQ 1 THEN BEGIN ; 1D time serie
126	; build array to remove from time serie
127	IF running GE 2 THEN BEGIN
128	lenght = (size(fld))[1]
129	fldrem = fltarr(running)
130	FOR t = 0, running-1 DO BEGIN
131	fldrem(t) = mean(fld(long(findgen(lenght/running)*running+t)))
132	ENDFOR
133	IF fld_flag EQ 1 THEN BEGIN
134	fldrem_t1 = fldrem
135	ENDIF
136	IF fld_flag EQ 2 THEN fldrem_t2 = fldrem
137	sc_ampl = sqrt((moment(fldrem))[1])
138	sc_ampl2 = max(fldrem)-min(fldrem)
139	fldrem = reform(fldrem#replicate(1, long(lenght/running)), lenght)
140
141	; Compute mean error-bar using running bootstrap window (width = 10 months)
142	temperr = stat_error(fld, 10L, /mean, niter = 999)
143	err_std = sqrt( total( (temperr-mean(fld))^2 )/n_elements(temperr) )
144
145	IF NOT keyword_set(silent) THEN BEGIN
146	print, ' 1D times-serie info : mean +- err ', (moment(fld))[0], err_std
147	print, ' Mean seasonal cycle amplitude stdev / max-min', sc_ampl, sc_ampl2
148	ENDIF
149
150	; Compute std error-bar using running bootstrap window (width in plt_def)
151	IF jpt GT boot_win(0) THEN BEGIN
152	temperr = stat_error(fld, boot_win(0))
153	err_std_1 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
154	ENDIF ELSE err_std_1 = 0
155	IF jpt GT boot_win(1) THEN BEGIN
156	temperr = stat_error(fld, boot_win(1))
157	err_std_2 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
158	ENDIF ELSE err_std_2 = 0
159	IF jpt GT boot_win(2) THEN BEGIN
160	temperr = stat_error(fld, boot_win(2))
161	err_std_3 = sqrt( total((temperr-sqrt((moment(fld - fldrem))[1]) )^2)/n_elements(temperr) )
162	ENDIF ELSE err_std_3 = 0
163
164	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
165	IF NOT keyword_set(silent) THEN print, ' '
166	ENDIF ELSE fldrem = mean(fld)
167
168	idx = where(fld EQ valmask)
169	fld = fld - fldrem
170	; keep mean seasonal cycle in common
171	mean_sc = fldrem
172
173	IF idx(0) NE -1 THEN BEGIN
174	fld(idx) = valmask
175	ENDIF
176
177	; compute time lag correlation
178	IF ioverchk EQ 2 AND c_correl EQ 1 THEN BEGIN
179	IF (size(fld))(1) EQ (size(fld_prev_t))(1) THEN BEGIN
180	lag_array = findgen(2*lag_correl+1)-lag_correl
181	lag_correlation = C_CORRELATE(fld, fld_prev_t, lag_array)
182	print, ' 1D time series lag correlation = ', lag_correlation
183	indx = where (lag_correlation EQ max(lag_correlation))
184	IF indx[0] NE -1 THEN BEGIN
185	print, ' Lag of max correlation =', lag_array(indx),lag_correlation(indx)
186	print, ' Lag correlation limit [-N,N] =', lag_correl
187	ENDIF
188	ENDIF
189	ENDIF
190
191	;
192	; apply running mean stdev
193	;
194	IF run_stddev GT 0 THEN BEGIN
195	print, ' Compute running std dev, window = ',run_stddev
196	print, ' '
197	full = 0
198	nyears = lenght/running
199	bootfoo = run_stddev/running
200	stat = spectra(fld, nyears, tukey, run_stddev/running, full, bootfoo)
201	stdarr = fld
202	stdarr[*] = !VALUES.F_NAN
203	size_data = (size(stat.run_std))[1]
204	stdarr[run_stddev/2:(run_stddev/2)+size_data-1]= stat.run_std
205	fld = stdarr
206	ENDIF
207
208	fld_prev_t = fld
209
210	ENDIF ELSE BEGIN ; hovmoeller
211	; build array to remove from time serie
212	IF NOT keyword_set(silent) THEN print, ' anomaly = remove average running mean (hovmoeller)'
213	CASE type OF
214	'xt': BEGIN
215	nx = (size(fld))[1]
216	lenght = (size(fld))[2]
217	fldrem = reform(fld, nx, running, lenght/running)
218	fldrem = total(fldrem, 3)/(lenght/running)
219	fldrem = fldrem[*]
220	fldrem = reform(fldrem#replicate(1, long(lenght/running)),nx, lenght)
221	END
222	'yt': BEGIN
223	ny = (size(fld))[1]
224	lenght = (size(fld))[2]
225	fldrem = reform(fld, ny, running, lenght/running)
226	fldrem = total(fldrem, 3)/(lenght/running)
227	fldrem = fldrem[*]
228	fldrem = reform(fldrem#replicate(1, long(lenght/running)),ny, lenght)
229	END
230	'zt': BEGIN
231	nz = (size(fld))[1]
232	lenght = (size(fld))[2]
233	fldrem = reform(fld, nz, running, lenght/running)
234	fldrem = total(fldrem, 3)/(lenght/running)
235	fldrem = fldrem[*]
236	fldrem = reform(fldrem#replicate(1, long(lenght/running)),nz, lenght)
237	END
238	'xyt': BEGIN
239	nx = (size(fld))[1]
240	ny = (size(fld))[2]
241	lenght = (size(fld))[3]
242	fldrem = reform(fld, nx*ny, running, lenght/running)
243	fldrem = total(fldrem, 3)/(lenght/running)
244	fldrem = fldrem[*]
245	fldrem = reform(fldrem#replicate(1, long(lenght/running)),nx*ny, lenght)
246	fldrem = reform(fldrem, nx, ny, lenght)
247	END
248	ELSE: stop, 'The type '+type+' is not handled in trends to build the seasonal cycle fldrem. Please change your cmd.plt'
249	ENDCASE
250
251	idx = where(fld EQ valmask)
252	fld = fld - fldrem
253	; keep mean seasonal cycle in common
254	mean_sc = fldrem
255	IF idx(0) NE -1 THEN fld(idx) = valmask
256	ENDELSE
257
258	END
259	'5': BEGIN ; Apply digital filter
260	; --------------------
261	IF strlen(trend_type) EQ 1 THEN BEGIN
262	print, ' *** specify filter width as <time_int>_<lower>_<upper>'
263	return, -1
264	ENDIF
265	IF (size(fld))[0] EQ 1 THEN BEGIN ; 1D time serie
266	print, ' 1D times-serie info : mean ', (moment(fld))[0]
267	print, ' 1D times-serie anomaly info : var/stddev ', (moment(fld))[1], sqrt((moment(fld))[1])
268	print, ' '
269	print, ' *** filter on 1D data not ready'
270
271	ENDIF ELSE BEGIN
272	print, ' *** filter on 2D data not ready'
273	return, -1
274	ENDELSE
275
276	; filter = digital_filter(0.15, 0.66666, 50, 15)
277	; fraction of niquisk freq (upper, lower
278	; limit indays))
279	; output_1d = convol(input_1d, filter)
280	END
281	'6': field_intl = 1; time integral done after
282	ELSE:
283	ENDCASE
284
285	; compute field time integral if required
286	IF n_elements(field_int) EQ 0 THEN field_int = 0
287	IF field_intl EQ 1 THEN BEGIN
288	; running integral ?
289	IF strlen(trend_type) GT 1 THEN $
290	int_win = long(strmid(trend_type, 1, strlen(trend_type)-1)) ELSE int_win = 0
291	CASE (size(fld))[0] OF
292	1: BEGIN
293	print, ' -> compute 1D data field time integral / lenght = ', strtrim(string(int_win), 2)
294	fldint = fld
295	lenght = (size(fld))[1]
296	IF int_win EQ 0 THEN BEGIN ; integral from start
297	FOR t = 0, lenght-1 DO BEGIN
298	fldint(t) = total(fld[0:t])
299	ENDFOR
300	ENDIF ELSE BEGIN ; running integral
301	FOR t = 0, int_win-1 DO BEGIN
302	fldint(t) = total(fld[0:t])
303	ENDFOR
304	FOR t = int_win,lenght-1 DO BEGIN
305	fldint(t) = total(fld[t-int_win+1:t])
306	ENDFOR
307	ENDELSE
308	fld = fldint
309	END
310	2: BEGIN
311	print, ' -> compute 2D data field time integral / lenght = ', strtrim(string(int_win), 2)
312	CASE type OF
313	'yt': BEGIN
314	ny = (size(fld))[1]
315	lenght = (size(fld))[2]
316	FOR i = 0, ny-1 DO BEGIN
317	zwrk = fld(i, *)
318	zint = zwrk
319	IF int_win EQ 0 THEN BEGIN ; integral from start
320	FOR t = 0, lenght-1 DO BEGIN
321	zint(t) = total(zwrk[0:t])
322	ENDFOR
323	ENDIF ELSE BEGIN ; running integral
324	FOR t = 0, int_win-1 DO BEGIN
325	zint(t) = total(zwrk[0:t])
326	ENDFOR
327	FOR t = int_win,lenght-1 DO BEGIN
328	zint(t) = total(zwrk[t-int_win+1:t])
329	ENDFOR
330	ENDELSE
331	fld(i, *) = zint
332	ENDFOR
333	END
334	'xt': BEGIN
335	nx = (size(fld))[1]
336	lenght = (size(fld))[2]
337	FOR i = 0, nx-1 DO BEGIN
338	zwrk = fld(i, *)
339	zint = zwrk
340	IF int_win EQ 0 THEN BEGIN ; integral from start
341	FOR t = 0, lenght-1 DO BEGIN
342	zint(t) = total(zwrk[0:t])
343	ENDFOR
344	ENDIF ELSE BEGIN ; running integral
345	FOR t = 0, int_win-1 DO BEGIN
346	zint(t) = total(zwrk[0:t])
347	ENDFOR
348	FOR t = int_win,lenght-1 DO BEGIN
349	zint(t) = total(zwrk[t-int_win+1:t])
350	ENDFOR
351	ENDELSE
352	fld(i, *) = zint
353	ENDFOR
354	END
355	ELSE: BEGIN
356	print, ' -> 2D data field type ', type, ' not ready'
357	END
358	ENDCASE
359	END
360	ENDCASE
361	ENDIF
362
363	; if write_data then write to ascii/cdf file
364
365	IF n_elements(write_data) EQ 0 THEN write_data = 0
366	IF write_data NE 0 THEN BEGIN
367	IF write_data GE 1 THEN BEGIN ; ascii
368	scale = 1.
369	get_lun, nuldat
370	IF field.origin EQ 'div' THEN $
371	name_exp = strjoin(strsplit(cmd_wrk.exp, '/', /EXTRACT), '-', /SINGLE) ELSE name_exp = cmd_wrk.exp
372	filename = name_exp+'_'+ cmd_wrk.var+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'.asc'
373	openw, nuldat, asciidir+filename
374	print, ' -> writing 1D data to ', asciidir+filename
375	print, ' '
376	printf, nuldat, fld*scale, format = '(f8.4)'
377	free_lun, nuldat
378	close, nuldat
379	ENDIF ELSE BEGIN ; netCDF (tcdf case)
380	flddesc = field
381	flddesc.data = fld
382	cdf_description = nc_build(cmd_wrk, flddesc, type, vargrid)
383	fldcdf = {data:fld, units:field.units, short_name:field.legend+'_'+legbox, long_name:field.legend+' ('+legbox+')', missing_value:valmask, direc:type}
384	file_out = cmd_wrk.var+'_'+strtrim(string(FORMAT = '(I2.2)', ABS(write_data)), 2)+'.nc'
385	pltcmd ='nc_put, fldcdf, file_out, ncdf_db = hom_idl'+cdf_description
386	printf, nulhis, strcompress(pltcmd)
387	res = execute(pltcmd)
388	ENDELSE
389	ENDIF
390
391	IF NOT keyword_set(silent) THEN print, ' Vargrid : ', vargrid
392	IF NOT keyword_set(silent) THEN print, ' Leaving trends.pro'
393
394	return, fld
395
396	END
397

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

Download in other formats:

Original Format