Context Navigation

trends.pro @ 27

Last change on this file since 27 was 21, checked in by kolasinski, 17 years ago
Bugfix in trends : corrections of the errors of the standard deviation of 1D time-series
File size: 14.7 KB

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

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

Download in other formats:

Original Format