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