source: trunk/SRC/ToBeReviewed/WIDGET/AUTOUR_de_XXX/scanfile.pro @ 219

;+
; @file_comments
;
;
; @categories
;
; 
; @param NAMEFILE
;
; 
; @keyword GRID
;
;
; @keyword _EXTRA
; Used to pass your keywords
; 
; @returns
; 
; 
; @uses
; 
; 
; @restrictions
; 
; 
; @examples
; 
; 
; @history
; 
; 
; @version 
; $Id$
;
; @todo
; seb : I don't know what to do with that...
; 
;-
;
;  liste des presupposes:
;       1) le fichier a lire est un fichier netcdf
;       2) le nom de ce fichier finit
;       par U.nc, V.nc, W.nc, T.nc ou F.nc, la lettre avant le
;       nc designant la grille a laquelle se rapporte la champ.
;       Si tel n''est pas la cas, le fichier est attribue a la grille
;       T.
;       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...','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
;       variable sera prise comme axe des temps. Rq: si plusieurs
;       variables verifient ces criteres on considere la premiere
;       variable
;       5) Cette variable axe des temps doit contenir l''attribut
;       'units'qui doit etre ecrit suivant la syntaxe:
;               "seconds since 0001-01-01 00:00:00" 
;               "hours since 0001-01-01 00:00:00" 
;               "days since 1979-01-01 00:59:59" 
;               "months since 1979-01-01 00:59:59" 
;               "years since 1979-01-01 00:59:59" 
;
; je crois que c''est tout!
;
;        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
;
  compile_opt idl2, strictarrsubs
;
@common
;------------------------------------------------------------
  res = -1
;------------------------------------------------------------
; filename
;------------------------------------------------------------
  fullname = isafile(filename = namefile, IODIRECTORY = iodir, _extra = ex)
;------------------------------------------------------------
; open file
;------------------------------------------------------------
  cdfid = ncdf_open(fullname)
;------------------------------------------------------------
; What contains the file?
;------------------------------------------------------------
  infile = ncdf_inquire(cdfid)  ;
;------------------------------------------------------------
; 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']
      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
      listgrid[empty] = vargrid
    ENDIF 
  ENDIF ELSE listgrid[*] = vargrid
;------------------------------------------------------------
; 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
      IF varid LT infile.nvars THEN BEGIN
        invar = ncdf_varinq(cdfid, varid) 
        varid = varid+1        
      ENDIF ELSE varid = 0
    endrep until (n_elements(invar.dim) EQ 1 AND invar.dim[0] EQ infile.recdim) OR (varid EQ 0)
    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(1>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(1>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 '+invar.name+'!C we create a fake calendar ...')
          fakecal = 1
          time = date0fk + lindgen(1>jpt)
        ENDIF ELSE BEGIN 
; we read the time axis
          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" ;
; time_counter:units = "days since 1979-01-01 00:00:00" ;
; 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]
          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)
          err = 0
          IF unite NE 'second' AND unite NE 'hour' AND unite NE 'day' $
             AND unite NE 'month' AND unite NE 'year' THEN BEGIN
            dummy = report('time units does not start with seconds/hours/days/months/years')
            err = 1
          ENDIF
          IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN
            dummy = report('attribut units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*')
            err = 1
          ENDIF
          IF err GT 0 THEN BEGIN 
            fakecal = 1
            time = date0fk + lindgen(1>jpt)
          ENDIF ELSE BEGIN 
            debut = str_sep(mots[2], '-')
;
; now we try to find the attribut called calendar...
; the the attribute "calendar" exists?
; 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
;            notused = report('Unknown calendar: '+value+', we use greg calendar.') 
                  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
;
            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(1>jpt) $
            ELSE time = long(time)
;
          ENDELSE
        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