Context Navigation

trends.pro @ 13

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

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

Download in other formats:

Original Format