source: trunk/procs/data_read.pro @ 26

Last change on this file since 26 was 26, checked in by kolasinski, 16 years ago

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