source: trunk/procs/data_read.pro @ 9

Last change on this file since 9 was 2, checked in by post_it, 17 years ago

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