source: trunk/ToBeReviewed/INIT/initncdf.pro @ 69

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:initncdf
;
; PURPOSE:initfile for Netcdf file. define all the grid parameters
;
; CATEGORY:
;
; CALLING SEQUENCE:initncdf, ncfilename
; 
; INPUTS:ncfilename: a string giving the name of the NetCdf file
;
; KEYWORD PARAMETERS:
;
;       GLAMBOUNDARY:a 2 elements vector, {lon1,lon2], the longitute
;       boundaries that should be used to visualize the data.
;         lon2 > lon1
;         lon2 - lon1 le 360
;       key_shift will be defined according to GLAMBOUNDARY.
;
;       /INVMASK: to inverse the mask: mask = 1-mask
;
;       IODIRECTORY;a string giving the name of iodirectory (see
;       isafile.pro for all possibilities). default value is common 
;       variable iodir
;
;       MASKNAME: a string giving the name of the variable in the file 
;       that contains the land/sea mask
;
;       MISSING_VALUE: to define (or redifine if the attribute is
;       already existing) the missing values used with USEASMASK
;       keyword
;
;       /start1: index the axis from 1 instead of 0 when using
;       /xyindex and/or zindex
;
;       USEASMASK: a string giving the name of the variable in the file 
;       that will be used to build the land/sea mask. In this case the
;       mask is based on the first record (if record dimension
;       exists). The mask is build according to :
;         1 the keyword missing_value if existing
;         2 the attribute 'missing_value' if existing
;         3 NaN values if existing
;
;       [XYZ]AXISNAME= a string giving the name of the variable in the file 
;       that contains the [xyz]axis. 
;       for X axis, default name must be 'x', 'longitude', 'nav_lon' or 'lon'.
;       for Y axis, default name must be 'y', 'latitude', 'nav_lat' or 'lat'.
;       for Z axis, default name must be 'z', 'level', 'lev', 'depth...'
;
;       [XYZ]MINMESH: to define the common variables i[xyz]minmesh
;       used to define the grid only in a zoomed part of the original
;       grid. Defaut values are 0L
;
;       [XYZ]MAXMESH: to define the common variables i[xyz]maxmesh
;       used to define the grid only in a zoomed part of the original
;       grid. Defaut values are jp[ijk]glo-1
;
;       /xyindex: to define the x/y axis with index instead of using
;       the values contained in X/YAXISNAME. 
;       x/yaxis = keyword_set(start1) + findgen(jpi/jpj) 
;       this forces key_onearth = 0
;
;       /zindex: to define the z axis with index instead of using
;       the values contained in ZAXISNAME. 
;       zaxis = keyword_set(start1) + findgen(jpk) 
;
; OUTPUTS:none, except the grid parameters of the common.pro
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:change the grid parameters of the common.pro
;
; RESTRICTIONS:
; the file must contain an x and an y axis. (1 ou 2 dimentional array)
;
; EXAMPLE: IDL> initncdf,'toto.nc',glam=[-180,180]
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      8 May 2002
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
PRO initncdf, ncfilein, XAXISNAME = xaxisname, YAXISNAME = yaxisname $
              , ZAXISNAME = zaxisname, MASKNAME = maskname $
              , INVMASK = invmask, USEASMASK = useasmask $
              , MISSING_VALUE = missing_value, START1 = start1 $
              , XYINDEX = xyindex, ZINDEX = zindex $
              , _EXTRA = ex
;
@common
;----------------------------------------------------------
; check the name of the file
  ncfile = isafile(FILENAME = ncfilein, IODIRECTORY = iodir, _extra = ex)
  if size(ncfile, /type) NE 7 then BEGIN
    print, 'initncdf cancelled'
    return
  endif
; if the file is stored on tape
  if !version.os_family EQ 'unix' then spawn, 'file '+ncfile+' > /dev/null' 
;----------------------------------------------------------
; open the file
  cdfid = ncdf_open(ncfile)
; what is inside the file
  inside = ncdf_inquire(cdfid)
;----------------------------------------------------------
; name of the variables
  namevar = strarr(inside.nvars)
  for varid = 0, inside.nvars-1 do begin
    invar = ncdf_varinq(cdfid, varid)
    namevar[varid] = strlowcase(invar.name)
  ENDFOR
;----------------------------------------------------------
; find the xaxis
  if keyword_set(xaxisname) then xaxisname = strlowcase(xaxisname) ELSE xaxisname = 'x'
  xvarid = where(namevar EQ xaxisname OR namevar EQ 'longitude' $
                 OR namevar EQ 'nav_lon' OR namevar EQ 'lon')
  xvarid = xvarid[0]
  if xvarid EQ -1 then begin
    print, 'the xaxis was not found, check the use of XAXISNAME keyword'
    stop
  endif
; get the size of xaxis
  xinq = ncdf_varinq(cdfid, xvarid)
  ncdf_diminq, cdfid, xinq.dim[0], blabla, jpifromx
; should we read or compute the xaxis?
  IF NOT keyword_set(xyindex) THEN BEGIN
; read the xaxis
    ncdf_varget, cdfid, xvarid, xaxis
; make sure of the shape of xaxis
    IF xinq.ndims GE 2 THEN BEGIN
      ncdf_diminq, cdfid, xinq.dim[1], blabla, jpjfromx
      xaxis = reform(xaxis, jpifromx, jpjfromx, /over)
    ENDIF
  ENDIF ELSE xaxis = keyword_set(start1) + findgen(jpifromx)
;----------------------------------------------------------
; find the yaxis
  if keyword_set(yaxisname) then yaxisname = strlowcase(yaxisname) ELSE yaxisname = 'y'
  yvarid = where(namevar EQ yaxisname OR namevar EQ 'latitude' OR namevar EQ 'nav_lat' OR namevar EQ 'lat')
  yvarid = yvarid[0]
  if yvarid EQ -1 then begin
    print, 'the yaxis was not found, check the use of YAXISNAME keyword'
    stop
  endif
; get the size of yaxis and check it is ok with the values found for x
  yinq = ncdf_varinq(cdfid, yvarid)
  IF xinq.ndims GE 2 THEN BEGIN
    ncdf_diminq, cdfid, yinq.dim[0], blabla, jpifromy
    ncdf_diminq, cdfid, yinq.dim[1], blabla, jpjfromy
    IF jpifromy NE jpifromx THEN BEGIN 
      print, 'xaxis and y axis do not have the same x dimension...'
    ENDIF
  ENDIF ELSE ncdf_diminq, cdfid, yinq.dim[0], blabla, jpjfromy
  IF n_elements(jpjfromx) NE 0 THEN BEGIN
    IF jpjfromy NE jpjfromx THEN BEGIN 
      print, 'xaxis and y axis do not have the same y dimension...'
    ENDIF    
  ENDIF
; should we read or compute the xaxis?
  IF NOT keyword_set(xyindex) THEN BEGIN
; read the yaxis
    ncdf_varget, cdfid, yvarid, yaxis
; make sure of the shape of xaxis
    IF xinq.ndims GE 2 THEN yaxis = reform(yaxis, jpifromy, jpjfromy, /over)
  ENDIF ELSE yaxis = keyword_set(start1) + findgen(jpjfromy)
;----------------------------------------------------------
; find the zaxis
  if keyword_set(zaxisname) then zaxisname = strlowcase(zaxisname) ELSE zaxisname = 'z'
  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')
  zvarid = zvarid[0]
  if zvarid EQ -1 AND inside.ndims GT 3 then begin
    print, 'initncdf: the zaxis was not found..., check the the use of ZAXISNAME keyword if you whant to find one...'
;     stop
  endif
; read the zaxis
  if zvarid NE -1 THEN ncdf_varget, cdfid, zvarid, zaxis
  IF keyword_set(zindex) THEN $
    zaxis = keyword_set(start1) + findgen(n_elements(zaxis))
;----------------------------------------------------------
; mask
  CASE 1 OF
    keyword_set(maskname):BEGIN
      mskid = (where(namevar EQ strlowcase(maskname)))[0]
      if mskid NE -1 THEN BEGIN 
        mskinq = ncdf_varinq(cdfid, mskid)
; is the mask variable containing the record dimension?
        withrcd = (where(mskinq.dim EQ inside.recdim))[0]
        IF withrcd NE -1 THEN BEGIN
; in order to read only the first record          
; we need to get the size of each dimension
          count = replicate(1L, mskinq.ndims)
          FOR d = 0, mskinq.ndims -1 DO BEGIN
            IF d NE withrcd THEN BEGIN
              ncdf_diminq, cdfid, mskinq.dim[d], name, size
              count[d] = size
            ENDIF
          ENDFOR
; read the variable for the first record          
          ncdf_varget, cdfid, mskid, tmask, count = count
        ENDIF ELSE ncdf_varget, cdfid, mskid, tmask
; check if we need to applay add_offset and scale factor         
        FOR a = 0, mskinq.natts-1 DO BEGIN
          attname = ncdf_attname(cdfid, mskid, a)     
          CASE strlowcase(attname) OF
            'add_offset':ncdf_attget, cdfid, mskid, attname, add_offset
            'scale_factor':ncdf_attget, cdfid, mskid, attname, scale_factor
            ELSE:
          ENDCASE
        ENDFOR
        IF n_elements(scale_factor) NE 0 THEN tmask = tmask*scale_factor
        IF n_elements(add_offset) NE 0 THEN tmask = tmask+add_offset
        if keyword_set(invmask) then tmask = 1-tmask
        tmask = byte(round(tmask))
      ENDIF ELSE tmask = -1       
    END
;..................
    keyword_set(useasmask):BEGIN
      mskid = (where(namevar EQ strlowcase(useasmask)))[0]
      if mskid NE -1 THEN BEGIN 
        mskinq = ncdf_varinq(cdfid, mskid)
; is the mask variable containing the record dimension?
        withrcd = (where(mskinq.dim EQ inside.recdim))[0]
        IF withrcd NE -1 THEN BEGIN
; in order to read only the first record          
; we need to get the size of each dimension
          count = replicate(1L, mskinq.ndims)
          FOR d = 0, mskinq.ndims -1 DO BEGIN
            IF d NE withrcd THEN BEGIN
              ncdf_diminq, cdfid, mskinq.dim[d], name, size
              count[d] = size
            ENDIF
          ENDFOR
; read the variable for the first record       
          ncdf_varget, cdfid, mskid, tmask, count = count
        ENDIF ELSE ncdf_varget, cdfid, mskid, tmask
; check if we need to applay add_offset and scale factor       
        FOR a = 0, mskinq.natts-1 DO BEGIN
          attname = ncdf_attname(cdfid, mskid, a)     
          CASE strlowcase(attname) OF
            'add_offset':ncdf_attget, cdfid, mskid, attname, add_offset
            'scale_factor':ncdf_attget, cdfid, mskid, attname, scale_factor
            'missing_value':IF n_elements(missing_value) EQ 0 THEN $
              ncdf_attget, cdfid, mskid, attname, missing_value
            ELSE:
          ENDCASE
        ENDFOR
        IF n_elements(scale_factor) NE 0 THEN tmask = tmask*scale_factor
        IF n_elements(add_offset) NE 0 THEN tmask = tmask+add_offset
        IF n_elements(missing_value) NE 0 THEN BEGIN
; we have to take care of the float accuracy
          CASE 1 OF
            missing_value GE 1.e6:tmask = tmask LT (missing_value-10)
            missing_value LE -1.e6:tmask = tmask GT (missing_value-10)
            abs(missing_value) LE 1.e-6:tmask = abs(tmask) GT 1.e-6
            ELSE:tmask = tmask NE missing_value
          ENDCASE
          if keyword_set(invmask) then tmask = 1-tmask
        ENDIF ELSE BEGIN
          tmask = finite(tmask)
          IF min(tmask) EQ 1 THEN BEGIN
            print, 'missing or nan values not found...'
            tmask = -1          
          ENDIF
        ENDELSE
      ENDIF ELSE tmask = -1       
    END
;..................
    ELSE:tmask  = -1
  ENDCASE
  ncdf_close, cdfid
;
; compute the grid
  if zvarid EQ -1 then BEGIN 
    computegrid, xaxis = xaxis, yaxis = yaxis $
      , mask = tmask, onearth = 1b - keyword_set(xyindex), _EXTRA = ex 
  ENDIF ELSE BEGIN 
    computegrid, xaxis = xaxis, yaxis = yaxis, zaxis = zaxis $
      , mask = tmask, onearth = 1b - keyword_set(xyindex), _EXTRA = ex
  ENDELSE 
  IF n_elements(time) EQ 0 THEN time = 0
  jpt = n_elements(time) 
;----------------------------------------------------------

  return
end