source: trunk/SRC/ReadWrite/ncdf_getaxis.pro @ 399

;+
;
; @file_comments
; get the x/y dimension Id and x/y axes from a netcdf file
;
; @categories
; Read NetCDF file
;
; @param fileid {in}{required}{type=scalar}
; the id of the netcdf file
;
; @param dimidx {out}{type=scalar (long)}
; id of the x dimension
;
; @param dimidy {out}{type=scalar (long)}
; id of the y dimension
;
; @param xaxis {out}{type=1D or 2D array}
; the x axis
;
; @param yaxis {out}{type=1D or 2D array}
; the y axis
;
; @keyword ROMSGRID {out}{type=scalar: 0 or 1}
; gives back if we are using a ROMS grid (1) or not (0)
;
; @keyword START1 {default=0}{type=scalar: 0 or 1}
; Index the axis from 1 instead of 0 when using /xyindex
;
; @keyword XDIMNAME {default='longitude', 'lon', 'x', 'longitude*', 'lon*', 'x*', '*longitude*', '*lon*' or '*x*'}{type=scalar string}
; A string giving the name of the x dimension or/and a named variable
; in which x dimension name is returned.
;
; @keyword YDIMNAME {default='latitude', 'lat', 'y', 'latitude*', 'lat*', 'y*', 'eta_*', '*latitude*', '*lat*', '*y*'}{type=scalar string}
; A string giving the name of the y dimension or/and a named variable
; in which y dimension name is returned.
;
; @keyword XAXISNAME {default='x', 'longitude', 'nav_lon', 'lon', 'lon_rho' or 'NbLongitudes'}{type=scalar string}
; A string giving the name of the variable in the file
; that contains the x axis or/and a named variable
; in which this variable name is returned.
;
; @keyword YAXISNAME {default='y', 'latitude', 'nav_lat','lat', 'lat_rho' or 'NbLatitudes'}{type=scalar string}
; A string giving the name of the variable in the file
; that contains the y axis or/and a named variable
; in which this variable name is returned.
;
; @keyword XYINDEX {default=0}{type=scalar: 0 or 1}
; 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)
;
; @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 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 _EXTRA
; Used to be able to call ncdf_getaxis with _extra
;
; @history
; March 2007: Sebastien Masson (smasson\@locean-ipsl.upmc.fr)
;
;
; @version
; $Id$
;-
PRO ncdf_getaxis, fileid, dimidx, dimidy, xaxis, yaxis $
                  , XAXISNAME = xaxisname, YAXISNAME = yaxisname $
                  , XDIMNAME = xdimname, YDIMNAME = ydimname $
                  , XYINDEX = xyindex, START1 = start1 $
                  , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $
                  , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $
                  , ROMSGRID = romsgrid, _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
@cm_general                   ;   needed for iodir
@cm_4mesh
;
  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
; 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_getaxis?', IODIRECTORY = iodir, _extra = ex)) $
  ELSE cdfid = fileid
; what is inside the file
  inside = ncdf_inquire(cdfid)
;------------------------------------------------------------
; name of all dimensions
  namedim = strarr(inside.ndims)
  for dimiq = 0, inside.ndims-1 do begin
    ncdf_diminq, cdfid, dimiq, tmpname, value
    namedim[dimiq] = strlowcase(tmpname)
  ENDFOR
;----------------------------------------------------------
; 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
; try to get the variable that contains 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' $
                  OR namevar EQ 'lon_rho' OR namevar EQ 'nblongitudes'))[0]
; no xaxis variable found, we will build a fake xaxis based on the size of the x dimension
; -> we must find the x dimension
  IF xvarid EQ -1 THEN BEGIN
    dummy = report(['xaxis variable was not found within the default names:' $
                    , '''longitude'', ''nav_lon'', ''lon'', ''lon_rho'', ''nblongitudes''' $
                    , ' we use a fake xaxis based on x dimension size (or use XAXISNAME keyword)'], /simple)
    xaxisname = 'Not Found'
; try to get the dimension corresponding to x
; roms file?
    dimidx = where(namedim EQ 'xi_rho' OR namedim EQ 'xi_u' OR namedim EQ 'xi_v' OR namedim EQ 'xi_psi')
    IF dimidx[0] EQ -1 THEN BEGIN
; we are looking for a x dimension with a name matching one of the following regular expression:
      if keyword_set(xdimname) then testname = strlowcase(xdimname) $
      ELSE testname = ['longitude', 'lon', 'x', 'longitude*', 'lon*', 'x*', '*longitude*', '*lon*', '*x*']
      cnt = -1
      ii = 0
      WHILE cnt NE 1 AND ii LT n_elements(testname) DO BEGIN
        dimidx = where(strmatch(namedim, testname[ii]) EQ 1, cnt)
        ii = ii+1
      ENDWHILE
      CASE cnt OF
        0:begin
          dummy = report(['none of the dimensions name matches one of the following regular expression:' $
                          , '''longitude'', ''lon'', ''x'', ''longitude*'', ''lon*'', ''x*'', ''*longitude*'', ''*lon*'', ''*x*''' $
                          , ' => we cannot find the x dimension, use XDIMNAME keyword'], /simple)
          stop
        END
        1:dimidx = dimidx[0]
        ELSE:begin
          dummy = report(['several (and not one unique) dimensions name matches the following regular expression:' $
                          , '''longitude'', ''lon'', ''x'', ''longitude*'', ''lon*'', ''x*'', ''*longitude*'', ''*lon*'', ''*x*''' $
                          , ' => we cannot find the x dimension, use XDIMNAME keyword'], /simple)
          stop
        ENDELSE
      ENDCASE
      romsgrid = 0b
    ENDIF ELSE romsgrid = 1b
  ENDIF ELSE BEGIN
    romsgrid = strmid(namevar[xvarid], 0, 4) EQ 'lon_'
    xinq = ncdf_varinq(cdfid, xvarid)
    xaxisname = xinq.name
    dimidx = xinq.dim[0]
    xoffset = lonarr(xinq.ndims)
    xcount = lonarr(xinq.ndims)
    FOR i = 0, xinq.ndims -1 DO BEGIN
      ncdf_diminq, cdfid, xinq.dim[i], blabla, tmpsz
      xcount[i] = tmpsz
    ENDFOR
    jpiglo = xcount[0]
    IF n_elements(xmaxmesh) NE 0 THEN BEGIN 
      IF xmaxmesh GE 0 THEN ixmaxmesh = long(xmaxmesh[0]) ELSE ixmaxmesh = xcount[0] -1 + long(xmaxmesh[0])
    ENDIF ELSE ixmaxmesh = xcount[0] - 1L
    ixmaxmesh  = 0 > ixmaxmesh < (xcount[0] - 1L) ; make sure ixmaxmesh is not too big
    xoffset[0] = ixminmesh
    xcount[0] = ixmaxmesh - ixminmesh + 1L
    jpifromx = xcount[0]
    jpi = jpifromx
    IF xinq.ndims GE 2 THEN BEGIN 
      jpjglo = xcount[1]
      IF n_elements(ymaxmesh) NE 0 THEN BEGIN 
        IF ymaxmesh GE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh = xcount[1] - 1L + long(ymaxmesh[0])
      ENDIF ELSE iymaxmesh = xcount[1] - 1L
      iymaxmesh  = 0 > iymaxmesh < (xcount[1] - 1L) ; make sure ixmaxmesh is not too big
      xoffset[1] = iyminmesh
      xcount[1] = iymaxmesh - iyminmesh + 1L
      jpjfromx = xcount[1]
      jpj = jpjfromx
    ENDIF
  ENDELSE
  IF arg_present(xdimname) THEN ncdf_diminq, cdfid, dimidx,  xdimname, tmp
;
  IF arg_present(xaxis) THEN BEGIN
; should we read or compute the xaxis?
    IF keyword_set(xyindex) OR xvarid EQ -1 THEN BEGIN
      xaxis = keyword_set(start1) + findgen(jpifromx)
      xyindex = 1
    ENDIF ELSE BEGIN
; read the xaxis
      ncdf_varget, cdfid, xvarid, xaxis, offset = xoffset, count = xcount
      ncdf_getatt, cdfid, xvarid, ADD_OFFSET = add_offset, SCALE_FACTOR = scale_factor
      IF scale_factor NE 1. THEN xaxis = temporary(xaxis) * scale_factor
      IF add_offset NE 0. THEN xaxis = temporary(xaxis) + add_offset
; make sure of the shape of xaxis
      IF n_elements(jpjfromx) NE 0 THEN xaxis = reform(xaxis, jpifromx, jpjfromx, /over)
    ENDELSE
  ENDIF

;----------------------------------------------------------
; find the yaxis
; try to get the variable that contains 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' $
                  OR namevar EQ 'lat_rho' OR namevar EQ 'nblatitudes'))[0]
  yvarid = yvarid[0]
; no yaxis variable found, we will build a fake yaxis based on the size of the y dimension
; -> we must find the y dimension
  if yvarid EQ -1 then begin
    dummy = report(['yaxis variable was not found within the default names:' $
                    , '''latitude'', ''nav_lat'', ''lat'', ''lat_rho'', ''nblatitudes''' $
                    , ' we use a fake yaxis based on y dimension size (or use YAXISNAME keyword)'], /simple)
    yaxisname = 'Not Found'
; try to get the dimension corresponding to y
; roms file?
    dimidy = where(namedim EQ 'eta_rho' OR namedim EQ 'eta_u' OR namedim EQ 'eta_v' OR namedim EQ 'eta_psi')
    IF dimidy[0] EQ -1 THEN BEGIN
; we are looking for a y dimension with a name matching one of the following regular expression:
      if keyword_set(ydimname) then testname = strlowcase(ydimname) $
      ELSE testname = ['latitude', 'lat', 'y', 'latitude*', 'lat*', 'y*', 'eta_*', '*latitude*', '*lat*', '*y*']
      cnt = -1
      ii = 0
      WHILE cnt NE 1 AND ii LT n_elements(testname) DO BEGIN
        dimidy = where(strmatch(namedim, testname[ii]) EQ 1, cnt)
        ii = ii+1
      ENDWHILE
      CASE cnt OF
        0:begin
          dummy = report(['none of the dimensions name matches one of the following regular expression:' $
                          , '''latitude'', ''lat'', ''y'', ''latitude*'', ''lat*'', ''y*'', ''eta_*'', ''*latitude*'', ''*lat*'', ''*y*''' $
                          , ' => we cannot find the y dimension, use YDIMNAME keyword'], /simple)
          stop
        END
        1:dimidy = dimidy[0]
        ELSE:begin
          dummy = report(['several (and not one unique) dimensions name matches the following regular expression:' $
                          , '''latitude'', ''lat'', ''y'', ''latitude*'', ''lat*'', ''y*'', ''eta_*'', ''*latitude*'', ''*lat*'', ''*y*''' $
                          , ' => we cannot find the y dimension, use YDIMNAME keyword'], /simple)
          stop
        ENDELSE
      ENDCASE
    ENDIF
  ENDIF ELSE BEGIN
    yinq = ncdf_varinq(cdfid, yvarid)
    yaxisname = yinq.name
    IF yinq.ndims GE 2 THEN dimidy = yinq.dim[1] ELSE dimidy = yinq.dim[0]
    yoffset = lonarr(yinq.ndims)
    ycount = lonarr(yinq.ndims)
    FOR i = 0, yinq.ndims -1 DO BEGIN
      ncdf_diminq, cdfid, yinq.dim[i], blabla, tmpsz
      ycount[i] = tmpsz
    ENDFOR
    idy = yinq.ndims GE 2
    jpjglo = ycount[idy]
    IF n_elements(ymaxmesh) NE 0 THEN BEGIN 
      IF ymaxmesh GE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh = ycount[idy] - 1L + long(ymaxmesh[0])
    ENDIF ELSE iymaxmesh = ycount[idy] - 1L
    iymaxmesh  = iymaxmesh < (ycount[idy] - 1L) ; make sure ixmaxmesh is not too big
    yoffset[idy] = iyminmesh
    ycount[idy] = iymaxmesh - iyminmesh + 1L
    jpjfromy = ycount[idy]
    jpj = jpjfromy
    IF yinq.ndims GE 2 THEN BEGIN 
      jpiglo = ycount[0]
      IF n_elements(xmaxmesh) NE 0 THEN BEGIN 
        IF xmaxmesh GE 0 THEN ixmaxmesh = long(xmaxmesh[0]) ELSE ixmaxmesh = ycount[0] - 1L + long(xmaxmesh[0])
      ENDIF ELSE ixmaxmesh = ycount[0] - 1L
      ixmaxmesh  = ixmaxmesh < (ycount[0] - 1L) ; make sure ixmaxmesh is not too big
      yoffset[0] = ixminmesh
      ycount[0] = ixmaxmesh - ixminmesh + 1L
      jpifromy = xcount[0]
      jpi = jpifromy
    ENDIF
  ENDELSE
  IF arg_present(ydimname) THEN ncdf_diminq, cdfid, dimidy,  ydimname, tmp
;
  IF arg_present(yaxis) THEN BEGIN
    IF n_elements(jpifromy) NE 0 THEN BEGIN
      IF n_elements(jpifromx) EQ 0 THEN ncdf_diminq, cdfid, dimidx, blabla, jpifromx
      IF jpifromy NE jpifromx THEN BEGIN
        dummy = report('x/y axes do not have the same x dimension...')
        stop
      ENDIF
    ENDIF
    IF n_elements(jpjfromx) NE 0 THEN BEGIN
      IF n_elements(jpjfromy) EQ 0 THEN ncdf_diminq, cdfid, dimidy, blabla, jpjfromy
      IF jpjfromy NE jpjfromx THEN BEGIN
        dummy = report(' x/y axes do not have the same y dimension...')
        stop
      ENDIF
    ENDIF
; should we read or compute the xaxis?
    IF keyword_set(xyindex) OR yvarid EQ -1 THEN BEGIN
      yaxis = keyword_set(start1) + findgen(jpjfromy)
    ENDIF ELSE BEGIN
; read the yaxis
      ncdf_varget, cdfid, yvarid, yaxis, offset = yoffset, count = ycount
      ncdf_getatt, cdfid, yvarid, ADD_OFFSET = add_offset, SCALE_FACTOR = scale_factor
      IF scale_factor NE 1. THEN yaxis = temporary(yaxis) * scale_factor
      IF add_offset NE 0. THEN yaxis = temporary(yaxis) + add_offset
; make sure of the shape of xaxis
      IF n_elements(jpifromy) NE 0 THEN yaxis = reform(yaxis, jpifromy, jpjfromy, /over)
    ENDELSE
  ENDIF

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

  return
END