Context Navigation

trends.pro @ 92

Last change on this file since 92 was 92, checked in by ericg, 16 years ago
Misc EG updates following merge with new nc_read.pro (r86)
File size: 15.4 KB

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

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

Download in other formats:

Original Format