source: trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro @ 74

Last change on this file since 74 was 74, checked in by smasson, 18 years ago

debug xxx and cie + clean data file + rm perldoc_idl

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 20.7 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME: read_ncdf
6;
7; PURPOSE:fonction de lecture pour fichier net_cdf.
8; Ce programme, est moins universel que ncdf_lec (il fait appelle au
9; variables declarees dans common.pro) mais il est du cop bcp plus
10; facile d''utilisation. Il prend en compte la declaration des
11; differents zoom qui ont ete definis (ixminmesh...premierx...) la
12; declaration de la variable key_shift... bref le resultat de
13; read_ncdf peut dorectement etre utilise dans plt...
14; C''est aussi ce programme qui est utilise par defaut dans mes
15; widgets pour la partie lecture.
16;
17; CATEGORY:lecture de fichiers NetCdf
18;
19; CALLING SEQUENCE:res = read_ncdf(name,debut[,fin])
20;
21; INPUTS: name: un string definissant le champ a lire.
22;         debut et fin: sont relatifs a l''axe des temps. Ce peut etre
23;         - 2 dates du type yyyymmdd et ds ce cas on selectionne les
24;         dates qui sont comprisent entre ces 2 dates.
25;         - 2 indices qui definissent entre quel et quel pas de temps
26;           on doit extraire la dimension temporelle.
27;         exp: ne sert a rien!
28;
29; KEYWORD PARAMETERS: utilisables hors du contexte des widgets
30;
31;        BOXZOOM: contient la boxzoom sur laquelle on doit faire la
32;        lecture
33;        FILENAME: string contennant le nom du fichier
34;        /INIT; to call automatically initncdf, filename and thus
35;        redefine all the grid parameters
36;        GRID='[UTVWF]' to specify the type of grid. Defaut is (1)
37;        based on the name of the file if the file ends by
38;        GRID[._][TUVFW].NC (not case sensible) or (2) T if case (1)
39;        is not found.
40;        IODIRECTORY;a string giving the name of iodirectory (see
41;        isafile.pro for all possibilities). default value is common
42;        variable iodir
43;        TIMESTEP:activer pour specifier que debut et fin font
44;        reference a des indices de l''axe du temps et non pas a des
45;        dates.
46;        TOUT: activer si on veut lire le ficher sur l''ensemble du
47;        domaine sans tenir compte du sous domaine definit par boxzoom
48;        ou lon1,lon2,lat1,lat2,vert1,vert2.
49;        NOSTRUCT: activer si on ne veut pas que read_ncdf reourne
50;        une structure mais uniquement le tableau se rapportant au
51;        champ.
52;        TIMEVAR: a string to define the name of the variable that
53;        contains the time axis. This keyword can be usefull if there
54;        is no unlimited dimension or if the time axis selected by defaut
55;        (the first 1D array with unlimited dimension) is not the good one
56;
57;
58; OUTPUTS:une stucture lisible par litchamp.pro ou un simple tableau
59; si /NOSTRUCT est active
60;
61; COMMON BLOCKS:common.pro
62;
63; SIDE EFFECTS:
64;
65; RESTRICTIONS:le champ doit avoir une dimension temporelle
66;
67; EXAMPLE:
68;
69; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
70;                      15/10/1999
71;-
72;------------------------------------------------------------
73;------------------------------------------------------------
74;------------------------------------------------------------
75FUNCTION read_ncdf, name, debut, fin, pour_etre_compatible, BOXZOOM = boxzoom, FILENAME = filename $
76                    , PARENTIN = parentin, TIMESTEP = timestep, TIMEVAR = timevar $
77                    , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, INIT = init $
78                    , GRID = grid, FBASE2TBASE = fbase2tbase, _EXTRA = ex
79;---------------------------------------------------------
80@cm_4mesh
81@cm_4data
82@cm_4cal
83  IF NOT keyword_set(key_forgetold) THEN BEGIN
84@updatenew
85@updatekwd
86  ENDIF
87;------------------------------------------------------------
88; we find the filename.
89;------------------------------------------------------------
90;  print,filename
91; is parent a valid widget ?
92  if keyword_set(parentin) then BEGIN
93    parent = long(parentin)
94    parent = parent*widget_info(parent, /managed)
95  ENDIF
96  filename = isafile(filename = filename, IODIRECTORY = iodir, _EXTRA = ex)
97;------------------------------------------------------------
98; ouverture du fichier nom
99;------------------------------------------------------------
100  if size(filename, /type) NE 7 then $
101    return, report('read_ncdf cancelled')
102  IF !version.OS_FAMILY EQ 'unix' THEN spawn, '\file '+filename+' > /dev/null'
103  cdfid = ncdf_open(filename)
104  contient = ncdf_inquire(cdfid)
105;------------------------------------------------------------
106; we check if the variable name exists in the file.
107;------------------------------------------------------------
108  if ncdf_varid(cdfid, name) EQ -1 then BEGIN
109    ncdf_close, cdfid
110    return, report('variable '+name+' !C not found in the file '+filename)
111  ENDIF
112  varcontient = ncdf_varinq(cdfid, name)
113;------------------------------------------------------------
114; shall we redefine the grid parameters
115;------------------------------------------------------------
116  if keyword_set(init) THEN initncdf, filename, _extra = ex
117;------------------------------------------------------------
118; check the time axis and the debut and fin dates
119;------------------------------------------------------------
120  if n_elements(debut) EQ 0 then begin
121    debut = 0
122    timestep = 1
123  endif
124  if keyword_set(timestep) then begin
125    firsttps = debut[0]
126    if n_elements(fin) NE 0 then lasttps = fin[0] ELSE lasttps = firsttps
127    jpt = lasttps-firsttps+1
128    time = julday(1, 1, 1) + lindgen(jpt)
129  ENDIF ELSE BEGIN
130    if keyword_set(parent) then BEGIN
131      widget_control, parent, get_uvalue = top_uvalue
132      filelist = extractatt(top_uvalue, 'filelist')
133      IF filelist[0] EQ 'many !' THEN filelist = filename
134      currentfile = (where(filelist EQ filename))[0]
135      time = (*(extractatt(top_uvalue, 'fileparameters'))[currentfile]).time_counter
136      date1 = date2jul(debut[0])
137      if n_elements(fin) NE 0 then date2 = date2jul(fin[0]) ELSE date2 = date1
138      firsttps = where(time EQ date1) & firsttps = firsttps[0]
139      lasttps = where(time EQ date2) & lasttps = lasttps[0]
140    ENDIF ELSE BEGIN
141      IF keyword_set(timevar) THEN BEGIN
142        timeid = ncdf_varid(cdfid, timevar)
143        IF timeid EQ -1 THEN BEGIN
144          ncdf_close, cdfid
145          return, report('the file '+filename+' as no variable ' + timevar $
146                         + '. !C Use the TIMESTEP keyword')
147        endif
148        timecontient = ncdf_varinq(cdfid, timeid)
149        contient.recdim = timecontient.dim[0]
150      ENDIF ELSE BEGIN
151; we find the infinite dimension
152        timedim = contient.recdim
153        if timedim EQ -1 then BEGIN
154          ncdf_close, cdfid
155          return, report('the file '+filename+' as no infinite dimension. !C Use TIMESTEP or TIMEVAR keyword')
156        endif
157; we find the FIRST time axis     
158        timeid = 0
159        repeat BEGIN       ; tant que l''on a pas trouve une variable qui n''a qu''
160                                ; une dimension: la dimension infinie
161          timecontient = ncdf_varinq(cdfid, timeid) ; que contient la variable
162          timeid = timeid+1
163        endrep until (n_elements(timecontient.dim) EQ 1 $
164                      AND timecontient.dim[0] EQ contient.recdim) $
165          OR timeid EQ contient.nvars+1
166;
167        if timeid EQ contient.nvars+1 then BEGIN
168          ncdf_close, cdfid
169          return, report('the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword')
170        endif
171        timeid = timeid-1
172      ENDELSE
173; we must found the time origin of the julian calendar used in the
174; time axis.
175; does the attribut units an dcalendar exist for the variable time axis?
176      if timecontient.natts EQ 0 then BEGIN
177        ncdf_close, cdfid
178        return, report('the variable '+timecontient.name+' has no attribut.!C Use the TIMESTEP keyword or add the attribut ''units'' to the variable')
179      endif
180      attnames = strarr(timecontient.natts)
181      for attiq = 0, timecontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, timeid, attiq)
182      if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN
183        ncdf_close, cdfid
184        return, report('Attribut ''units'' not found for the variable '+timecontient.name+' !C Use the TIMESTEP keyword')
185      ENDIF
186;
187; now we try to find the attribut called calendar...
188; the the attribute "calendar" exists?
189; If no, we suppose that the calendar is gregorian calendar
190;
191      if (where(attnames EQ 'calendar'))[0] NE -1 then BEGIN
192        ncdf_attget, cdfid, timeid, 'calendar', value
193        value = string(value)
194        CASE value OF
195          'noleap':key_caltype = 'noleap'
196          '360d':key_caltype = '360d'
197          'greg':IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
198          ELSE:BEGIN
199;            notused = report('Unknown calendar: '+value+', we use greg calendar.')
200            key_caltype = 'greg'
201          END
202        ENDCASE
203      ENDIF ELSE BEGIN
204;        notused = report('Unknown calendar, we use '+key_caltype+' calendar.')
205        IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
206      ENDELSE
207;
208; now we take acre of units attribut
209      ncdf_attget, cdfid, timeid, 'units', value
210;
211; time_counter:units = "seconds since 0001-01-01 00:00:00" ;
212; time_counter:units = "hours since 0001-01-01 00:00:00" ;
213; time_counter:units = "days since 1979-01-01 00:00:00" ;
214; time_counter:units = "months since 1979-01-01 00:00:00" ;
215; time_counter:units = "years since 1979-01-01 00:00:00" ;
216;
217; we decript the "units" attribut to find the time origin
218      value = strtrim(strcompress(string(value)), 2)
219      mots = str_sep(value, ' ')
220      unite = mots[0]
221      depart = str_sep(mots[2], '-')
222      ncdf_varget, cdfid, timeid, time
223      time = double(time)
224      unite = strlowcase(unite)
225      IF strpos(unite, 's', strlen(unite)-1) NE -1 THEN unite = strmid(unite, 0, strlen(unite)-1)
226      IF strpos(unite, 'julian_') NE -1 THEN unite = strmid(unite, 7)
227      case unite of
228        'second':time = julday(depart[1], depart[2], depart[0])+time/86400.d
229        'hour':time = julday(depart[1], depart[2], depart[0])+time/24.d
230        'day':time = julday(depart[1], depart[2], depart[0])+time
231        'month':BEGIN
232          if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m
233            time = julday(depart[1], depart[2], depart[0])+round(time*30) $
234            ELSE for t = 0, n_elements(time)-1 DO $
235            time[t] = julday(depart[1]+time[t], depart[2], depart[0])
236        END
237        'year':BEGIN
238          if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y
239            time = julday(depart[1], depart[2], depart[0])+round(time*365) $
240            ELSE for t = 0, n_elements(time)-1 do $
241            time[t] = julday(depart[1], depart[2], depart[0]+time[t])
242        END
243        ELSE:BEGIN
244          ncdf_close, cdfid
245          return, report('The "units" attribu of the time axis must be something like: !C "seconds since 0001-01-01 ..." !C "days since 1979-01-01 ..." !C "months since 1979-01-01 ..." !C "years since 1979-01-01 ..."')
246        end
247      ENDCASE
248      date1 = date2jul(debut[0])
249      if n_elements(fin) NE 0 then date2 = date2jul(fin[0]) ELSE date2 = date1
250      time = double(time)
251      firsttps = where(time GE date1) & firsttps = firsttps[0]
252      if firsttps EQ -1 THEN BEGIN
253        ncdf_close, cdfid
254        return, report('date 1: '+strtrim(jul2date(date1), 1)+' is not found in the time axis.')
255      ENDIF
256      lasttps = where(time LE date2)
257      if lasttps[0] EQ -1 THEN BEGIN
258        ncdf_close, cdfid
259        return, report('the time axis as no date before date 2: '+strtrim(jul2date(date2), 1))
260      endif
261      lasttps = lasttps[n_elements(lasttps)-1]
262      if lasttps LT firsttps then BEGIN
263        ncdf_close, cdfid
264        return, report('the time axis as no dates between date1 and  date 2: '+strtrim(jul2date(date1), 1)+' '+strtrim(jul2date(date2), 1))
265      endif
266    ENDELSE
267    time = time[firsttps:lasttps]
268    jpt = lasttps-firsttps+1
269  ENDELSE
270;------------------------------------------------------------
271; nom de la grille a laquelle se rapporte le champ
272;------------------------------------------------------------
273  IF keyword_set(grid) THEN vargrid = strupcase(grid) ELSE BEGIN
274    vargrid = 'T'               ; default definition
275    IF finite(glamu[0]) EQ 1 THEN BEGIN
276      pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_']
277      gdtype = ['T', 'U', 'V', 'W', 'F']
278      fnametest = strupcase(filename)
279      FOR i = 0, n_elements(pattern)-1 DO BEGIN
280        FOR j = 0, n_elements(gdtype)-1 DO BEGIN
281          substr = pattern[i]+gdtype[j]
282          pos = strpos(fnametest, substr)
283          IF pos NE -1 THEN $
284             vargrid = strmid(fnametest, pos+strlen(substr)-1, 1)
285        ENDFOR
286      ENDFOR
287    ENDIF
288  ENDELSE
289;---------------------------------------------------------------
290; call the init function ?
291;---------------------------------------------------------------
292
293;---------------------------------------------------------------
294; redefinition du domaine
295;---------------------------------------------------------------
296  if keyword_set(tout) then begin
297    nx = jpi
298    ny = jpj
299    nz = jpk
300    firstx = 0
301    firsty = 0
302    firstz = 0
303    lastx = jpi-1
304    lasty = jpj-1
305    lastz = jpk-1
306    case strupcase(vargrid) of
307      'T':mask = tmask
308      'U':mask = umask()
309      'V':mask = vmask()
310      'W':mask = tmask
311      'F':mask = fmask()
312    endcase
313  ENDIF ELSE BEGIN
314    if keyword_set(boxzoom) then BEGIN
315      Case 1 Of
316        N_Elements(Boxzoom) Eq 1:bte = [lon1, lon2, lat1, lat2, 0., boxzoom[0]]
317        N_Elements(Boxzoom) Eq 2:bte = [lon1, lon2, lat1, lat2, boxzoom[0], boxzoom[1]]
318        N_Elements(Boxzoom) Eq 4:bte = [Boxzoom, vert1, vert2]
319        N_Elements(Boxzoom) Eq 5:bte = [Boxzoom[0:3], 0, Boxzoom[4]]
320        N_Elements(Boxzoom) Eq 6:bte = Boxzoom
321        Else: BEGIN
322          ncdf_close, cdfid
323          return, report('Wrong Definition of Boxzoom')
324        end
325      ENDCASE
326      savedbox = 1b
327      saveboxparam, 'boxparam4rdncdf.dat'
328      domdef, bte, GRIDTYPE = ['T', vargrid], _extra = ex
329    ENDIF
330    grille, mask, glam, gphi, gdep, nx, ny, nz, firstx, firsty, firstz, lastx, lasty, lastz
331    undefine, glam & undefine, gphi & ; on libere un peu de memoire!
332  ENDELSE
333;---------------------------------------------------------------------
334; on initialise les ixmindta, iymindta au besoin
335;---------------------------------------------------------------------
336  if n_elements(jpidta) EQ 0 THEN jpidta = jpiglo
337  if n_elements(jpjdta) EQ 0 THEN jpjdta = jpjglo
338  if n_elements(jpkdta) EQ 0 THEN jpkdta = jpkglo
339  if n_elements(ixmindta) EQ 0 THEN ixmindta = 0
340  if n_elements(ixmaxdta) EQ 0 then ixmaxdta = jpidta-1
341  if ixmindta EQ -1 THEN ixmindta = 0
342  IF ixmaxdta EQ -1 then ixmaxdta = jpidta-1
343  if n_elements(iymindta) EQ 0 THEN iymindta = 0
344  IF n_elements(iymaxdta) EQ 0 then iymaxdta = jpjdta-1
345  if iymindta EQ -1 THEN iymindta = 0
346  IF iymaxdta EQ -1 then iymaxdta = jpjdta-1
347  if n_elements(izmindta) EQ 0 THEN izmindta = 0
348  IF n_elements(izmaxdta) EQ 0 then izmaxdta = jpkdta-1
349  if izmindta EQ -1 THEN izmindta = 0
350  IF izmaxdta EQ -1 then izmaxdta = jpkdta-1
351;----------------------------------------------------------------
352; on va lire le fichier
353;---------------------------------------------------------------
354  if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1]
355  key_stride = 1l > long(key_stride)
356  key_shift = long(testvar(var = key_shift))
357;
358  IF n_elements(key_yreverse) EQ 0 THEN key_yreverse = 0
359  IF keyword_set(key_yreverse) THEN BEGIN
360    tmp = jpj-1-firsty
361    firsty = jpj-1-lasty
362    lasty = tmp
363  ENDIF
364;
365  IF keyword_set(fbase2tbase) THEN BEGIN
366    case strupcase(vargrid) of
367      'U':BEGIN
368        IF NOT keyword_set(key_periodic) THEN BEGIN
369          firstx = firstx+1
370          lastx = lastx+1
371        ENDIF
372      END
373      'V':BEGIN
374        firsty = firsty+1
375        lasty = lasty+1
376      END
377      'F':BEGIN
378        firsty = firsty+1
379        lasty = lasty+1
380        IF NOT keyword_set(key_periodic) THEN BEGIN
381          firstx = firstx+1
382          lastx = lastx+1
383        ENDIF
384      END
385      ELSE:
386    endcase
387  ENDIF
388;
389  IF keyword_set(fbase2tbase) AND keyword_set(key_periodic) $
390    AND (strupcase(vargrid) EQ 'U' OR strupcase(vargrid) EQ 'F') THEN key_shift =  key_shift-1
391;
392;---------------------------------------------------------------------
393;---------------------------------------------------------------------
394@read_ncdf_varget
395;---------------------------------------------------------------------
396;---------------------------------------------------------------------
397;
398  IF keyword_set(fbase2tbase) AND keyword_set(key_periodic) $
399    AND (strupcase(vargrid) EQ 'U' OR strupcase(vargrid) EQ 'F') THEN key_shift =  key_shift+1
400;---------------------------------------------------------------------
401; on definit les variables globales rattachees a la variable
402;---------------------------------------------------------------------
403; varname
404  varname = name
405; varunit
406  if varcontient.natts NE 0 then begin
407    attnames = strarr(varcontient.natts)
408    for attiq = 0, varcontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, name, attiq)
409    lowattnames = strlowcase(attnames)
410;
411    found = (where(lowattnames EQ 'units'))[0]
412    IF found NE -1 then ncdf_attget, cdfid, name, attnames[found], value ELSE value = ''
413    varunit = strtrim(string(value), 2)
414;
415    found = (where(lowattnames EQ 'add_offset'))[0]
416    if found NE -1 then ncdf_attget, cdfid, name, attnames[found], add_offset ELSE add_offset = 0.
417;
418    found = (where(lowattnames EQ 'scale_factor'))[0]
419    if found NE -1 then ncdf_attget, cdfid, name, attnames[found], scale_factor ELSE scale_factor = 1.
420;
421    missing_value = 'no'
422    found = (where(lowattnames EQ '_fillvalue'))[0]
423    if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value
424    found = (where(lowattnames EQ 'missing_value'))[0]
425    if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value
426;
427  ENDIF ELSE BEGIN
428    varunit = ''
429    add_offset = 0.
430    scale_factor = 1.
431    missing_value = 'no'
432  ENDELSE
433; vardate
434; on construit une belle date lisible en fonction du langage specifie !!!
435  year = long(debut[0])/10000
436  month = (long(debut[0])/100) MOD 100
437  day = (long(debut[0]) MOD 100)
438  vardate = string(format = '(C(CMoA))', 31*(month-1))+' '+strtrim(day, 1)+', '+strtrim(year, 1)
439  varexp = file_basename(filename)
440
441; we apply reverse
442  if keyword_set(key_yreverse) then res = reverse(temporary(res),  2)
443  if keyword_set(key_zreverse) AND (size(res))[0] EQ 3 AND jpt EQ 1 then res = reverse(temporary(res), 3)
444  if keyword_set(key_zreverse) AND (size(res))[0] EQ 4 THEN res = reverse(temporary(res), 3)
445
446; on applique la valeur valmask sur les points terre
447  if NOT keyword_set(cont_nofill) then begin
448    valmask = 1e20
449    case 1 of
450      varcontient.ndims eq 2:BEGIN ;xy array
451        mask = mask[*, *, 0]
452        earth = where(mask EQ 0)
453      END
454      varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array
455        earth = where(mask EQ 0)
456      END
457      varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array
458        mask = mask[*, *, 0]
459        earth = where(mask EQ 0)
460        if earth[0] NE -1 then BEGIN
461          earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt)
462        END
463      END
464      varcontient.ndims eq 4:BEGIN ;xyzt array
465        earth = where(mask EQ 0)
466        if earth[0] NE -1 then BEGIN
467          earth = earth#replicate(1, jpt)+replicate(nx*ny*nz, n_elements(earth))#lindgen(jpt)
468        END
469      END
470    endcase
471  ENDIF ELSE earth = -1
472; we look for  missing_value
473  IF size(missing_value, /type) NE 7 then BEGIN
474    IF size(missing_value, /type) EQ 1 THEN BEGIN
475      IF isnumber(string(missing_value), tmp) EQ 1 THEN missing_value = tmp
476    ENDIF
477;    if missing_value NE valmask then begin
478    if abs(missing_value) LT 1e6 then missing = where(res EQ missing_value) $
479    ELSE missing = where(abs(res) gt abs(missing_value)/10.)
480;    ENDIF ELSE missing = -1
481  ENDIF ELSE missing = -1
482; on applique les add_offset, scale_factor et missing_value
483  if scale_factor NE 1 then res = temporary(res)*scale_factor
484  if add_offset NE 0 then res = temporary(res)+add_offset
485  if missing[0] NE -1 then res[temporary(missing)] = !values.f_nan
486  if earth[0] NE -1 then res[temporary(earth)] = 1e20
487;---------------------------------------------------------------------
488  ncdf_close, cdfid
489;---------------------------------------------------------------------
490  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4rdncdf.dat'
491  if keyword_set(nostruct) then return, res $
492  ELSE BEGIN
493    IF keyword_set(key_forgetold) THEN BEGIN
494      return, {arr:res, grid:vargrid, unit:varunit, experiment:varexp, name:varname}
495    ENDIF ELSE BEGIN
496      return, {tab:res, grille:vargrid, unite:varunit, experience:varexp, nom:varname}
497    ENDELSE
498  ENDELSE
499END
500
501
502
503
504
505
Note: See TracBrowser for help on using the repository browser.