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 | ;------------------------------------------------------------ |
---|
75 | FUNCTION 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 |
---|
499 | END |
---|
500 | |
---|
501 | |
---|
502 | |
---|
503 | |
---|
504 | |
---|
505 | |
---|