source: trunk/SRC/ReadWrite/ncdf_gettime.pro

;+
;
; @file_comments
; get the time axis from a netcdf_file and transforms it in
; Julian days of IDL.
;
; @categories
; Read NetCDF file
;
; @param filename {in}{required}{type=scalar string}
; the name of the ncdf_file
;
; @param cdfid {in}{optional}{type=scalar}
; the ID of the ncdf_file, if the file is already open.
; if not provided, ncdf_gettime open the file defined by filename
;
; @keyword TIMEVAR {type=string}
; It define the name of the variable that
; contains the time axis. This keyword can be useful if there
; is no unlimited dimension or if the time axis selected by default
; (the first 1D array with unlimited dimension) is not the good one.
;
; @keyword CALLER {required}{type=string}
; Used to specify the error messages. Give the name of the calling
; procedure. It can be only '<pro>read_ncdf</pro>' or '<pro>scanfile</pro>'.
;
; @keyword CALTYPE
; Used to specify (or overwrite) the type of calendar that should be
; used. We accept only 'noleap', '360d', 'greg' and 'gregorian'. By
; default, we use the type of calendar defined in the attribute
; calendar. if not, we define it as gregorian.
;
; @keyword ERR
; Set this keyword to a named variable in which the value of the error
; message will be returned
;
; @keyword _EXTRA
; _EXTRA to be able to call ncdf_getmask with _extra keyword
;
; @returns
; a double 1D array of IDL Julian days
; In case of error return -1 if the time dimension was not found
;               or return -jpt if it as been found that the time dimension size is jpt
;
; @restrictions
; the calendar variable must have the units attribute
; following the syntax below:
;
; time_counter:units = "seconds since 0001-01-01 00:00:00" ;
; time_counter:units = "minutes 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" ;
;
; @history
; August 2007: Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;-
FUNCTION ncdf_gettime, filename, cdfid, CALTYPE = caltype $
                       , TIMEVAR = timevar, CALLER = caller, ERR = err, _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
  @cm_4cal                      ;   needed for key_caltype
;
  IF n_elements(cdfid) EQ 0 THEN BEGIN
    cdfid = ncdf_open(isafile(filename, title = 'which file must be open by ncdf_gettime?', IODIRECTORY = iodir, _extra = ex))
    tobeclosed = 1
  ENDIF
  IF n_elements(caller) EQ 0 THEN caller = 'read_ncdf'
  inq = ncdf_inquire(cdfid)
;----------------------------------------------------
;`find the variable containing the time axis
;----------------------------------------------------
; we get its name through the keyword timevar
  IF keyword_set(timevar) THEN BEGIN
    timeid = ncdf_varid(cdfid, timevar)
    IF timeid EQ -1 THEN BEGIN  ;   the variable is not found
      CASE caller OF
        'read_ncdf':err = 'No variable ' + timevar + 'found in '+filename+'. Use the TIMESTEP keyword'
        'scanfile':err = 'No variable ' + timevar + 'found in '+filename+'. We create a fake calendar'
      ENDCASE
      IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
      return, -1
    ENDIF
    timeinq = ncdf_varinq(cdfid, timeid)
    inq.recdim = timeinq.dim[0]
    ncdf_diminq, cdfid, inq.recdim, timedimname, jpt
  ENDIF ELSE BEGIN
; we try to find the time axis automatically
; we look for the infinite dimension
    IF inq.recdim EQ -1 THEN BEGIN
      CASE caller OF
        'read_ncdf':err = 'the file '+filename+' as no infinite dimension. !C Use TIMESTEP or TIMEVAR keyword'
        'scanfile':err = 'the file '+filename+' as no infinite dimension. We create a fake calendar'
      ENDCASE
      IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
      return, -1
    ENDIF
    ncdf_diminq, cdfid, inq.recdim, timedimname, jpt
; we look for the variable containing the time axis
; we look for the first variable having for only dimension inq.recdim
    timeid = 0
    REPEAT BEGIN                           ; As long as we have not find a variable having only one dimension: the infinite one
      timeinq = ncdf_varinq(cdfid, timeid) ; that the variable contain.
      tstok = n_elements(timeinq.dim) EQ 1 AND timeinq.dim[0] EQ inq.recdim
      timeid = timeid+1
    ENDREP UNTIL tstok OR timeid EQ inq.nvars
    IF NOT tstok THEN BEGIN
      CASE caller OF
        'read_ncdf':err = 'the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword'
        'scanfile':err = 'the file '+filename+' has no time axis.!C we create a fake calendar ...'
      ENDCASE
      IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
      return, -jpt
    ENDIF
    timeid = timeid-1
  ENDELSE
;----------------------------------------------------
; look for attribute units and calendar to know how to compte the calendar
;----------------------------------------------------
; no attribute for time variable...
  IF timeinq.natts EQ 0 then begin
    CASE caller OF
      'read_ncdf':err = 'the variable '+timeinq.name+' has no attribute.!C Use the TIMESTEP keyword'
      'scanfile':err = 'the variable '+timeinq.name+' has no attribute.!C we create a fake calendar ...'
    ENDCASE
    IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
    return, -jpt
  ENDIF
; get attributes names
  attnames = strarr(timeinq.natts)
  for attiq = 0, timeinq.natts-1 do attnames[attiq] = strlowcase(ncdf_attname(cdfid, timeid, attiq))
; do we find units attribute?
  IF (where(attnames EQ 'units'))[0] EQ -1 then BEGIN
    CASE caller OF
      'read_ncdf':err = 'Attribute ''units'' not found for the variable '+timeinq.name+' !C Use the TIMESTEP keyword'
      'scanfile':err = 'Attribute ''units'' not found for the variable '+timeinq.name+'!C we create a fake calendar ...'
    ENDCASE
    IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
    return, -jpt
  ENDIF
; Is attribute "calendar" existing?
; If no, we suppose that the calendar is gregorian calendar
;
  ncdf_getatt, cdfid, timeid, calendar = calendar, units = value
  IF keyword_set(caltype) THEN calendar = strlowcase(caltype)
  CASE calendar OF
    'noleap':key_caltype = 'noleap'
    '365day':key_caltype = 'noleap'
    '365days':key_caltype = 'noleap'
    '365_day':key_caltype = 'noleap'
    '365_days':key_caltype = 'noleap'
    '360d':key_caltype = '360d'
    '360day':key_caltype = '360d'
    '360days':key_caltype = '360d'
    '360_day':key_caltype = '360d'
    '360_days':key_caltype = '360d'
    'standard':key_caltype = 'greg'
    'greg':key_caltype = 'greg'
    'gregorian':key_caltype = 'greg'
    ELSE:BEGIN
      dummy = report( ['unknown calendar type: '+calendar, 'we use gregorian calendar'] )
      key_caltype = 'greg'
    END
  ENDCASE
;----------------------------------------------------
; decode units attribute
;----------------------------------------------------
;
  value = strlowcase(value)
  IF value NE 'true julian day' THEN BEGIN
; time_counter:units = "seconds since 0001-01-01 00:00:00" ;
; time_counter:units = "minutes 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" ;
;
; we decript the "units" attribute to find the time origin
    words = str_sep(value, ' ')
    units = words[0]
    IF strpos(units, 's', strlen(units)-1) NE -1 THEN units = strmid(units, 0, strlen(units)-1)
    IF strpos(units, 'julian_') NE -1 THEN units = strmid(units, 7)
    IF units NE 'second' AND units NE 'minute' AND units NE 'hour' AND units NE 'day' $
       AND units NE 'month' AND units NE 'year' THEN BEGIN
      CASE caller OF
        'read_ncdf':err = 'time units does not start with seconds/minutes/hours/days/months/years !C Use the TIMESTEP keyword'
        'scanfile':err = 'time units does not start with seconds/minutes/hours/days/months/years !C we create a fake calendar ...'
      ENDCASE
      IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
      return, -jpt
    ENDIF
    IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN
      CASE caller OF
        'read_ncdf':err = 'attribute units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*!C Use the TIMESTEP keyword'
        'scanfile':err = 'attribute units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*!C we create a fake calendar ...'
      ENDCASE
      IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
      return, -jpt
    ENDIF
    start = str_sep(words[2], '-')
    IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2} ([0-9]){1,2}:([0-9]){1,2}:([0-9]){1,2}', /boolean) EQ 0 THEN BEGIN
      print, 'hh:mm:ss of the calendar origin not found... -> start calendar at 00:00:00'
      start = [start, '0', '0', '0']
    ENDIF ELSE start = [start, str_sep(words[3], ':')]
  ENDIF ELSE units = value
;----------------------------------------------------
; compute time axis
;----------------------------------------------------
  ncdf_varget, cdfid, timeid, time
  time = double(time)
  case units OF
    'true julian day':
    'second':time = julday(start[1], start[2], start[0], start[3], start[4], start[5])+time/86400.d
    'minute':time = julday(start[1], start[2], start[0], start[3], start[4], start[5])+time/1440.d
    'hour':time = julday(start[1], start[2], start[0], start[3], start[4], start[5])+time/24.d
    'day':time = julday(start[1], start[2], start[0], start[3], start[4], start[5])+time
    'month':BEGIN
      if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m
         time = julday(start[1], start[2], start[0], start[3], start[4], start[5]) + round(time*30) $
      ELSE time = julday(start[1]+fix(time), replicate(start[2], jpt), replicate(start[0], jpt), start[3], start[4], start[5])
    END
    'year':BEGIN
      if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y
         time = julday(start[1], start[2], start[0], start[3], start[4], start[5]) + round(time*365) $
      ELSE time = julday(replicate(start[1], jpt), replicate(start[2], jpt), start[0]+fix(time), start[3], start[4], start[5])
    END
  ENDCASE
  time = double(time)
;
  IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
  return, time
END