Index: trunk/SRC/ReadWrite/ncdf_getaxis.pro
===================================================================
--- trunk/SRC/ReadWrite/ncdf_getaxis.pro	(revision 221)
+++ trunk/SRC/ReadWrite/ncdf_getaxis.pro	(revision 221)
@@ -0,0 +1,223 @@
+;+
+; @file_comments
+; get the x/y dimension Id and x/y axes from a netcdf file
+;
+; @categories
+; Reading
+;
+; @param cdfid {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 XDIRNAME {default='longitude', 'lon', 'x', 'longitude*', 'lon*', 'x*', '*longitude*', '*lon*' or '*x*'}{type=scalar string}
+; A string giving the name of the x dimension
+;      
+; @keyword YDIRNAME {default='latitude', 'lat', 'y', 'latitude*', 'lat*', 'y*', 'eta_*', '*latitude*', '*lat*', '*y*'}{type=scalar string}
+; A string giving the name of the y dimension
+;      
+; @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 [xyz]axis. 
+;      
+; @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 [xyz]axis. 
+;
+; @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 XYINDEX {default=0}{type=scalar: 0 or 1}
+;
+; @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, cdfid, dimidx, dimidy, xaxis, yaxis $
+                  , XAXISNAME = xaxisname, YAXISNAME = yaxisname $
+                  , XDIMNAME = xdimname, YDIMNAME = ydimname $
+                  , XYINDEX = xyindex, START1 = start1, ROMSGRID = romsgrid, _extra = ex
+;
+  compile_opt idl2, strictarrsubs
+;
+
+; 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)
+; 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
+    ENDIF
+    romsgrid = 0b   
+  ENDIF ELSE BEGIN 
+    romsgrid = strmid(namevar[xvarid], 0, 4) EQ 'lon_'
+    xinq = ncdf_varinq(cdfid, xvarid)
+    dimidx = xinq.dim[0]
+    IF xinq.ndims GE 2 THEN ncdf_diminq, cdfid, xinq.dim[1], blabla, jpjfromx
+  ENDELSE
+;
+  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
+; 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)
+; 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)
+    IF yinq.ndims GE 2 THEN BEGIN
+      ncdf_diminq, cdfid, yinq.dim[0], blabla, jpifromy
+      IF jpifromy NE jpifromx THEN BEGIN 
+        dummy = report('x/y axes do not have the same x dimension...')
+        stop
+      ENDIF
+      dimidy = yinq.dim[1]
+    ENDIF ELSE dimidy = yinq.dim[0]
+  ENDELSE
+;
+  IF arg_present(yaxis) THEN BEGIN
+    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
+; make sure of the shape of xaxis
+      IF n_elements(jpifromy) NE 0 THEN yaxis = reform(yaxis, jpifromy, jpjfromy, /over)
+    ENDELSE
+  ENDIF
+
+  return
+END
