;+
;
; @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 'read_ncdf' or 'scanfile'.
;
;
; @keyword CALTYPE
; Used to specify (or orverwrite) the type of calendar that should be
; used. We accept only 'noleap', '360d', 'greg' and 'gregorian'. By
; default, we use the type of calendar defied in the attibute
; 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 = "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
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.
timeid = timeid+1
ENDREP UNTIL (n_elements(timeinq.dim) EQ 1 AND timeinq.dim[0] EQ inq.recdim) $
OR timeid EQ inq.nvars+1
IF timeid EQ inq.nvars+1 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 '+fullname+' 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'
'360d':key_caltype = '360d'
'greg':key_caltype = 'greg'
'gregorian':key_caltype = 'greg'
ELSE:BEGIN
dummy = report( ['unknow 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 = "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 '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/hours/days/months/years !C Use the TIMESTEP keyword'
'scanfile':err = 'time units does not start with seconds/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], '-')
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], 0, 0, 0)+time/86400.d
'hour':time = julday(start[1], start[2], start[0], 0, 0, 0)+time/24.d
'day':time = julday(start[1], start[2], start[0], 0, 0, 0)+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], 0, 0, 0) + round(time*30) $
ELSE time = julday(start[1]+fix(time), replicate(start[2], jpt), replicate(start[0], jpt), 0, 0, 0)
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], 0, 0, 0) + round(time*365) $
ELSE time = julday(replicate(start[1], jpt), replicate(start[2], jpt), start[0]+fix(time), 0, 0, 0)
END
ENDCASE
time = double(time)
;
IF keyword_set(tobeclosed) THEN ncdf_close, cdfid
return, time
END