Changeset 69 for trunk/ToBeReviewed/WIDGET/AUTOUR_de_XXX/scanfile.pro

Timestamp:

05/11/06 12:35:53 (18 years ago)

Author:

smasson

Message:

debug + new xxx

File:

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

trunk/ToBeReviewed/WIDGET/AUTOUR_de_XXX/scanfile.pro

@@ -9 +9 @@
 ;       3) ce fichier contient une dimension infinie qui doit etre
 ;       celle qui se rapporte au temps et au mois 2 autres dimensions
-;       dont les noms sont 'x','lon','longitude' et 'y','lat' ou
-;       'latitude' ou bien en majuscule.
+;       dont les noms sont 'x','lon...','xi_...' et 'y','lat...' ou
+;       'eta_...' ou bien en majuscule.
 ;       4) il doit exiter ds ce fichier une unique variable n''ayant
 ;       qu''une dimension et etant la dimension temporelle. cette
@@ -26 +26 @@
 ; je crois que c''est tout!
 ;
-;
-;------------------------------------------------------------
-FUNCTION scanfile, nomficher, _extra = ex
+;        GRID='[UTVWF]' to specify the type of grid. Defaut is (1)
+;        based on the name of the file if the file ends by
+;        GRID[._][TUVFW].NC (not case sensible) or (2) T if case (1)
+;        is not found.
+;
+;------------------------------------------------------------
+FUNCTION scanfile, namefile, GRID = GRID, _extra = ex
 @common
 ;------------------------------------------------------------
-; bidouille pour utiliser les mots cles (on passe par des variables
-; declarees ds un common)
-;------------------------------------------------------------
   res = -1
 ;------------------------------------------------------------
-; choix du nom du fichier
-;------------------------------------------------------------
-  nom = isafile(filename = nomficher, IODIRECTORY = iodir, _extra = ex)
-;------------------------------------------------------------
-; ouverture du fichier nom
-;------------------------------------------------------------
-  cdfid = ncdf_open(nom)
-;------------------------------------------------------------
-; que contient le fichier??
-;------------------------------------------------------------
-  contient = ncdf_inquire(cdfid);
+; filename
+;------------------------------------------------------------
+  fullname = isafile(filename = namefile, IODIRECTORY = iodir, _extra = ex)
+;------------------------------------------------------------
+; open file
+;------------------------------------------------------------
+  cdfid = ncdf_open(fullname)
+;------------------------------------------------------------
+; What contains the file?
+;------------------------------------------------------------
+  infile = ncdf_inquire(cdfid)  ;
 ; find vargrid ...
-  vargrid = 'T'                 ; default definition
-  pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_']
-  gdtype = ['T', 'U', 'V', 'W', 'F']
-  fnametest = strupcase(nom)
-  FOR i = 0, n_elements(pattern)-1 DO BEGIN
-    FOR j = 0, n_elements(gdtype)-1 DO BEGIN
-      substr = pattern[i]+gdtype[j]
-      pos = strpos(fnametest, substr)
-      IF pos NE -1 THEN $
-        vargrid = strmid(fnametest, pos+strlen(substr)-1, 1)
-    ENDFOR
+  IF keyword_set(grid) THEN vargrid = strupcase(grid) ELSE BEGIN
+    vargrid = 'T'               ; default definition
+    IF finite(glamu[0]) EQ 1 THEN BEGIN
+      pattern = ['GRID.', 'GRID_', 'GRID', 'UPID_', '30ID_']
+      gdtype = ['T', 'U', 'V', 'W', 'F']
+      fnametest = strupcase(fullname)
+      FOR i = 0, n_elements(pattern)-1 DO BEGIN
+        FOR j = 0, n_elements(gdtype)-1 DO BEGIN
+          substr = pattern[i]+gdtype[j]
+          pos = strpos(fnametest, substr)
+          IF pos NE -1 THEN $
+             vargrid = strmid(fnametest, pos+strlen(substr)-1, 1)
+        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
-;------------------------------------------------------------
-; nom des differentes variables
-;------------------------------------------------------------
-; on ne garde les noms de variable uniquement si cette variable
-; contient au moins les dimensions qui doivent etre appelles x et y
-; et la dimension infinie
-  nomdim   = strarr(contient.ndims)
-  for dimiq = 0, contient.ndims-1 do begin
-    ncdf_diminq, cdfid, dimiq, name, value ; nom et valeur de la dimension
-    nomdim[dimiq] = strlowcase(name)
-  ENDFOR
-  indexdimx = where(nomdim EQ 'x' OR nomdim EQ 'lon' OR nomdim EQ 'longitude' OR nomdim EQ 'xt_i7_156' OR nomdim EQ 'xi_rho' OR nomdim EQ 'xi_u' OR nomdim EQ 'xi_v' OR nomdim EQ 'xi_psi')
-  indexdimx = indexdimx[0]
-  if indexdimx EQ -1 then begin
-    print, 'one of the dimensions must have the name: ''x'' or ''lon'' or ''longitude'' or ''xt_i7_156'''
+; we are looking for a x dimension...
+  dimidx = where(namedim EQ 'x' OR strmid(namedim, 0, 3) EQ 'lon' OR strmid(namedim, 0, 3) EQ 'xi_' OR namedim EQ 'xt_i7_156')
+  dimidx = dimidx[0]
+  if dimidx EQ -1 then begin
+    print, 'one of the dimensions must have the name: ''x'' or ''lon...'' or ''xi_...'' or ''xt_i7_156'''
     stop
   endif
-  indexdimy = where(nomdim EQ 'y' OR nomdim EQ 'lat' OR nomdim EQ 'latitude' OR nomdim EQ 'yt_j6_75' OR nomdim EQ 'eta_rho' OR nomdim EQ 'eta_u' OR nomdim EQ 'eta_v' OR nomdim EQ 'eta_psi')
-  indexdimy = indexdimy[0]
-  if indexdimy EQ -1 then begin
-    print, 'one of the dimensions must have the name: ''y'' or ''lat'' or ''latitude'' or ''yt_j6_75'''
+; we are looking for a y dimension...
+  dimidy = where(namedim EQ 'y' OR strmid(namedim, 0, 3) EQ 'lat' OR strmid(namedim, 4) EQ 'eta_' OR namedim EQ 'yt_j6_75')
+  dimidy = dimidy[0]
+  if dimidy EQ -1 then begin
+    print, 'one of the dimensions must have the name: ''y'' or ''lat...'' or ''eta_...'' or ''yt_j6_75'''
     stop
   endif
-;
-  namevar = strarr(contient.nvars)
-  for varid = 0, contient.nvars-1 do begin
-    varcontient = ncdf_varinq(cdfid, varid) ; que contient la variable
-    if (where(varcontient.dim EQ indexdimx))[0] NE -1 AND $
-      (where(varcontient.dim EQ indexdimy))[0] NE -1 AND $
-      (where(varcontient.dim EQ contient.recdim))[0] NE -1 THEN namevar[varid] = varcontient.name 
+;------------------------------------------------------------
+; 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))
-;------------------------------------------------------------
-; on recupere l''axe des temps
-;------------------------------------------------------------
-; on recupere le nom de la variable contenant l''axe des temps
-  varid = 0
-  repeat BEGIN ; tant que l''on a pas trouve une variable qui n''a qu''
-                                ; une dimension: la dimension infinie
-    varcontient = ncdf_varinq(cdfid, varid) ; que contient la variable
-    varid = varid+1
-  endrep until n_elements(varcontient.dim) EQ 1 AND varcontient.dim[0] EQ contient.recdim ;
-  varid = varid-1
-;
-; we want to know which attrributes are attached to the time variable...
-;
-  if varcontient.natts EQ 0 then BEGIN
-    ncdf_close, cdfid
-    return, report('the variable '+varcontient.name+' has no attribut.!C we need attribut ''units'' for the calendar ...')
-  endif
-  attnames = strarr(varcontient.natts)
-  for attiq = 0, varcontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, varid, attiq)
-  if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN
-    ncdf_close, cdfid
-    return, report('Attribut ''units'' not found for the variable '+varid.name+'we need it for the calendar...')
-  endif
+  listgrid = replicate(vargrid, n_elements(namevar))
+;------------------------------------------------------------
+; time axis
+;------------------------------------------------------------
+  date0fk = date2jul(19000101)
+  IF infile.recdim EQ -1 THEN BEGIN 
+    jpt = 1
+    time = date0fk
+    fakecal = 1
+  ENDIF ELSE BEGIN
+    ncdf_diminq, cdfid, infile.recdim, timedimname, jpt
+; we look for the variable containing the time axis
+; we look for the first variable having for only dimension infile.recdim
+    varid = 0
+    repeat BEGIN
+      invar = ncdf_varinq(cdfid, varid) 
+      varid = varid+1
+    endrep until n_elements(invar.dim) EQ 1 AND invar.dim[0] EQ infile.recdim
+    varid = varid-1
+;
+    CASE 1 OF
+      varid EQ -1:BEGIN 
+        dummy = report('the file '+fullname+' has no time axis.!C we create a fake calendar ...')
+        fakecal = 1
+        time = date0fk + lindgen(jpt)
+      END
+      invar.natts EQ 0:BEGIN 
+        dummy = report('the variable '+invar.name+' has no attribut.!C we create a fake calendar ...')
+        fakecal = 1
+        time = date0fk + lindgen(jpt)
+      END
+      ELSE:BEGIN
+;
+; we want to know which attributes are attached to the time variable...
+;
+        attnames = strarr(invar.natts)
+        for attiq = 0, invar.natts-1 do attnames[attiq] = ncdf_attname(cdfid, varid, attiq)
+        if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN
+          dummy = report('Attribut ''units'' not found for the variable '+varid.name+'!C we create a fake calendar ...')
+          fakecal = 1
+          time = date0fk + lindgen(jpt)
+        ENDIF ELSE BEGIN 
 ; on lit l''axe des temps
-  ncdf_varget, cdfid, varid, time
-  time = double(time)
-  ncdf_attget, cdfid, varid, 'units', value
+          ncdf_varget, cdfid, varid, time
+          time = double(time)
+          ncdf_attget, cdfid, varid, 'units', value
 ; time_counter:units = "seconds since 0001-01-01 00:00:00" ;
 ; time_counter:units = "hours since 0001-01-01 00:00:00" ;
@@ -126 +151 @@
 ; time_counter:units = "months since 1979-01-01 00:00:00" ;
 ; time_counter:units = "years since 1979-01-01 00:00:00" ;
-  value = strtrim(strcompress(string(value)), 2)
-  mots = str_sep(value, ' ')
-  unite = mots[0]
-  debut = str_sep(mots[2], '-')
+          value = strtrim(strcompress(string(value)), 2)
+          mots = str_sep(value, ' ')
+          unite = mots[0]
+          debut = str_sep(mots[2], '-')
 ;
 ; now we try to find the attribut called calendar...
@@ -135 +160 @@
 ; 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
-          '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'
-          END
-        ENDCASE
-      ENDIF ELSE BEGIN
+                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
+            IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg'
+          ENDELSE
 ;
 ; ATTENTION il faut recuperer l''attribut calendar et ajuster time en
@@ -156 +182 @@
 ;
 ; on passe time en jour julien d''idl 
-; on suppose qu''on utilise le vrai calendrier.
-;
-  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/(long(24)*3600)
-    'hour':time = julday(debut[1], debut[2], debut[0])+time/(long(24))
-    '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
-  time = long(time)
-;------------------------------------------------------------
-  return, {filename:nom, time_counter:time, listvar:namevar, listgrid:strupcase(listgrid), calendartype:key_caltype}
+;
+          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)])) 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
+    ENDCASE
+  ENDELSE
+;------------------------------------------------------------
+  ncdf_close, cdfid
+;------------------------------------------------------------
+  return, {filename:fullname, time_counter:time, listvar:namevar $
+           , listgrid:strupcase(listgrid), caltype:key_caltype $
+           , fakecal:date0fk*fakecal}
 end