Changeset 172 for trunk/SRC/ToBeReviewed/WIDGET/AUTOUR_de_XXX/scanfile.pro

Timestamp:

09/11/06 09:11:26 (18 years ago)

Author:

smasson

Message:

bugfix + manage roms outputs

File:

• trunk/SRC/ToBeReviewed/WIDGET/AUTOUR_de_XXX/scanfile.pro (modified) (4 diffs)

trunk/SRC/ToBeReviewed/WIDGET/AUTOUR_de_XXX/scanfile.pro

@@ -89 +89 @@
 ;------------------------------------------------------------
   infile = ncdf_inquire(cdfid)  ;
-; find vargrid ...
-  IF keyword_set(grid) THEN vargrid = strupcase(grid) ELSE BEGIN
-    vargrid = 'T'               ; default definition
-    IF finite(glamu[0]) EQ 1 THEN BEGIN
+;------------------------------------------------------------
+; name of all dimensions
+;------------------------------------------------------------
+  namedim = strarr(infile.ndims)
+  for dimiq = 0, infile.ndims-1 do begin
+    ncdf_diminq, cdfid, dimiq, tmpname, value 
+    namedim[dimiq] = strlowcase(tmpname)
+  ENDFOR
+; 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:
+    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'])
+        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'])
+        stop
+      ENDELSE
+    ENDCASE
+  ENDIF
+; 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:
+    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'])
+        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 x dimension'])
+        stop
+      ENDELSE
+    ENDCASE
+  ENDIF
+;------------------------------------------------------------
+; name of all variables
+;------------------------------------------------------------
+; we keep only the variables containing at least x, y and time dimension (if existing...)
+  namevar = strarr(infile.nvars)
+  for varid = 0, infile.nvars-1 do begin
+    invar = ncdf_varinq(cdfid, varid) ; what contains the variable?
+    if (inter(invar.dim, dimidx))[0] NE -1 AND $
+       (inter(invar.dim, dimidy))[0] NE -1 AND $
+       ((where(invar.dim EQ infile.recdim))[0] NE -1 OR infile.recdim EQ -1) $
+    THEN namevar[varid] = invar.name 
+  ENDFOR
+  namevar = namevar[where(namevar NE '')]
+;------------------------------------------------------------
+; find vargrid for each selected variable...
+;------------------------------------------------------------
+  listgrid = strarr(n_elements(namevar))
+; default definitions
+  IF keyword_set(grid) THEN vargrid = strupcase(grid) ELSE vargrid = 'T'
+; look for values of vargrid for each variable
+  IF finite(glamu[0]) EQ 1 AND NOT keyword_set(grid) THEN BEGIN
+; for each variable, look if we in one of the case corresponding to ROMS conventions?
+    FOR i = 0, n_elements(namevar)-1 do begin
+      invar = ncdf_varinq(cdfid, namevar[i])
+      tmpnm = namedim[invar.dim]
+; are we in one of the case corresponding to ROMS conventions?
+      CASE 1 OF
+        tmpnm[2 <(invar.ndims-1)] EQ 's_w':vargrid = 'W'
+        tmpnm[0] EQ 'xi_rho' AND tmpnm[1] EQ 'eta_rho':listgrid[i] = 'T'
+        tmpnm[0] EQ 'xi_u'   AND tmpnm[1] EQ 'eta_u'  :listgrid[i] = 'U'
+        tmpnm[0] EQ 'xi_v'   AND tmpnm[1] EQ 'eta_v'  :listgrid[i] = 'V'
+        tmpnm[0] EQ 'xi_psi' AND tmpnm[1] EQ 'eta_psi':listgrid[i] = 'F'
+        tmpnm[0] EQ 'xi_rho' AND tmpnm[1] EQ 'eta_v'  :listgrid[i] = 'V'
+        tmpnm[0] EQ 'xi_u'   AND tmpnm[1] EQ 'eta_rho':listgrid[i] = 'U'
+        tmpnm[0] EQ 'xi_u'   AND tmpnm[1] EQ 'eta_v'  :listgrid[i] = 'F'
+        ELSE: 
+      ENDCASE 
+    ENDFOR
+    empty = where(listgrid EQ '')
+    IF empty[0] NE -1 THEN BEGIN   
+; could we define the grid type from the file name??
       pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_']
       gdtype = ['T', 'U', 'V', 'W', 'F']
@@ -104 +203 @@
         ENDFOR
       ENDFOR
-    ENDIF
-  ENDELSE
-;------------------------------------------------------------
-; name of all dimensions
-;------------------------------------------------------------
-  namedim = strarr(infile.ndims)
-  for dimiq = 0, infile.ndims-1 do begin
-    ncdf_diminq, cdfid, dimiq, tmpname, value 
-    namedim[dimiq] = strlowcase(tmpname)
-  ENDFOR
-; we are looking for a x dimension with a name matching one of the following regular expression:
-  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'])
-      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'])
-      stop
-    END
-  ENDCASE
-; we are looking for a y dimension with a name matching one of the following regular expression:
-  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'])
-      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 x dimension'])
-      stop
-    END
-  ENDCASE
-;------------------------------------------------------------
-; name of all variables
-;------------------------------------------------------------
-; we keep only the variables containing at least x, y and time dimension (if existing...)
-  namevar = strarr(infile.nvars)
-  for varid = 0, infile.nvars-1 do begin
-    invar = ncdf_varinq(cdfid, varid) ; what contains the variable?
-    if (where(invar.dim EQ dimidx))[0] NE -1 AND $
-       (where(invar.dim EQ dimidy))[0] NE -1 AND $
-       ((where(invar.dim EQ infile.recdim))[0] NE -1 OR infile.recdim EQ -1) $
-    THEN namevar[varid] = invar.name 
-  ENDFOR
-  namevar = namevar[where(namevar NE '')]
-  listgrid = replicate(vargrid, n_elements(namevar))
+      listgrid[empty] = vargrid
+    ENDIF 
+  ENDIF ELSE listgrid[*] = vargrid
 ;------------------------------------------------------------
 ; time axis
@@ -229 +261 @@
           mots = str_sep(value, ' ')
           unite = mots[0]
-          debut = str_sep(mots[2], '-')
+          err = 0
+          IF unite NE 'seconds' AND unite NE 'hours' AND unite NE 'days' $
+             AND unite NE 'months' AND unite NE 'years' THEN BEGIN
+            dummy = report('time units does not start with seconds/hours/days/months/years')
+            err = 1
+          ENDIF
+          err = err + 1 - stregex(value, '[^ ]* since ([0-9]){4}-([0-9]){2}-([0-9]){2}.*', /boolean)
+          IF err GT 0 THEN BEGIN 
+            dummy = report('attribut units of time has not the good format: [^ ]* since ([0-9]){4}-([0-9]){2}-([0-9]){2}.*')
+            fakecal = 1
+            time = date0fk + lindgen(jpt)
+          ENDIF ELSE BEGIN 
+            debut = str_sep(mots[2], '-')
 ;
 ; now we try to find the attribut called calendar...
@@ -235 +279 @@
 ; If no, we suppose that the calendar is gregorian calendar
 ;
-          if (where(attnames EQ 'calendar'))[0] NE -1 then BEGIN
-            ncdf_attget, cdfid, varid, 'calendar', value
-            value = string(value)
-            CASE value OF
-              'noleap':key_caltype = 'noleap'
-              '360d':key_caltype = '360d'
-              'greg':IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
-              ELSE:BEGIN
+            if (where(attnames EQ 'calendar'))[0] NE -1 then BEGIN
+              ncdf_attget, cdfid, varid, 'calendar', value
+              value = string(value)
+              CASE value OF
+                'noleap':key_caltype = 'noleap'
+                '360d':key_caltype = '360d'
+                'greg':IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
+                ELSE:BEGIN
 ;            notused = report('Unknown calendar: '+value+', we use greg calendar.') 
-                key_caltype = 'greg'
+                  key_caltype = 'greg'
+                END
+              ENDCASE
+            ENDIF ELSE BEGIN
+;        notused = report('Unknown calendar, we use '+key_caltype+' calendar.') 
+              IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
+            ENDELSE
+;
+; BEWARE we have to get back the calendar attribute and ajust time by consequence...
+;
+;
+; We pass time in IDL julian days
+;
+            unite = strlowcase(unite)
+            IF strpos(unite, 's', strlen(unite)-1) NE -1 THEN unite = strmid(unite, 0, strlen(unite)-1)
+            IF strpos(unite, 'julian_') NE -1 THEN unite = strmid(unite, 7)
+            case unite of
+              'second':time = julday(debut[1], debut[2], debut[0])+time/86400.d
+              'hour':time = julday(debut[1], debut[2], debut[0])+time/24.d
+              'day':time = julday(debut[1], debut[2], debut[0])+time
+              'month':BEGIN 
+                if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m
+                   time = julday(debut[1], debut[2], debut[0])+round(time*30) $
+                ELSE for t = 0, n_elements(time)-1 DO $
+                   time[t] = julday(debut[1]+time[t], debut[2], debut[0])
+              END
+              'year':BEGIN
+                if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y
+                   time = julday(debut[1], debut[2], debut[0])+round(time*365) $
+                ELSE for t = 0, n_elements(time)-1 do $
+                   time[t] = julday(debut[1], debut[2], debut[0]+time[t])
               END
             ENDCASE
-          ENDIF ELSE BEGIN
-;        notused = report('Unknown calendar, we use '+key_caltype+' calendar.') 
-            IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
+;
+; high frequency calendar: more than one element per day
+            IF max(histogram([long(time-time[0])])) GT 1 THEN fakecal = 1 ELSE fakecal = 0
+            date0fk = date2jul(19000101)
+            IF keyword_set(fakecal) THEN time = date0fk+lindgen(jpt) $
+            ELSE time = long(time)
+;
           ENDELSE
-;
-; BEWARE we have to get back the calendar attribute and ajust time by consequence...
-;
-;
-; We pass time in IDL julian days
-;
-          unite = strlowcase(unite)
-          IF strpos(unite, 's', strlen(unite)-1) NE -1 THEN unite = strmid(unite, 0, strlen(unite)-1)
-          IF strpos(unite, 'julian_') NE -1 THEN unite = strmid(unite, 7)
-          case unite of
-            'second':time = julday(debut[1], debut[2], debut[0])+time/86400.d
-            'hour':time = julday(debut[1], debut[2], debut[0])+time/24.d
-            'day':time = julday(debut[1], debut[2], debut[0])+time
-            'month':BEGIN 
-              if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m
-                 time = julday(debut[1], debut[2], debut[0])+round(time*30) $
-              ELSE for t = 0, n_elements(time)-1 DO $
-                 time[t] = julday(debut[1]+time[t], debut[2], debut[0])
-            END
-            'year':BEGIN
-              if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y
-                 time = julday(debut[1], debut[2], debut[0])+round(time*365) $
-              ELSE for t = 0, n_elements(time)-1 do $
-                 time[t] = julday(debut[1], debut[2], debut[0]+time[t])
-            END
-          ENDCASE
-;
-; high frequency calendar: more than one element per day
-          IF max(histogram([long(time-time[0])])) GT 1 THEN fakecal = 1 ELSE fakecal = 0
-          date0fk = date2jul(19000101)
-          IF keyword_set(fakecal) THEN time = date0fk+lindgen(jpt) $
-          ELSE time = long(time)
-;
         ENDELSE
       END