1 | ;------------------------------------------------------------ |
---|
2 | ;------------------------------------------------------------ |
---|
3 | ;------------------------------------------------------------ |
---|
4 | ;+ |
---|
5 | ; NAME:initncdf |
---|
6 | ; |
---|
7 | ; PURPOSE:initfile for Netcdf file. define all the grid parameters |
---|
8 | ; |
---|
9 | ; CATEGORY: |
---|
10 | ; |
---|
11 | ; CALLING SEQUENCE:initncdf, ncfilename |
---|
12 | ; |
---|
13 | ; INPUTS:ncfilename: a string giving the name of the NetCdf file |
---|
14 | ; |
---|
15 | ; KEYWORD PARAMETERS: |
---|
16 | ; |
---|
17 | ; GLAMBOUNDARY:a 2 elements vector, {lon1,lon2], the longitute |
---|
18 | ; boundaries that should be used to visualize the data. |
---|
19 | ; lon2 > lon1 |
---|
20 | ; lon2 - lon1 le 360 |
---|
21 | ; key_shift will be defined according to GLAMBOUNDARY. |
---|
22 | ; |
---|
23 | ; /INVMASK: to inverse the mask: mask = 1-mask |
---|
24 | ; |
---|
25 | ; IODIRECTORY;a string giving the name of iodirectory (see |
---|
26 | ; isafile.pro for all possibilities). default value is common |
---|
27 | ; variable iodir |
---|
28 | ; |
---|
29 | ; MASKNAME: a string giving the name of the variable in the file |
---|
30 | ; that contains the land/sea mask |
---|
31 | ; |
---|
32 | ; MISSING_VALUE: to define (or redifine if the attribute is |
---|
33 | ; already existing) the missing values used with USEASMASK |
---|
34 | ; keyword |
---|
35 | ; |
---|
36 | ; /start1: index the axis from 1 instead of 0 when using |
---|
37 | ; /xyindex and/or zindex |
---|
38 | ; |
---|
39 | ; USEASMASK: a string giving the name of the variable in the file |
---|
40 | ; that will be used to build the land/sea mask. In this case the |
---|
41 | ; mask is based on the first record (if record dimension |
---|
42 | ; exists). The mask is build according to : |
---|
43 | ; 1 the keyword missing_value if existing |
---|
44 | ; 2 the attribute 'missing_value' if existing |
---|
45 | ; 3 NaN values if existing |
---|
46 | ; |
---|
47 | ; [XYZ]AXISNAME= a string giving the name of the variable in the file |
---|
48 | ; that contains the [xyz]axis. |
---|
49 | ; for X axis, default name must be 'x', 'longitude', 'nav_lon' or 'lon'. |
---|
50 | ; for Y axis, default name must be 'y', 'latitude', 'nav_lat' or 'lat'. |
---|
51 | ; for Z axis, default name must be 'z', 'level', 'lev', 'depth...' |
---|
52 | ; |
---|
53 | ; [XYZ]MINMESH: to define the common variables i[xyz]minmesh |
---|
54 | ; used to define the grid only in a zoomed part of the original |
---|
55 | ; grid. Defaut values are 0L |
---|
56 | ; |
---|
57 | ; [XYZ]MAXMESH: to define the common variables i[xyz]maxmesh |
---|
58 | ; used to define the grid only in a zoomed part of the original |
---|
59 | ; grid. Defaut values are jp[ijk]glo-1 |
---|
60 | ; |
---|
61 | ; /xyindex: to define the x/y axis with index instead of using |
---|
62 | ; the values contained in X/YAXISNAME. |
---|
63 | ; x/yaxis = keyword_set(start1) + findgen(jpi/jpj) |
---|
64 | ; this forces key_onearth = 0 |
---|
65 | ; |
---|
66 | ; /zindex: to define the z axis with index instead of using |
---|
67 | ; the values contained in ZAXISNAME. |
---|
68 | ; zaxis = keyword_set(start1) + findgen(jpk) |
---|
69 | ; |
---|
70 | ; OUTPUTS:none, except the grid parameters of the common.pro |
---|
71 | ; |
---|
72 | ; COMMON BLOCKS:common.pro |
---|
73 | ; |
---|
74 | ; SIDE EFFECTS:change the grid parameters of the common.pro |
---|
75 | ; |
---|
76 | ; RESTRICTIONS: |
---|
77 | ; the file must contain an x and an y axis. (1 ou 2 dimentional array) |
---|
78 | ; |
---|
79 | ; EXAMPLE: IDL> initncdf,'toto.nc',glam=[-180,180] |
---|
80 | ; |
---|
81 | ; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr) |
---|
82 | ; 8 May 2002 |
---|
83 | ;- |
---|
84 | ;------------------------------------------------------------ |
---|
85 | ;------------------------------------------------------------ |
---|
86 | ;------------------------------------------------------------ |
---|
87 | PRO initncdf, ncfilein, XAXISNAME = xaxisname, YAXISNAME = yaxisname $ |
---|
88 | , ZAXISNAME = zaxisname, MASKNAME = maskname $ |
---|
89 | , INVMASK = invmask, USEASMASK = useasmask $ |
---|
90 | , MISSING_VALUE = missing_value, START1 = start1 $ |
---|
91 | , XYINDEX = xyindex, ZINDEX = zindex $ |
---|
92 | , _EXTRA = ex |
---|
93 | ; |
---|
94 | ; |
---|
95 | compile_opt idl2, strictarrsubs |
---|
96 | ; |
---|
97 | @common |
---|
98 | ;---------------------------------------------------------- |
---|
99 | ; check the name of the file |
---|
100 | ncfile = isafile(FILENAME = ncfilein, IODIRECTORY = iodir, _extra = ex) |
---|
101 | if size(ncfile, /type) NE 7 then BEGIN |
---|
102 | print, 'initncdf cancelled' |
---|
103 | return |
---|
104 | endif |
---|
105 | ; if the file is stored on tape |
---|
106 | if !version.os_family EQ 'unix' then spawn, 'file '+ncfile+' > /dev/null' |
---|
107 | ;---------------------------------------------------------- |
---|
108 | ; open the file |
---|
109 | cdfid = ncdf_open(ncfile) |
---|
110 | ; what is inside the file |
---|
111 | inside = ncdf_inquire(cdfid) |
---|
112 | ;---------------------------------------------------------- |
---|
113 | ; name of the variables |
---|
114 | namevar = strarr(inside.nvars) |
---|
115 | for varid = 0, inside.nvars-1 do begin |
---|
116 | invar = ncdf_varinq(cdfid, varid) |
---|
117 | namevar[varid] = strlowcase(invar.name) |
---|
118 | ENDFOR |
---|
119 | ;---------------------------------------------------------- |
---|
120 | ; find the xaxis |
---|
121 | if keyword_set(xaxisname) then xaxisname = strlowcase(xaxisname) ELSE xaxisname = 'x' |
---|
122 | xvarid = where(namevar EQ xaxisname OR namevar EQ 'longitude' $ |
---|
123 | OR namevar EQ 'nav_lon' OR namevar EQ 'lon') |
---|
124 | xvarid = xvarid[0] |
---|
125 | if xvarid EQ -1 then begin |
---|
126 | print, 'the xaxis was not found, check the use of XAXISNAME keyword' |
---|
127 | stop |
---|
128 | endif |
---|
129 | ; get the size of xaxis |
---|
130 | xinq = ncdf_varinq(cdfid, xvarid) |
---|
131 | ncdf_diminq, cdfid, xinq.dim[0], blabla, jpifromx |
---|
132 | ; should we read or compute the xaxis? |
---|
133 | IF NOT keyword_set(xyindex) THEN BEGIN |
---|
134 | ; read the xaxis |
---|
135 | ncdf_varget, cdfid, xvarid, xaxis |
---|
136 | ; make sure of the shape of xaxis |
---|
137 | IF xinq.ndims GE 2 THEN BEGIN |
---|
138 | ncdf_diminq, cdfid, xinq.dim[1], blabla, jpjfromx |
---|
139 | xaxis = reform(xaxis, jpifromx, jpjfromx, /over) |
---|
140 | ENDIF |
---|
141 | ENDIF ELSE xaxis = keyword_set(start1) + findgen(jpifromx) |
---|
142 | ;---------------------------------------------------------- |
---|
143 | ; find the yaxis |
---|
144 | if keyword_set(yaxisname) then yaxisname = strlowcase(yaxisname) ELSE yaxisname = 'y' |
---|
145 | yvarid = where(namevar EQ yaxisname OR namevar EQ 'latitude' OR namevar EQ 'nav_lat' OR namevar EQ 'lat') |
---|
146 | yvarid = yvarid[0] |
---|
147 | if yvarid EQ -1 then begin |
---|
148 | print, 'the yaxis was not found, check the use of YAXISNAME keyword' |
---|
149 | stop |
---|
150 | endif |
---|
151 | ; get the size of yaxis and check it is ok with the values found for x |
---|
152 | yinq = ncdf_varinq(cdfid, yvarid) |
---|
153 | IF xinq.ndims GE 2 THEN BEGIN |
---|
154 | ncdf_diminq, cdfid, yinq.dim[0], blabla, jpifromy |
---|
155 | ncdf_diminq, cdfid, yinq.dim[1], blabla, jpjfromy |
---|
156 | IF jpifromy NE jpifromx THEN BEGIN |
---|
157 | print, 'xaxis and y axis do not have the same x dimension...' |
---|
158 | ENDIF |
---|
159 | ENDIF ELSE ncdf_diminq, cdfid, yinq.dim[0], blabla, jpjfromy |
---|
160 | IF n_elements(jpjfromx) NE 0 THEN BEGIN |
---|
161 | IF jpjfromy NE jpjfromx THEN BEGIN |
---|
162 | print, 'xaxis and y axis do not have the same y dimension...' |
---|
163 | ENDIF |
---|
164 | ENDIF |
---|
165 | ; should we read or compute the xaxis? |
---|
166 | IF NOT keyword_set(xyindex) THEN BEGIN |
---|
167 | ; read the yaxis |
---|
168 | ncdf_varget, cdfid, yvarid, yaxis |
---|
169 | ; make sure of the shape of xaxis |
---|
170 | IF xinq.ndims GE 2 THEN yaxis = reform(yaxis, jpifromy, jpjfromy, /over) |
---|
171 | ENDIF ELSE yaxis = keyword_set(start1) + findgen(jpjfromy) |
---|
172 | ;---------------------------------------------------------- |
---|
173 | ; find the zaxis |
---|
174 | if keyword_set(zaxisname) then zaxisname = strlowcase(zaxisname) ELSE zaxisname = 'z' |
---|
175 | zvarid = where(namevar EQ 'nav_lev' or namevar EQ zaxisname OR namevar EQ 'level' OR namevar EQ 'lev' OR strmid(namevar, 0, 5) EQ 'depth') |
---|
176 | zvarid = zvarid[0] |
---|
177 | if zvarid EQ -1 AND inside.ndims GT 3 then begin |
---|
178 | print, 'initncdf: the zaxis was not found..., check the the use of ZAXISNAME keyword if you whant to find one...' |
---|
179 | ; stop |
---|
180 | endif |
---|
181 | ; read the zaxis |
---|
182 | if zvarid NE -1 THEN ncdf_varget, cdfid, zvarid, zaxis |
---|
183 | IF keyword_set(zindex) THEN $ |
---|
184 | zaxis = keyword_set(start1) + findgen(n_elements(zaxis)) |
---|
185 | ;---------------------------------------------------------- |
---|
186 | ; mask |
---|
187 | CASE 1 OF |
---|
188 | keyword_set(maskname):BEGIN |
---|
189 | mskid = (where(namevar EQ strlowcase(maskname)))[0] |
---|
190 | if mskid NE -1 THEN BEGIN |
---|
191 | mskinq = ncdf_varinq(cdfid, mskid) |
---|
192 | ; is the mask variable containing the record dimension? |
---|
193 | withrcd = (where(mskinq.dim EQ inside.recdim))[0] |
---|
194 | IF withrcd NE -1 THEN BEGIN |
---|
195 | ; in order to read only the first record |
---|
196 | ; we need to get the size of each dimension |
---|
197 | count = replicate(1L, mskinq.ndims) |
---|
198 | FOR d = 0, mskinq.ndims -1 DO BEGIN |
---|
199 | IF d NE withrcd THEN BEGIN |
---|
200 | ncdf_diminq, cdfid, mskinq.dim[d], name, size |
---|
201 | count[d] = size |
---|
202 | ENDIF |
---|
203 | ENDFOR |
---|
204 | ; read the variable for the first record |
---|
205 | ncdf_varget, cdfid, mskid, tmask, count = count |
---|
206 | ENDIF ELSE ncdf_varget, cdfid, mskid, tmask |
---|
207 | ; check if we need to applay add_offset and scale factor |
---|
208 | FOR a = 0, mskinq.natts-1 DO BEGIN |
---|
209 | attname = ncdf_attname(cdfid, mskid, a) |
---|
210 | CASE strlowcase(attname) OF |
---|
211 | 'add_offset':ncdf_attget, cdfid, mskid, attname, add_offset |
---|
212 | 'scale_factor':ncdf_attget, cdfid, mskid, attname, scale_factor |
---|
213 | ELSE: |
---|
214 | ENDCASE |
---|
215 | ENDFOR |
---|
216 | IF n_elements(scale_factor) NE 0 THEN tmask = tmask*scale_factor |
---|
217 | IF n_elements(add_offset) NE 0 THEN tmask = tmask+add_offset |
---|
218 | if keyword_set(invmask) then tmask = 1-tmask |
---|
219 | tmask = byte(round(tmask)) |
---|
220 | ENDIF ELSE tmask = -1 |
---|
221 | END |
---|
222 | ;.................. |
---|
223 | keyword_set(useasmask):BEGIN |
---|
224 | mskid = (where(namevar EQ strlowcase(useasmask)))[0] |
---|
225 | if mskid NE -1 THEN BEGIN |
---|
226 | mskinq = ncdf_varinq(cdfid, mskid) |
---|
227 | ; is the mask variable containing the record dimension? |
---|
228 | withrcd = (where(mskinq.dim EQ inside.recdim))[0] |
---|
229 | IF withrcd NE -1 THEN BEGIN |
---|
230 | ; in order to read only the first record |
---|
231 | ; we need to get the size of each dimension |
---|
232 | count = replicate(1L, mskinq.ndims) |
---|
233 | FOR d = 0, mskinq.ndims -1 DO BEGIN |
---|
234 | IF d NE withrcd THEN BEGIN |
---|
235 | ncdf_diminq, cdfid, mskinq.dim[d], name, size |
---|
236 | count[d] = size |
---|
237 | ENDIF |
---|
238 | ENDFOR |
---|
239 | ; read the variable for the first record |
---|
240 | ncdf_varget, cdfid, mskid, tmask, count = count |
---|
241 | ENDIF ELSE ncdf_varget, cdfid, mskid, tmask |
---|
242 | ; check if we need to applay add_offset and scale factor |
---|
243 | FOR a = 0, mskinq.natts-1 DO BEGIN |
---|
244 | attname = ncdf_attname(cdfid, mskid, a) |
---|
245 | CASE strlowcase(attname) OF |
---|
246 | 'add_offset':ncdf_attget, cdfid, mskid, attname, add_offset |
---|
247 | 'scale_factor':ncdf_attget, cdfid, mskid, attname, scale_factor |
---|
248 | 'missing_value':IF n_elements(missing_value) EQ 0 THEN $ |
---|
249 | ncdf_attget, cdfid, mskid, attname, missing_value |
---|
250 | ELSE: |
---|
251 | ENDCASE |
---|
252 | ENDFOR |
---|
253 | IF n_elements(scale_factor) NE 0 THEN tmask = tmask*scale_factor |
---|
254 | IF n_elements(add_offset) NE 0 THEN tmask = tmask+add_offset |
---|
255 | IF n_elements(missing_value) NE 0 THEN BEGIN |
---|
256 | ; we have to take care of the float accuracy |
---|
257 | CASE 1 OF |
---|
258 | missing_value GE 1.e6:tmask = tmask LT (missing_value-10) |
---|
259 | missing_value LE -1.e6:tmask = tmask GT (missing_value-10) |
---|
260 | abs(missing_value) LE 1.e-6:tmask = abs(tmask) GT 1.e-6 |
---|
261 | ELSE:tmask = tmask NE missing_value |
---|
262 | ENDCASE |
---|
263 | if keyword_set(invmask) then tmask = 1-tmask |
---|
264 | ENDIF ELSE BEGIN |
---|
265 | tmask = finite(tmask) |
---|
266 | IF min(tmask) EQ 1 THEN BEGIN |
---|
267 | print, 'missing or nan values not found...' |
---|
268 | tmask = -1 |
---|
269 | ENDIF |
---|
270 | ENDELSE |
---|
271 | ENDIF ELSE tmask = -1 |
---|
272 | END |
---|
273 | ;.................. |
---|
274 | ELSE:tmask = -1 |
---|
275 | ENDCASE |
---|
276 | ncdf_close, cdfid |
---|
277 | ; |
---|
278 | ; compute the grid |
---|
279 | if zvarid EQ -1 then BEGIN |
---|
280 | computegrid, xaxis = xaxis, yaxis = yaxis $ |
---|
281 | , mask = tmask, onearth = 1b - keyword_set(xyindex), _EXTRA = ex |
---|
282 | ENDIF ELSE BEGIN |
---|
283 | computegrid, xaxis = xaxis, yaxis = yaxis, zaxis = zaxis $ |
---|
284 | , mask = tmask, onearth = 1b - keyword_set(xyindex), _EXTRA = ex |
---|
285 | ENDELSE |
---|
286 | IF n_elements(time) EQ 0 THEN time = 0 |
---|
287 | jpt = n_elements(time) |
---|
288 | ;---------------------------------------------------------- |
---|
289 | |
---|
290 | return |
---|
291 | end |
---|