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

;+
;
; @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 ad 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
;
; @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 :
;    1 the keyword missing_value if existing
;    2 the attribute 'missing_value' if existing
;    3 NaN values if existing
;
; @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')
;
; @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 $
                     , _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
  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(fileid) 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]
    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, mask, count = count
    ENDIF ELSE ncdf_varget, cdfid, mskid, mask
; check if we need to applay add_offset and scale factor
    ncdf_getatt, cdfid, mskid, add_offset = add, scale_factor = scl, missing_value = miss
    IF n_elements(missing_value) NE 0 THEN miss = missing_value

    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(miss) NE 0 THEN BEGIN
; we have to take care of the float accuracy
        CASE 1 OF
          miss GE 1.e6:mask = mask LT (miss-10)
          miss LE -1.e6:mask = mask GT (miss+10)
          abs(miss) LE 1.e-6:mask = abs(mask) GT 1.e-6
          ELSE:mask = mask NE miss
        ENDCASE
      ENDIF ELSE BEGIN
        mask = finite(mask)
        IF min(mask) EQ 1 THEN BEGIN
          ras = report( 'missing or nan values not found...')
          mask = -1
        ENDIF
      ENDELSE
    ENDIF

    mask = byte(round(mask))
    if keyword_set(invmask) then mask = 1b-mask

  ENDIF ELSE mask = -1

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

  return, mask
END