source: trunk/LECTURE/read_ncdf_nomask.pro @ 2

Last change on this file since 2 was 2, checked in by opalod, 22 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 15.4 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;        BOITE: contient la boite sur laquelle on doit faire la
32;        lecture
33;        FILENAME: string contennant le nom du fichier
34;        IODIRECTORY: string contennant le nom du repertoire ds lequel
35;        se trouve le fichier
36;        TIMESTEP:activer pour specifier que debut et fin font
37;        reference a des indices de l''axe du temps et non pas a des
38;        dates.
39;        TOUT: activer si on veut lire le ficher sur l''ensemble du
40;        domaine sans tenir compte du sous domaine definit par boite
41;        ou lon1,lon2,lat1,lat2,prof1,prof2.
42;        NOSTRUCT: activer si on ne veut pas que read_ncdf retourne
43;        une structure mais uniquement le tableau se rapportant au
44;        champ.
45;
46; OUTPUTS:une stucture lisible par litchamp.pro ou un simple tableau
47; si /NOSTRUCT est active
48;
49; COMMON BLOCKS:common.pro
50;
51; SIDE EFFECTS:
52;
53; RESTRICTIONS:le champ doit avoir une dimension temporelle
54;
55; EXAMPLE:
56;
57; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
58;                      15/10/1999
59;-
60;------------------------------------------------------------
61;------------------------------------------------------------
62;------------------------------------------------------------
63FUNCTION read_ncdf,name,debut,fin, pour_etre_compatible, BOITE=boite, FILENAME = filename $
64             , IODIRECTORY = iodirectory, PARENTIN = parentin, TIMESTEP = timestep $
65             , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, _EXTRA = ex
66@common
67;------------------------------------------------------------
68; we find the filename.
69;------------------------------------------------------------
70; is parent a valid widget ?
71   if keyword_set(parentin) then BEGIN
72      parent = long(parentin)
73      parent = parent*widget_info(parent,/managed)
74   ENDIF
75   if n_elements(iodirectory) EQ 0 AND n_elements(iodir) NE 0 then iodirectory = iodir
76   filename = isafile(filename = filename, IODIRECTORY = iodirectory)
77;------------------------------------------------------------
78; ouverture du fichier nom
79;------------------------------------------------------------
80   if size(filename, /type) NE 7 then $
81    return, report('read_ncdf cancelled')
82   cdfid=ncdf_open(filename)
83   contient=ncdf_inquire(cdfid)
84;------------------------------------------------------------
85; we check if the variable name exists in the file.
86;------------------------------------------------------------
87   if ncdf_varid(cdfid,name) EQ -1 then $
88    return, report('variable '+name+' !C not found in the file '+filename)
89;
90   varcontient=ncdf_varinq(cdfid,name)
91;------------------------------------------------------------
92; check the time axis and the debut and fin dates
93;------------------------------------------------------------
94   if n_elements(debut) EQ 0 then begin
95      debut = 0
96      timestep = 1
97   endif
98   if keyword_set(timestep) then begin
99      premiertps = debut
100      if n_elements(fin) NE 0 then derniertps = fin ELSE derniertps = premiertps
101      jpt = derniertps-premiertps+1
102      time = lindgen(jpt)
103   ENDIF ELSE BEGIN
104      date1 = juldate(debut, /vraidate)
105      if n_elements(fin) NE 0 then date2 = juldate(fin, /vraidate) ELSE date2 = date1
106      if keyword_set(parent) then BEGIN
107         widget_control, parent, get_uvalue = top_uvalue
108         filelist = extractatt(top_uvalue,'filelist')
109         currentfile = (where(filelist EQ filename))[0]
110         time = (*(extractatt(top_uvalue, 'fileparameters'))[currentfile]).time_counter
111         premiertps = where(time EQ date1) & premiertps = premiertps[0]
112         derniertps = where(time EQ date2) & derniertps = derniertps[0]
113      ENDIF ELSE BEGIN
114; we find the infinite dimension
115         timedim = contient.recdim
116         if timedim EQ -1 then $
117          return, report('the file '+filename+' as no infinite dimension. !C Use the TIMESTEP keyword')
118; we find the FIRST time axis     
119         timeid = 0
120         repeat BEGIN           ; tant que l''on a pas trouve une variable qui n''a qu''
121                                ; une dimension: la dimension infinie
122            timecontient=ncdf_varinq(cdfid,timeid) ; que contient la variable
123            timeid = timeid+1
124         endrep until (n_elements(timecontient.dim) EQ 1 AND timecontient.dim[0] EQ contient.recdim) $
125          OR timeid EQ contient.nvars+1
126;
127         if timeid EQ contient.nvars+1 then $
128          return, report('the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword')
129; we must found the time origin of the julian calendar used in the
130; time axis.
131; does the attribut units exist for the variable time axis?
132         timeid = timeid-1
133         if timecontient.natts EQ 0 then $
134          return, report('the variable '+timecontient.name+' has no attribut.!C Use the TIMESTEP keyword or add the attribut ''units'' to the variable')
135         attnames = strarr(timecontient.natts)
136         for attiq=0,timecontient.natts-1 do attnames[attiq]=ncdf_attname(cdfid,timeid,attiq)
137         if (where(attnames EQ 'units'))[0] EQ -1 then $
138          return, report('Attribut ''units'' not found for the variable '+timecontient.name+' !C Use the TIMESTEP keyword')
139         ncdf_attget,cdfid,timeid,'units',value
140;
141; time_counter:units = "seconds since 0001-01-01 00:00:00" ;
142; time_counter:units = "hours since 0001-01-01 00:00:00" ;
143; time_counter:units = "days since 1979-01-01 00:00:00" ;
144; time_counter:units = "months since 1979-01-01 00:00:00" ;
145; time_counter:units = "years since 1979-01-01 00:00:00" ;
146;
147; we decript the "units" attribut to find the time origin
148         value = strtrim(strcompress(string(value)), 2)
149         mots = str_sep(value, ' ')
150         unite = mots[0]
151         depart = str_sep(mots[2], '-')
152         ncdf_varget, cdfid, timeid, time
153         time = long(time)
154         case strlowcase(unite) of
155            'seconds':time = julday(depart[1], depart[2], depart[0])+time/(long(24)*3600)
156            'hours':time = julday(depart[1], depart[2], depart[0])+time/long(24)
157            'days':time = julday(depart[1], depart[2], depart[0])+time
158            'months':BEGIN
159               for t = 0, n_elements(time)-1  do begin
160                  time[t] = julday(depart[1]+time[t], depart[2], depart[0])
161               endfor
162            END
163            'years':BEGIN
164               for t = 0, n_elements(time)-1  do begin
165                  time[t] = julday(depart[1], depart[2], depart[0]+time[t])
166               endfor
167            END
168            ELSE: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 ..."')
169         ENDCASE
170         premiertps = where(time GE date1) & premiertps = premiertps[0]
171         if premiertps EQ -1 THEN $
172          return, report('date 1: '+strtrim(date1, 1)+' is not found in the time axis.')
173         derniertps = where(time LE date2)
174         if derniertps[0] EQ -1 THEN $
175          return, report('the time axis as no date before date 2: '+strtrim(date2, 1))
176         derniertps = derniertps[n_elements(derniertps)-1]
177      ENDELSE
178      time = time[premiertps:derniertps]
179      jpt = derniertps-premiertps+1
180   ENDELSE
181;------------------------------------------------------------
182; nom de la grille a laquelle se rapporte le champ
183;------------------------------------------------------------
184   vargrid = strupcase(strmid(filename, strlen(filename)-9))
185   if vargrid EQ 'GRID.T.NC' OR vargrid EQ 'GRID.U.NC' $
186    OR vargrid EQ 'GRID.V.NC' OR vargrid EQ 'GRID.F.NC' $
187    OR vargrid eq 'GRID.W.NC' $
188    OR vargrid EQ 'GRID_T.NC' OR vargrid EQ 'GRID_U.NC' $
189    OR vargrid EQ 'GRID_V.NC' OR vargrid EQ 'GRID_F.NC' $
190    OR vargrid eq 'GRID_W.NC' then $
191    vargrid = strupcase(strmid(filename, strlen(filename)-4, 1)) ELSE vargrid='T'
192;---------------------------------------------------------------
193; redefinition du domaine
194;---------------------------------------------------------------
195   if keyword_set(tout) then begin
196      nx=jpi
197      ny=jpj
198      nz=jpk
199      premierx = 0
200      premiery = 0
201      premierz = 0
202      dernierx = jpi-1
203      derniery = jpj-1
204      dernierz = jpk-1
205      case strupcase(vargrid) of
206         'T':mask = tmask
207         'U':mask = umask()
208         'V':mask = vmask()
209         'W':mask = tmask
210         'F':mask = fmask()
211      endcase
212   ENDIF ELSE BEGIN
213      if keyword_set(boite) then BEGIN
214         oldboite = [lon1, lon2, lat1, lat2, prof1, prof2]
215         Case 1 Of
216            N_Elements(Boite) Eq 1:bte=[lon1, lon2, lat1, lat2, 0.,boite[0]]
217            N_Elements(Boite) Eq 2:bte=[lon1, lon2, lat1, lat2, boite[0],boite[1]]
218            N_Elements(Boite) Eq 4:bte=[Boite, prof1, prof2]
219            N_Elements(Boite) Eq 5:bte=[Boite[0:3], 0, Boite[4]]
220            N_Elements(Boite) Eq 6:bte=Boite
221            Else: return, report('Mauvaise Definition de Boite')
222         ENDCASE
223         domdef, bte,GRILLE=['T', vargrid], _extra = ex
224      ENDIF
225      grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz
226      undefine, glam & undefine, gphi & ; on libere un peu de memoire!
227   ENDELSE
228;---------------------------------------------------------------------
229; on initialise les ixmindta, iymindta au besoin
230;---------------------------------------------------------------------
231   if n_elements(jpidta) EQ 0 THEN jpidta = jpiglo
232   if n_elements(jpjdta) EQ 0 THEN jpjdta = jpjglo
233   if n_elements(jpkdta) EQ 0 THEN jpkdta = jpkglo
234   if n_elements(ixmindta) EQ 0 THEN ixmindta = 0
235   if n_elements(ixmaxdta) EQ 0 then ixmaxdta = jpidta-1
236   if ixmindta EQ -1 THEN ixmindta = 0
237   IF ixmaxdta EQ -1 then ixmaxdta = jpidta-1
238   if n_elements(iymindta) EQ 0 THEN iymindta = 0
239   IF n_elements(iymaxdta) EQ 0 then iymaxdta = jpjdta-1
240   if iymindta EQ -1 THEN iymindta = 0
241   IF iymaxdta EQ -1 then iymaxdta = jpjdta-1
242   if n_elements(izmindta) EQ 0 THEN izmindta = 0
243   IF n_elements(izmaxdta) EQ 0 then izmaxdta = jpkdta-1
244   if izmindta EQ -1 THEN izmindta = 0
245   IF izmaxdta EQ -1 then izmaxdta = jpkdta-1
246;----------------------------------------------------------------
247; on va lire le fichier
248;---------------------------------------------------------------
249   if n_elements(key_stride) LE 2 then key_stride = [1, 1, 1]
250   key_stride = 1l > long(key_stride)
251   jpitotal = long(ixmaxmesh-ixminmesh+1)
252   key_shift = long(testvar(var = key_shift))
253;
254   key = long(key_shift MOD jpitotal)
255   if key LT 0 then key = key+jpitotal
256   ixmin = ixminmesh-ixmindta
257;
258@read_ncdf_varget
259;---------------------------------------------------------------------
260; on definit les variables globales rattachees a la variable
261;---------------------------------------------------------------------
262; varname
263   varname = name
264; varunit
265   if varcontient.natts NE 0 then begin
266      attnames = strarr(varcontient.natts)
267      for attiq=0,varcontient.natts-1 do attnames[attiq]=strlowcase(ncdf_attname(cdfid,name,attiq))
268;
269      if (where(attnames EQ 'units'))[0] EQ -1 then value = '' $
270      ELSE ncdf_attget,cdfid,name,'units',value
271      varunit = strtrim(string(value), 2)
272;
273      if (where(attnames EQ 'add_offset'))[0] EQ -1 then add_offset = 0. $
274      ELSE ncdf_attget,cdfid,name,'add_offset',add_offset
275;
276      if (where(attnames EQ 'scale_factor'))[0] EQ -1 then scale_factor = 1. $
277      ELSE ncdf_attget,cdfid,name,'scale_factor',scale_factor
278;
279      if (where(attnames EQ 'missing_value'))[0] EQ -1 then missing_value = 'no' $
280      ELSE ncdf_attget,cdfid,name,'missing_value',missing_value
281;
282   ENDIF ELSE BEGIN
283      varunit = ''
284      add_offset = 0.
285      scale_factor = 1.
286      missing_value = 'no'
287   ENDELSE
288; vardate
289; on construit une belle date lisible en fonction du langage specifie !!!
290   nothing = juldate(debut,/vrai) ; bidouille pour recuperer year, month, day
291   if keyword_set(language) then BEGIN
292      if language eq 'gb' THEN $
293       vardate = string(format='(C(CMoA))',31*(month-1))+' '+strtrim(day, 1)+' '+strtrim(year, 1) $
294      ELSE vardate = strtrim(day, 1)+' '+string(format='(C(CMoA))',31*(month-1))+' '+strtrim(year, 1)
295   ENDIF ELSE vardate = strtrim(day, 1)+' '+string(format='(C(CMoA))',31*(month-1))+' '+strtrim(year, 1)
296; varexp ???
297; on applique la valeur valmask sur les points terre
298   if NOT keyword_set(cont_nofill) then begin
299      valmask = 1e20
300      case 1 of
301         varcontient.ndims eq 2:BEGIN ;xy array
302            mask = mask[*, *, 0]
303            earth = where(mask EQ 0)
304            if earth[0] NE -1 then res[earth] = 1e20
305         END
306         varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array
307            earth = where(mask EQ 0)
308            if earth[0] NE -1 then res[earth] = 1e20
309         END
310         varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array
311            mask = mask[*, *, 0]
312            earth = where(mask EQ 0)
313            if earth[0] NE -1 then BEGIN
314               earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt)
315               res[earth] = 1e20
316            END
317         END
318         varcontient.ndims eq 4:BEGIN ;xyzt array
319            earth = where(mask EQ 0)
320            if earth[0] NE -1 then BEGIN
321               earth = earth#replicate(1, jpt)+replicate(nx*ny*nz, n_elements(earth))#lindgen(jpt)
322               res[earth] = 1e20
323            END
324         END
325      endcase
326   endif
327; on applique les add_offset, scale_factor et missing_value
328   if size(missing_value, /type) NE 7 then BEGIN
329      if missing_value NE valmask then begin
330         if abs(missing_value) LT 1e6 then missing = where(res EQ missing_value) $
331         ELSE missing = where(abs(res) gt abs(missing_value)/10.)
332         if missing[0] NE -1 then res[temporary(missing)] = !values.f_nan
333      endif
334   ENDIF
335   if scale_factor NE 1 then res = res*scale_factor
336   if add_offset NE 0 then res = res+add_offset
337;
338   if keyword_set(key_reverse) then res = reverse(temporary(res),  2)
339;---------------------------------------------------------------------
340   ncdf_close,cdfid
341;---------------------------------------------------------------------
342   if keyword_set(oldboite) then domdef, oldboite,GRILLE=['T', vargrid]
343   if keyword_set(nostruct) then return, res $
344   ELSE return, {tab:res, grille:vargrid, unite:varunit, experience:varexp, nom:varname}
345
346end
347
348
349
350
351
352
Note: See TracBrowser for help on using the repository browser.