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

;+
;
; @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 _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 $
                  , ROMSGRID=romsgrid, _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
  @cm_general                      ;   needed for iodir
;
; 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]
    IF xinq.ndims GE 2 THEN ncdf_diminq, cdfid, xinq.dim[1], blabla, jpjfromx
  ENDELSE
  IF arg_present(xdimname) THEN ncdf_diminq, cdfid, dimidx,  xdimname, jpifromx
;
  IF arg_present(xaxis) THEN BEGIN
    ncdf_diminq, cdfid, dimidx, blabla, jpifromx
; 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
      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 BEGIN
      ncdf_diminq, cdfid, yinq.dim[0], blabla, jpifromy
      dimidy = yinq.dim[1]
    ENDIF ELSE dimidy = yinq.dim[0]
  ENDELSE
  IF arg_present(ydimname) THEN ncdf_diminq, cdfid, dimidy,  ydimname, jpjfromy
;
  IF arg_present(yaxis) THEN BEGIN
    IF n_elements(jpifromy) NE 0 THEN BEGIN
      IF jpifromy NE jpifromx THEN BEGIN
        dummy = report('x/y axes do not have the same x dimension...')
        stop
      ENDIF
    ENDIF
    ncdf_diminq, cdfid, dimidy, blabla, jpjfromy
    IF n_elements(jpjfromx) NE 0 THEN BEGIN
      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
      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