source: trunk/SRC/ReadWrite/ncdf_getmask.pro @ 470

;+
;
; @file_comments
; get the land/sea mask array from a NetCDF file
;
; @categories
; Read NetCDF file
;
; @param fileid {in}{required}{type=salar string or long}
; if fileid is a scalar string then it is the name of the file (with
; the full path) to be opened (in that case, the file will be opened
; and closed within ncdf_getmask).
;
; if fileid is a scalar then it is the id of the file return by a call
; to ncdf_open outside of ncdf_getmask (in that case, the file will
; NOT be opened and closed within ncdf_getmask)
;
; @keyword ADDSCL_BEFORE {default=0}{type=scalar: 0 or 1}
; put 1 to apply add_offset and scale factor on data before looking for
; missing values when using USEASMASK keyword
;
; @keyword INVMASK {default=0}{type=scalar: 0 or 1}
; Inverse the land/sea mask (that should have 0/1 values for land/sea): mask = 1-mask
;
; @keyword MASKNAME {type=string}
; A string giving the name of the variable in the file
; that contains the land/sea mask
;
; @keyword MISSING_VALUE {type=scalar}
; To define (or redefine if the attribute is already existing) the
; missing values used with USEASMASK keyword. Note that this value is
; not used if TESTOP keyword is given and contains 2 words.  
;
; @keyword USEASMASK {type=scalar string}
; 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 operator defined by TESTOP
; keyword (default NE) and the testing values defined as 
;   1) the second word of TESTOP if existing
;   2) MISSING_VALUE keyword
;   3) attribute missing_value or _fillvalue of the variable USEASMASK
;   4) !Values.f_nan (can be used only with NE and EQ operators)
;
; @keyword TESTOP {default='NE'} {type=scalar string, for example 'GT 0.5'}
; a string describing the type of test that will be done to define the
; mask. The test is performed on the variable specified by USEASMASK
; keyword.
;
; TESTOP can contain 1 or 2 words. The first word is the operator
; definition: "EQ" "NE" "GE" "GT" "LE" "LT" (default is NE). The
; second word define the testing value. If TESTOP contains only 1
; word, then the test value is denifed by
;   1) MISSING_VALUE keyword
;   2) attribute missing_value or _fillvalue of the variable USEASMASK
;   3) !Values.f_nan (can be used only with NE and EQ operators)
;
; @keyword XMINMESH {default=0L}{type=scalar}
;       Define common (cm_4mesh) variables ixminmesh used to define the localization
;       of the first point of the grid along the x direction in a zoom of the original grid
;
; @keyword YMINMESH {default=0L}{type=scalar}
;       Define common (cm_4mesh) variables iyminmesh used to define the localization
;       of the first point of the grid along the y direction in a zoom of the original grid
;
; @keyword ZMINMESH {default=0L}{type=scalar}
;       Define common (cm_4mesh) variables izminmesh used to define the localization
;       of the first point of the grid along the z direction in a zoom of the original grid
;
; @keyword XMAXMESH {default=jpiglo-1}{type=scalar}
;       Define common (cm_4mesh) variables ixmaxmesh used to define the localization
;       of the last point of the grid along the x direction in a zoom of the original grid
;       Note that if XMAXMESH < 0 then ixmaxmesh is defined as ixmaxmesh = jpiglo -1 + xmaxmesh
;
; @keyword YMAXMESH {default=jpjglo-1}{type=scalar}
;       Define common (cm_4mesh) variables iymaxmesh used to define the localization
;       of the last point of the grid along the y direction in a zoom of the original grid
;       Note that if YMAXMESH < 0 then iymaxmesh is defined as iymaxmesh = jpjglo -1 + ymaxmesh
;
; @keyword ZMAXMESH {default=jpkglo-1}{type=scalar}
;       Define common (cm_4mesh) variables izmaxmesh used to define the localization
;       of the last point of the grid along the z direction in a zoom of the original grid
;       Note that if ZMAXMESH < 0 then izmaxmesh is defined as izmaxmesh = jpkglo -1 + maxmesh
;
; @keyword 
; _EXTRA to be able to call ncdf_getmask with _extra keyword
;
; @returns
; the land/sea mask 2D or 3D array or -1 in case of error or mask absence
;
; @examples
;
;   IDL> mask = ncdf_getmask('HadISST1_1m_187001_200702_sst_reg1m.nc',useasmask = 'sst', missing_value = -1.00000e+30)
;
;   IDL> mask = ncdf_getmask('meshmaskORCA2.nc', maskname = 'tmask')
;
;   IDL> mask = ncdf_getmask('t106.nc', useasmask = 'SLM', testop = 'le 0.5')
;
; @history
; August 2007: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
FUNCTION ncdf_getmask, fileid, ADDSCL_BEFORE = addscl_before $
                       , MASKNAME = maskname, USEASMASK = useasmask $
                       , MISSING_VALUE = missing_value, INVMASK = invmask $
                       , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
                       , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
                       , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $
                       , TESTOP = testop, _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
  @cm_general                   ;   needed for iodir
;
  IF n_elements(xminmesh) NE 0 THEN ixminmesh = 0L > long(xminmesh[0]) ELSE ixminmesh  = 0l
  IF n_elements(yminmesh) NE 0 THEN iyminmesh = 0L > long(yminmesh[0]) ELSE iyminmesh  = 0l
  IF n_elements(zminmesh) NE 0 THEN izminmesh = 0L > long(zminmesh[0]) ELSE izminmesh  = 0l
;
  IF NOT (keyword_set(maskname) OR keyword_set(useasmask)) AND keyword_set(romsgrid) THEN maskname = 'mask_rho'
  IF NOT (keyword_set(maskname) OR keyword_set(useasmask)) THEN return, -1
;----------------------------------------------------------
; should we open the file?
  IF size(fileid, /type) EQ 7 THEN $
     cdfid = ncdf_open(isafile(fileid, title = 'which file must be open by ncdf_getmask?', IODIRECTORY = iodir, _extra = ex)) $
  ELSE cdfid = fileid
; what is inside the file
  inq = ncdf_inquire(cdfid)
; name of the variables
  namevar = strarr(inq.nvars)
  for varid = 0, inq.nvars-1 do begin
    invar = ncdf_varinq(cdfid, varid)
    namevar[varid] = strlowcase(invar.name)
  ENDFOR
;----------------------------------------------------------
  CASE 1 OF
    keyword_set(maskname):mskid = (where(namevar EQ strlowcase(maskname)))[0]
    keyword_set(useasmask):mskid = (where(namevar EQ strlowcase(useasmask)))[0]
  ENDCASE
;
  if mskid NE -1 THEN BEGIN
    mskinq = ncdf_varinq(cdfid, mskid)
; is the mask variable containing the record dimension?
    withrcd = (where(mskinq.dim EQ inq.recdim))[0]
; in order to read only the first record
; we need to get the size of each dimension
    offset = lonarr(mskinq.ndims)
    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
    
    IF n_elements(xmaxmesh) NE 0 THEN BEGIN 
      IF xmaxmesh GE 0 THEN ixmaxmesh = long(xmaxmesh[0]) ELSE ixmaxmesh = count[0] - 1L + long(xmaxmesh[0])
    ENDIF ELSE ixmaxmesh = count[0] - 1L
    ixmaxmesh  = 0 > ixmaxmesh < (count[0] - 1L) ; make sure ixmaxmesh is not too big
    jpiglo = count[0]
    offset[0] = ixminmesh
    count[0] = ixmaxmesh - ixminmesh + 1L
    jpi = count[0]

    IF n_elements(ymaxmesh) NE 0 THEN  BEGIN 
      IF ymaxmesh GE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh = count[1] - 1L + long(ymaxmesh[0])
    ENDIF ELSE iymaxmesh = count[1] - 1L 
    iymaxmesh  = 0 > iymaxmesh < (count[1] - 1L) ; make sure ixmaxmesh is not too big
    jpjglo = count[1]
    offset[1] = iyminmesh
    count[1] = iymaxmesh - iyminmesh + 1L
    jpj = count[1]

    IF (mskinq.ndims EQ 3 AND withrcd EQ -1) OR (mskinq.ndims EQ 4) THEN BEGIN
      IF n_elements(zmaxmesh) NE 0 THEN  BEGIN 
        IF zmaxmesh GE 0 THEN izmaxmesh = long(zmaxmesh[0]) ELSE izmaxmesh = count[2] - 1L + long(zmaxmesh[0])
      ENDIF ELSE izmaxmesh = count[2] - 1L
      izmaxmesh  = 0 > izmaxmesh < (count[2] - 1L) ; make sure ixmaxmesh is not too big
      jpkglo = count[2]
      offset[2] = izminmesh
      count[2] = izmaxmesh - izminmesh + 1L
      jpk = count[2]
    ENDIF

; read the variable for the first record
    ncdf_varget, cdfid, mskid, mask, count = count,  offset = offset
    
; get add_offset, scale factor and missing value attributes
    ncdf_getatt, cdfid, mskid, add_offset = add, scale_factor = scl, missing_value = miss
; do we apply add_offset and scale factor ?
    IF keyword_set(addscl_before) THEN BEGIN
      IF scl NE 1 THEN mask = mask * scl
      IF add NE 0 THEN mask = mask + add
    ENDIF

    IF keyword_set(useasmask)  THEN BEGIN

      IF n_elements(missing_value) NE 0 THEN miss = missing_value
      IF size(miss, /type) EQ 7 THEN miss = !values.f_nan

      IF NOT keyword_set(testop) THEN testop = 'NE'
      tmp = strsplit(testop, ' ', /extract)
      op = strupcase(tmp[0])
      IF op EQ 'EQ' THEN BEGIN
        op = 'NE'
        invmask = 1b - keyword_set(invmask) 
      ENDIF
      IF n_elements(tmp) EQ 1 THEN testval = miss ELSE testval = float(tmp[1]) 
      IF finite(testval) EQ 0 THEN BEGIN
        IF op NE 'NE' THEN mask = report('NaN test value can be used only with EQ or NE operator') ELSE mask = finite(mask)
      ENDIF ELSE BEGIN
        CASE op OF
          'GE':mask = mask GE testval
          'GT':mask = mask GT testval
          'LE':mask = mask LE testval
          'LT':mask = mask LT testval
          'NE':BEGIN 
; we have to take care of the float accuracy
            CASE 1 OF
              testval GE  1.e6:mask = mask LT (testval - 10)
              testval LE -1.e6:mask = mask GT (testval + 10)
              abs(testval) LE 1.e-6:mask = abs(mask) GT 1.e-6
              ELSE:mask = mask NE testval
            ENDCASE
          END
        ENDCASE
      ENDELSE

    ENDIF

    IF mask[0] NE -1 THEN BEGIN 
      mask = byte(round(mask))
      if keyword_set(invmask) then mask = 1b - mask
    ENDIF

  ENDIF ELSE mask = -1

  IF size(fileid, /type) EQ 7 THEN ncdf_close, cdfid

  return, mask
END