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 | |
---|