Context Navigation

data_read.pro @ 41

Last change on this file since 41 was 41, checked in by ericg, 16 years ago
Misc modifs, incl. read_from_grid in yfx, debug_w prints, removal of mesh_lmdz and small bug corrections
File size: 19.1 KB

Line
1	;---------------------------
2	; Reading Data
3	;---------------------------
4	FUNCTION data_read, cmd, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex
5
6	@common
7	@com_eg
8
9
10	; same data as previous read ?
11
12	IF debug_w THEN print, ' '
13	IF debug_w THEN print, ' ENTER data_read...'
14	IF debug_w THEN print, ' cmd at top of data_read = ', cmd
15
16	read = 1
17
18	IF cmd.var EQ cmd_prev.var AND cmd.grid EQ cmd_prev.grid AND cmd.plt EQ cmd_prev.plt AND cmd.timave EQ cmd_prev.timave AND cmd.date1 EQ cmd_prev.date1 THEN BEGIN
19	IF hotyp NE '-' THEN BEGIN
20	; time serie
21	IF cmd.spec EQ cmd_prev.spec THEN BEGIN
22	fldr = field
23	read = 0
24	ENDIF ELSE read = 1
25
26	ENDIF ELSE BEGIN
27	read = 1
28	IF strpos(cmd.plt, '#') NE -1 THEN BEGIN
29	IF cmd.plt EQ cmd_prev.plt THEN BEGIN ; not same zoom
30	fldr = field
31	read = 0
32	ENDIF ELSE begin
33	read = 1
34	ENDELSE
35	ENDIF
36	ENDELSE
37	ENDIF ELSE read = 1
38
39	read = 1
40
41	IF read EQ 1 THEN BEGIN
42
43	cmd_prev = cmd
44
45	; difference from cmd.exp ? division from exp ?
46
47	diff_from_exp = 0
48	div_from_exp = 0
49
50	; argvar = str_sep(cmd.exp,'-')
51
52	IF strpos(cmd.exp, '-') NE -1 THEN BEGIN
53	argvar = strsplit(cmd.exp,'-', /EXTRACT)
54	ENDIF ELSE IF strpos(cmd.exp, '/') NE -1 THEN BEGIN
55	argvar = strsplit(cmd.exp,'/', /EXTRACT)
56	ENDIF ELSE argvar = cmd.exp
57
58	IF n_elements(argvar) EQ 2 THEN BEGIN
59	IF strpos(cmd.exp, '-') NE -1 THEN BEGIN
60	diff_from_exp = 1
61	exp_init = cmd.exp
62	cmd.exp = argvar[0]
63	ENDIF
64	IF strpos(cmd.exp, '/') NE -1 THEN BEGIN
65	div_from_exp = 1
66	exp_init = cmd.exp
67	cmd.exp = argvar[0]
68	ENDIF
69	ENDIF
70	varexp = cmd.exp
71
72	; define data base
73
74	ncdf_db = def_dbase(cmd.exp)
75
76	IF debug_w THEN print, ' varexp defined in data_read = ', varexp
77
78	; define file name
79
80	def_file_name, cmd, ncdf_db, file_name, delta_t1
81
82	; get file if remote
83
84	get_file, file_name, ncdf_db
85
86	; read grid from file if require (@<grid> option)
87
88	IF read_grid_from_file EQ 1 THEN BEGIN
89	mesh_from_file, cmd.grid, file_name, ncdf_db, cmd.var
90	key_shift_map = key_shift
91	file_name2 = file_name
92	ncdf_db2 = ncdf_db
93	ENDIF
94
95	; define data domain
96
97	IF hotyp NE '-' THEN BEGIN ; time serie
98	; get hovmoeller box + stride
99	box_plot = def_box(cmd.plt, dimplot, legbox, time_stride)
100	IF debug_w THEN print, ' box_plot for domdef in data_read = ', box_plot
101	CASE n_elements(box_plot) OF
102	4: domdef, box_plot[0:3], /MEMEINDICES
103	6: domdef, box_plot[0:5], /MEMEINDICES
104	ENDCASE
105	; time grid for hovmoeller
106	time1 = 1+delta_t1
107	IF cmd.out EQ 'cdf' THEN nb_cycles = 1
108	timearr = compute_time(cmd.timave, cmd.date1, cmd.spec)
109	jpt = timearr.count
110	time2 = jpt+delta_t1
111	time = timearr.scale
112	niveau = 1
113	; print, ' Data_read, hotyp/jpt = ', hotyp, jpt
114	; print, ' Data_read, box = ', box_plot[0:3]
115	ENDIF ELSE BEGIN ; single date
116	jpt = 1
117	CASE strmid(cmd.timave, strlen(cmd.timave)-1, 1) OF
118	'y': time1 = delta_t1+1
119	'm': BEGIN
120	CASE strmid(cmd.timave, strlen(cmd.timave)-2, 1) OF
121	'm': time1 = delta_t1+long(strmid(cmd.date1, 0, 2)) ; mean month
122	ELSE: time1 = delta_t1+1
123	ENDCASE
124	END
125	'd': time1 = delta_t1+1
126	ELSE: time1 = delta_t1+1
127	ENDCASE
128	time2 = time1
129	time = time1
130	print, ' Setting domdef in data_read, plt typ =', plttyp
131	CASE plttyp OF
132	'plt': BEGIN
133	box_plot = def_box(cmd.plt, 2, legbox, time_stride)
134	IF strmid(cmd.var, 0,2) EQ 'vo' THEN BEGIN ; get right zindex for vertical
135	CASE vert_type OF
136	'level':domdef,[ box_plot[0:3],vert_mean],/zindex, /MEMEINDICES
137	'z':domdef,[ box_plot[0:3],vert_mean], /MEMEINDICES
138	ELSE: domdef, box_plot[0:3], /MEMEINDICES
139	ENDCASE
140	ENDIF ELSE BEGIN
141	domdef, box_plot[0:3], /MEMEINDICES
142	ENDELSE
143	END
144	'pltz': BEGIN
145	box_plot = def_box(cmd.plt, 2, legbox, time_stride)
146	domdef, box_plot[0:3], /MEMEINDICES
147	END
148	'plt1d': BEGIN
149	box_plot = def_box(cmd.plt, 2, legbox, time_stride)
150	domdef, box_plot[0:3], /MEMEINDICES
151	; force all_data read
152	all_data = 1
153	END
154	ELSE: domdef, /MEMEINDICES
155	ENDCASE
156	ENDELSE
157
158	; print, jpt, time1, time2, hotyp
159
160	; save time in common fld_att
161
162	time1_r = time1
163	time2_r = time2
164
165	; call specific read routine
166
167	CASE STRMID(cmd.var, 0, 2) OF
168	'@@': BEGIN
169	IF debug_w THEN print, 'keyword_set(all_data) : ', keyword_set(all_data)
170	fldr = macro_read(file_name, cmd.var, ncdf_db, TIME_1 = time1, TIME_2 = time2, ALL_DATA = all_data, _extra = ex)
171	IF stddev_diff EQ 1 THEN fldr.origin = 'diff'
172	END
173	ELSE: BEGIN
174	CASE plttyp OF
175	'yfx': BEGIN ; scatter plot y=f(x)
176	idx = strpos(cmd.var, '=f(')
177	var1 = strmid(cmd.var, 0, idx)
178	fldr1 = nc_read(file_name, var1, ncdf_db, $
179	TIME_1 = time1, TIME_2 = time2, ALL_DATA = all_data, _extra = ex)
180	file_name1 = file_name
181	ncdf_db1 = ncdf_db
182	CASE cmd.grid OF
183	'U': vargrid1 = 'U'
184	'V': vargrid1 = 'V'
185	ELSE: vargrid1 = 'T'
186	ENDCASE
187	IF strpos(cmd.var, '=f(next)') EQ -1 THEN BEGIN ; type 1 y=f(x) on 1 line
188	var2 = strmid(cmd.var, idx+3, strlen(cmd.var)-idx-4)
189	fldr2 = nc_read(file_name, var2, ncdf_db, $
190	TIME_1 = time1, TIME_2 = time2, ALL_DATA = all_data, _extra = ex)
191	datyp2 = datyp
192	cmd2 = cmd
193	sw_diffg = 0
194	vargrid2 = vargrid1
195	ENDIF ELSE BEGIN ; type 2 y=f(next) on 2 lines
196	cmd2 = decode_cmd(cmdline_main, idx_main+1)
197	cmd2.timave = cmd.timave
198	cmd2.date1 = cmd.date1
199	cmd2.spec = cmd.spec
200	datyp2 = def_dptyp(cmd2)
201	cmd2_back = cmd2
202	; other grid ?
203	grid2_test = cmd2.grid
204	sw_diffg = 0
205	IF grid2_test NE cmd.grid THEN BEGIN
206	sw_diffg = 1
207	jptb = jpt
208	IF debug_w THEN print, ' (in data_read) cmd2.grid <> cmd.grid ', cmd2.grid, cmd.grid
209	IF debug_w THEN print, ' exec def_grid, cmd2_back in data_read : '
210	grid2_full = cmd2_back.grid
211	def_grid, cmd2_back
212	jpt = jptb
213	ENDIF
214	fldr2 = data_read(cmd2, datyp2.hotyp, datyp2.plttyp, '1', '0', ALL_DATA = all_data, _extra = ex)
215	var2 = cmd2.var
216	print, ' Variable 2 ', var2, ' read on vargrid ', vargrid
217	jptb = jpt
218	if debug_w THEN print, ' yfx grids different ? / sw_diffg = ', sw_diffg
219	IF sw_diffg EQ 1 THEN BEGIN
220	IF debug_w THEN print, ' exec def_grid, cmd1_back in data_read : '
221	grid1_full = cmd1_back.grid
222	def_grid, cmd1_back
223	IF read_grid_from_file EQ 1 THEN BEGIN
224	mesh_from_file, cmd1_back.grid, file_name1, ncdf_db1, cmd1_back.var
225	key_shift_map = key_shift
226	ENDIF
227	ENDIF
228	jpt = jptb
229	CASE cmd2.grid OF
230	'U': vargrid2 = 'U'
231	'V': vargrid2 = 'V'
232	ELSE: vargrid2 = 'T'
233	ENDCASE
234	ENDELSE
235	; consistency check
236	IF fldr1.dim NE fldr2.dim THEN BEGIN
237	print, ' *** dimension inconsistency in y=f(x) : ', cmd.var
238	return, -1
239	ENDIF
240	; custom stucture for scatter plot
241	fldr = {name: var1, data: fldr1.data, legend: fldr1.legend, $
242	units: fldr1.units, origin: fldr1.origin, dim: fldr1.dim, $
243	name2: var2, data2: fldr2.data, legend2: fldr2.legend, $
244	units2: fldr2.units, origin2: fldr2.origin, dim2: fldr2.dim}
245	; diff not possible yet
246	IF diff_from_exp EQ 1 OR STRMID(cmd.spec, 0, 2) EQ 'd:' THEN BEGIN
247	diff_from_exp = 0
248	cmd.spec = '-'
249	print, ' data_read : difference in scatter plot not ready'
250	ENDIF
251	; re-organise cmd.var
252	cmd.var = var1
253	cmd.var2 = var2
254
255	END
256	ELSE: BEGIN
257	IF vecplot EQ 1 THEN BEGIN ; vectors case
258	idx = strpos(cmd.var, ',')
259	var1 = strmid(cmd.var, 1, idx-1)
260	var2 = strmid(cmd.var, idx+1, strlen(cmd.var)-idx-2)
261	; print, var1, var2
262	;old version file_nam = strmid(file_name, 0, strlen(file_name)-4)
263	file_nam = strmid(base_file_name+base_suffix, 0, strlen(base_file_name+base_suffix)-1)
264	vargrid = strmid(cmd.grid, 0, 1)
265	CASE STRMID(var1, 0, 2) OF
266	'@@': fldr1 = macro_read(file_name, var1, ncdf_db, TIME_1 = time1, TIME_2 = time2, _extra = ex)
267	ELSE: fldr1 = nc_read(file_nam+strmid(cmd.grid, 0, 1)+suff_domain+'.nc', var1, ncdf_db, $
268	TIME_1 = time1, TIME_2 = time2, ALL_DATA = all_data, _extra = ex)
269	ENDCASE
270	vargrid = strmid(cmd.grid, 1, 1)
271	CASE STRMID(var2, 0, 2) OF
272	'@@': fldr2 = macro_read(file_name, var2, ncdf_db, TIME_1 = time1, TIME_2 = time2, _extra = ex)
273	ELSE: fldr2 = nc_read(file_nam+strmid(cmd.grid, 1, 1)+suff_domain+'.nc', var2, ncdf_db, $
274	TIME_1 = time1, TIME_2 = time2, ALL_DATA = all_data, _extra = ex)
275	ENDCASE
276	; consistency check
277	IF fldr1.dim NE fldr2.dim THEN BEGIN
278	print, ' *** dimension inconsistency in vector plot : ', cmd.var
279	return, -1
280	ENDIF
281	; custom stucture for vector plot
282	fldr = {name: var1, data: fldr1.data, legend: fldr1.legend, $
283	units: fldr1.units, origin: fldr1.origin, dim: fldr1.dim, $
284	name2: var2, data2: fldr2.data, legend2: fldr2.legend, $
285	units2: fldr2.units, origin2: fldr2.origin, dim2: fldr2.dim}
286	IF diff_from_exp EQ 1 OR STRMID(cmd.spec, 0, 2) EQ 'd:' THEN BEGIN
287	diff_from_exp = 0
288	cmd.spec = '-'
289	print, ' data_read : difference in vector plot not ready'
290	ENDIF
291	; re-organise cmd.var
292	cmd.var = var1
293	; cmd.grid = strmid(cmd.grid, 2, 1)
294	cmd.var2 = var2
295
296
297	ENDIF ELSE BEGIN ; general case
298	fldr = nc_read(file_name, cmd.var, ncdf_db, $
299	TIME_1 = time1, TIME_2 = time2, ALL_DATA = all_data, _extra = ex)
300
301	; perform mean if hovmoeller diagram is required
302	; It is absolutely necessary for the division case with
303	; a hovmoeller diagram (before performing the division)
304	IF plttyp EQ 'pltt' AND div_from_exp EQ 1 THEN BEGIN
305	mask_z, fldr.data, cmd, boite_pltt, dimplot, legz
306	z2d = checkfield(fldr.data, plttyp, TYPE = hotyp, BOXZOOM = boite_pltt, DIREC = direc, _extra = ex)
307	fldr = { name : fldr.name, data : z2d, legend : fldr.legend, units : fldr.units, $
308	origin : fldr.origin, dim : size(z2d, /N_DIMENSIONS), direc : hotyp}
309	print, 'Averaging made before the call to SAXO plt routines !!!!'
310	ENDIF
311	; fix
312	; idx = where(fldr.data EQ 1.00000e+20)
313	; IF idx(0) NE -1 THEN fldr.data(idx) = 0.
314	ENDELSE
315	END
316	ENDCASE
317	END
318	ENDCASE
319
320	; if tkeavt, take log
321
322	; CASE cmd.var OF
323	; 'votkeavt': BEGIN
324	; idx = where(fldr.data LE valmask/10)
325	; idxm = where(fldr.data EQ valmask)
326	; fldr.data(idx) = alog10(fldr.data(idx))
327	; ; on T grid
328	; vargrid = 'T'
329	; fldr.data = (fldr.data+shift(fldr.data, 0, 0, 1))*.5
330	; fldr.data(idxm) = valmask
331	; END
332	; ELSE:
333	; ENDCASE
334
335
336	; density projection
337	IF splot EQ 1 THEN BEGIN
338	IF cmd.timave EQ '1y' AND 1 THEN BEGIN
339	print, 'data_read: Constructing annual mean from monthly means...'
340	cmd.timave = '1m'
341	year_s = cmd.date1
342	cumulative = 0.
343	FOR i = 1, 12 DO BEGIN
344	IF i EQ 1 THEN cmd.date1 = year_s + '01'
345	IF i EQ 2 THEN cmd.date1 = year_s + '02'
346	IF i EQ 3 THEN cmd.date1 = year_s + '03'
347	IF i EQ 4 THEN cmd.date1 = year_s + '04'
348	IF i EQ 5 THEN cmd.date1 = year_s + '05'
349	IF i EQ 6 THEN cmd.date1 = year_s + '06'
350	IF i EQ 7 THEN cmd.date1 = year_s + '07'
351	IF i EQ 8 THEN cmd.date1 = year_s + '08'
352	IF i EQ 9 THEN cmd.date1 = year_s + '09'
353	IF i EQ 10 THEN cmd.date1 = year_s + '10'
354	IF i EQ 11 THEN cmd.date1 = year_s + '11'
355	IF i EQ 12 THEN cmd.date1 = year_s + '12'
356	print, ' Date1 = ', cmd.date1
357	sfild = fldr
358	bin_sigma, cmd, sfild ;;;;;;;;;;;;;; prob
359	cumulative = cumulative + sfild.data
360	ENDFOR
361	fldr = sfild
362	fldr.data = cumulative/12.
363	cmd.timave = '1y'
364	cmd.date1 = year_s
365	ENDIF ELSE BEGIN
366	IF cmd.timave EQ '1m' AND strmid(cmd.plt, 0, 2) EQ 'st' AND really_1m_st EQ 1 THEN BEGIN
367	sfild = fldr
368	timedim = (size(sfild.data))[(size(sfild.data))[0]]
369	fldrd = fltarr((size(sfild.data))[1],(size(sfild.data))[2], (sig_max - sig_min)/sig_del + 1, timedim)
370	FOR i = 0, timedim-1 DO BEGIN
371	print, 'data_read: Performing monthly bining. Indices: first, last, current = 0,', timedim-1, i
372	sfild2 = {name: sfild.name, data: sfild.data(, , *, i), legend: sfild.legend, units: sfild.units, origin: sfild.origin, dim: sfild.dim-1}
373	bin_sigma, cmd, sfild2
374	fldrd(, , *, i)= sfild2.data
375	ENDFOR
376	fldr = {name: sfild.name, data:fldrd, legend: sfild.legend, units: sfild.units, origin: sfild.origin, dim: sfild.dim}
377	ENDIF ELSE BEGIN
378	sfild = fldr
379	bin_sigma, cmd, sfild
380	fldr = sfild
381	ENDELSE
382	ENDELSE
383	ENDIF
384
385	; modify data info if needed (actual data modification done in; trends.pro called by pltt.pro)
386
387	CASE plttyp OF
388	'pltt': BEGIN
389	IF run_stddev EQ 0 THEN BEGIN
390	CASE strmid(cmd.trend, 0, 1) OF
391	'1': fldr.origin = 'diff'
392	'2': fldr.origin = 'diff'
393	'3': fldr.origin = 'diff'
394	'4': fldr.origin = 'diff'
395	ELSE:
396	ENDCASE
397	ENDIF
398	END
399	ELSE:
400	ENDCASE
401
402	IF n_elements(fldr.data) EQ 1 THEN return, fldr
403
404	; Difference 2 cases : from cmd.exp or from cmd.spec
405
406
407	IF diff_from_exp EQ 1 THEN BEGIN
408
409	cmd.exp = argvar[1]
410	print, ' remove : ', cmd.exp
411
412	; read field2
413	field2 = data_read(cmd, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex)
414	IF n_elements(field2.data) EQ 1 THEN return, field2
415
416	; perform difference
417
418	diff = fldr.data-field2.data
419	IF (where(fldr.data EQ valmask))[0] NE -1 THEN $
420	diff[where(fldr.data EQ valmask)] = valmask
421	fldr.data = diff
422	fldr.origin = 'diff'
423	cmd.exp = exp_init
424	varexp = cmd.exp
425
426	ENDIF ELSE IF div_from_exp EQ 1 THEN BEGIN ; Division : 1 case
427
428	cmd.exp = argvar[1]
429	print, ' divide by : ', cmd.exp
430
431	div = fldr.data
432
433	; read field2
434
435	field2 = data_read(cmd, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex)
436	IF n_elements(field2.data) EQ 1 THEN return, field2
437
438	IF plttyp EQ 'pltt' THEN BEGIN
439	mask_z, field2.data, cmd, boite_pltt, dimplot, legz
440	z2d = checkfield(field2.data, plttyp, TYPE = hotyp, BOXZOOM = boite_pltt, DIREC = direc, _extra = ex)
441	field2 = { name : field2.name, data : z2d, legend : field2.legend, units : field2.units, $
442	origin : field2.origin, dim : size(z2d, /N_DIMENSIONS), direc : hotyp}
443	print, 'Averaging made before the call to SAXO plt routines !!!!'
444	ENDIF
445
446	; perform division
447	IF (where(field2.data EQ 0.0))[0] NE -1 THEN BEGIN
448	print, 'Be careful some data are null. They will be masked ! '
449	idx0 = where(field2.data EQ 0.0)
450	field2.data(idx0) = valmask
451	fldr.data(idx0) = valmask
452	div(idx0) = valmask
453	ENDIF
454
455	idx = where(field2.data NE valmask)
456	div(idx) = fldr.data(idx) / field2.data(idx)
457
458	fldr.data = div
459	fldr.origin = 'div'
460	cmd.exp = exp_init
461	varexp = cmd.exp
462
463	ENDIF ELSE BEGIN
464
465	; If difference, decode second field, read and perform difference
466
467	IF STRMID(cmd.spec, 0, 2) EQ 'd:' THEN BEGIN
468
469	cmd2 = cmd
470	cmdl = STRMID(cmd.spec, 2, strlen(cmd.spec)-2)
471
472	print, ' remove : ', cmdl
473
474	; decode field2
475
476	argvar = str_sep(cmdl, '/')
477
478	cmd2.exp = argvar[0]
479	cmd2.timave = argvar[1]
480	cmd2.date1 = argvar[2]
481	cmd2.spec = '-'
482	IF n_elements(argvar) EQ 4 THEN cmd2.spec = argvar[3]
483
484	; read field2
485
486	field2 = data_read(cmd2, hotyp, plttyp, dimplot, iover, ALL_DATA = all_data, _extra = ex)
487	IF n_elements(field2.data) EQ 1 THEN return, field2
488
489	; perform difference
490
491	diff = fldr.data-field2.data
492	IF (where(fldr.data EQ valmask))[0] NE -1 THEN $
493	diff[where(fldr.data EQ valmask)] = valmask
494	fldr.data = diff
495	fldr.origin = 'diff'
496	cmd.exp = cmd.exp+' - '+cmd2.exp+ ' ('+cmd2.timave+' '+cmd2.date1+')'
497	varexp = cmd.exp
498
499	ENDIF
500
501	ENDELSE
502
503	ENDIF ELSE BEGIN
504
505	print, ''
506	print, ' Data already read (previous field)'
507	print, ''
508
509	ENDELSE
510
511	IF debug_w THEN print, ' varexp before exit in data_read = ', varexp
512	IF debug_w THEN print, ' ...EXIT data_read'
513
514	return, fldr
515
516	END

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

Download in other formats:

Original Format