NetCDF4 Python

"""
Performs conversions of netCDF time coordinate data to/from datetime objects.
It is based on the original netcdftime.py file from netcdf4 python library.
Changes added by:
Ramiro Checa-Garcia <r.checagarcia@gmail.com>
Description of changes:
The original file of netcdf4 has been changed to include months units. Apparently
given some possible ambiguities on the exact date-time when the units are
"months since" this is not included on the original netcdftime.py code. Here,
I added this possibility for units but the user should be sure about the correct
intrepretation of the results.
"""importnumpyasnpimportmathimportnumpyimportrefromdatetimeimportdatetimeasreal_datetimetry:fromitertoolsimportizipaszipexceptImportError:# python 3.xpassfromdatetimeimportdatetime# RCHG:# Note that this has been changed. To simply replace# the original netcdftime.py this line should be:# from ._datetime import datetime microsec_units=['microseconds','microsecond','microsec','microsecs']millisec_units=['milliseconds','millisecond','millisec','millisecs']sec_units=['second','seconds','sec','secs','s']min_units=['minute','minutes','min','mins']hr_units=['hour','hours','hr','hrs','h']day_units=['day','days','d']month_units=['month','months','mon','mons']_units=microsec_units+millisec_units+sec_units+min_units+hr_units+day_units+month_units_calendars=['standard','gregorian','proleptic_gregorian','noleap','julian','all_leap','365_day','366_day','360_day']__version__='1.4.1'# Adapted from http://delete.me.uk/2005/03/iso8601.htmlISO8601_REGEX=re.compile(r"(?P<year>[+-]?[0-9]{1,4})(-(?P<month>[0-9]{1,2})(-(?P<day>[0-9]{1,2})"r"(((?P<separator1>.)(?P<hour>[0-9]{1,2}):(?P<minute>[0-9]{1,2})(:(?P<second>[0-9]{1,2})(\.(?P<fraction>[0-9]+))?)?)?"r"((?P<separator2>.?)(?P<timezone>Z|(([-+])([0-9]{1,2}):([0-9]{1,2}))))?)?)?)?")TIMEZONE_REGEX=re.compile("(?P<prefix>[+-])(?P<hours>[0-9]{1,2}):(?P<minutes>[0-9]{1,2})")defJulianDayFromDate(date,calendar='standard'):"""
creates a Julian Day from a 'datetime-like' object. Returns the fractional
Julian Day (resolution approx 0.1 second).
if calendar='standard' or 'gregorian' (default), Julian day follows Julian
Calendar on and before 1582-10-5, Gregorian calendar after 1582-10-15.
if calendar='proleptic_gregorian', Julian Day follows gregorian calendar.
if calendar='julian', Julian Day follows julian calendar.
Algorithm:
Meeus, Jean (1998) Astronomical Algorithms (2nd Edition). Willmann-Bell,
Virginia. p. 63
"""# based on redate.py by David Finlayson.# check if input was scalar and change return accordinglyisscalar=Falsetry:date[0]except:isscalar=Truedate=np.atleast_1d(np.array(date))year=np.empty(len(date),dtype=np.int32)month=year.copy()day=year.copy()hour=year.copy()minute=year.copy()second=year.copy()microsecond=year.copy()fori,dinenumerate(date):year[i]=d.yearmonth[i]=d.monthday[i]=d.dayhour[i]=d.hourminute[i]=d.minutesecond[i]=d.secondmicrosecond[i]=d.microsecond# Convert time to fractions of a dayday=day+hour/24.0+minute/1440.0+(second+microsecond/1.e6)/86400.0# Start Meeus algorithm (variables are in his notation)month_lt_3=month<3month[month_lt_3]=month[month_lt_3]+12year[month_lt_3]=year[month_lt_3]-1# MC - assure array# A = np.int64(year / 100)A=(year/100).astype(np.int64)# MC# jd = int(365.25 * (year + 4716)) + int(30.6001 * (month + 1)) + \# day - 1524.5jd=365.*year+np.int32(0.25*year+2000.)+np.int32(30.6001*(month+1))+ \
day+1718994.5# optionally adjust the jd for the switch from# the Julian to Gregorian Calendar# here assumed to have occurred the day after 1582 October 4ifcalendarin['standard','gregorian']:# MC - do not have to be contiguous dates# if np.min(jd) >= 2299170.5:# # 1582 October 15 (Gregorian Calendar)# B = 2 - A + np.int32(A / 4)# elif np.max(jd) < 2299160.5:# # 1582 October 5 (Julian Calendar)# B = np.zeros(len(jd))# else:# print(date, calendar, jd)# raise ValueError(# 'impossible date (falls in gap between end of Julian calendar and beginning of Gregorian calendar')ifnp.any((jd>=2299160.5)&(jd<2299170.5)):# missing days in Gregorian calendarraiseValueError('impossible date (falls in gap between end of Julian calendar and beginning of Gregorian calendar')B=np.zeros(len(jd))# 1582 October 5 (Julian Calendar)ii=np.where(jd>=2299170.5)[0]# 1582 October 15 (Gregorian Calendar)ifii.size>0:B[ii]=2-A[ii]+np.int32(A[ii]/4)elifcalendar=='proleptic_gregorian':B=2-A+np.int32(A/4)elifcalendar=='julian':B=np.zeros(len(jd))else:raiseValueError('unknown calendar, must be one of julian,standard,gregorian,proleptic_gregorian, got %s'%calendar)# adjust for Julian calendar if necessaryjd=jd+B# Add a small offset (proportional to Julian date) for correct re-conversion.# This is about 45 microseconds in 2000 for Julian date starting -4712.# (pull request #433).eps=np.finfo(float).epseps=np.maximum(eps*jd,eps)jd+=epsifisscalar:returnjd[0]else:returnjddef_NoLeapDayFromDate(date):"""
creates a Julian Day for a calendar with no leap years from a datetime
instance. Returns the fractional Julian Day (resolution approx 0.1 second).
"""year=date.yearmonth=date.monthday=date.dayhour=date.hourminute=date.minutesecond=date.secondmicrosecond=date.microsecond# Convert time to fractions of a dayday=day+hour/24.0+minute/1440.0+(second+microsecond/1.e6)/86400.0# Start Meeus algorithm (variables are in his notation)if(month<3):month=month+12year=year-1jd=int(365.*(year+4716))+int(30.6001*(month+1))+ \
day-1524.5returnjddef_AllLeapFromDate(date):"""
creates a Julian Day for a calendar where all years have 366 days from
a 'datetime-like' object.
Returns the fractional Julian Day (resolution approx 0.1 second).
"""year=date.yearmonth=date.monthday=date.dayhour=date.hourminute=date.minutesecond=date.secondmicrosecond=date.microsecond# Convert time to fractions of a dayday=day+hour/24.0+minute/1440.0+(second+microsecond/1.e6)/86400.0# Start Meeus algorithm (variables are in his notation)if(month<3):month=month+12year=year-1jd=int(366.*(year+4716))+int(30.6001*(month+1))+ \
day-1524.5returnjddef_360DayFromDate(date):"""
creates a Julian Day for a calendar where all months have 30 daysfrom
a 'datetime-like' object.
Returns the fractional Julian Day (resolution approx 0.1 second).
"""year=date.yearmonth=date.monthday=date.dayhour=date.hourminute=date.minutesecond=date.secondmicrosecond=date.microsecond# Convert time to fractions of a dayday=day+hour/24.0+minute/1440.0+(second+microsecond/1.e6)/86400.0jd=int(360.*(year+4716))+int(30.*(month-1))+dayreturnjddefDateFromJulianDay(JD,calendar='standard'):"""
returns a 'datetime-like' object given Julian Day. Julian Day is a
fractional day with a resolution of approximately 0.1 seconds.
if calendar='standard' or 'gregorian' (default), Julian day follows Julian
Calendar on and before 1582-10-5, Gregorian calendar after 1582-10-15.
if calendar='proleptic_gregorian', Julian Day follows gregorian calendar.
if calendar='julian', Julian Day follows julian calendar.
The datetime object is a 'real' datetime object if the date falls in
the Gregorian calendar (i.e. calendar='proleptic_gregorian', or
calendar = 'standard'/'gregorian' and the date is after 1582-10-15).
Otherwise, it's a 'phony' datetime object which is actually an instance
of netcdftime.datetime.
Algorithm:
Meeus, Jean (1998) Astronomical Algorithms (2nd Edition). Willmann-Bell,
Virginia. p. 63
"""# based on redate.py by David Finlayson.julian=np.array(JD,dtype=float)ifnp.min(julian)<0:raiseValueError('Julian Day must be positive')dayofwk=np.atleast_1d(np.int32(np.fmod(np.int32(julian+1.5),7)))# get the day (Z) and the fraction of the day (F)# add 0.000005 which is 452 ms in case of jd being after# second 23:59:59 of a day we want to round to the next day see issue #75Z=np.atleast_1d(np.int32(np.round(julian)))F=np.atleast_1d(julian+0.5-Z).astype(np.float64)ifcalendarin['standard','gregorian']:# MC# alpha = int((Z - 1867216.25)/36524.25)# A = Z + 1 + alpha - int(alpha/4)alpha=np.int32(((Z-1867216.)-0.25)/36524.25)A=Z+1+alpha-np.int32(0.25*alpha)# check if dates before oct 5th 1582 are in the arrayind_before=np.where(julian<2299160.5)[0]iflen(ind_before)>0:A[ind_before]=Z[ind_before]elifcalendar=='proleptic_gregorian':# MC# alpha = int((Z - 1867216.25)/36524.25)# A = Z + 1 + alpha - int(alpha/4)alpha=np.int32(((Z-1867216.)-0.25)/36524.25)A=Z+1+alpha-np.int32(0.25*alpha)elifcalendar=='julian':A=Zelse:raiseValueError('unknown calendar, must be one of julian,standard,gregorian,proleptic_gregorian, got %s'%calendar)B=A+1524# MC# C = int((B - 122.1)/365.25)# D = int(365.25 * C)C=np.atleast_1d(np.int32(6680.+((B-2439870.)-122.1)/365.25))D=np.atleast_1d(np.int32(365*C+np.int32(0.25*C)))E=np.atleast_1d(np.int32((B-D)/30.6001))# Convert to dateday=np.clip(B-D-np.int64(30.6001*E)+F,1,None)nday=B-D-123dayofyr=nday-305ind_nday_before=np.where(nday<=305)[0]iflen(ind_nday_before)>0:dayofyr[ind_nday_before]=nday[ind_nday_before]+60# MC# if E < 14:# month = E - 1# else:# month = E - 13# if month > 2:# year = C - 4716# else:# year = C - 4715month=E-1month[month>12]=month[month>12]-12year=C-4715year[month>2]=year[month>2]-1year[year<=0]=year[year<=0]-1# a leap year?leap=np.zeros(len(year),dtype=dayofyr.dtype)leap[year%4==0]=1ifcalendar=='proleptic_gregorian':leap[(year%100==0)&(year%400!=0)]=0elifcalendarin['standard','gregorian']:leap[(year%100==0)&(year%400!=0)&(julian<2299160.5)]=0inc_idx=np.where((leap==1)&(month>2))[0]dayofyr[inc_idx]=dayofyr[inc_idx]+leap[inc_idx]# Subtract the offset from JulianDayFromDate from the microseconds (pull# request #433).eps=np.finfo(float).epseps=np.maximum(eps*julian,eps)hour=np.clip((F*24.).astype(np.int64),0,23)F-=hour/24.minute=np.clip((F*1440.).astype(np.int64),0,59)# this is an overestimation due to added offset in JulianDayFromDatesecond=np.clip((F-minute/1440.)*86400.,0,None)microsecond=(second%1)*1.e6# remove the offset from the microsecond calculation.microsecond=np.clip(microsecond-eps*86400.*1e6,0,999999)# convert year, month, day, hour, minute, second to int32year=year.astype(np.int32)month=month.astype(np.int32)day=day.astype(np.int32)hour=hour.astype(np.int32)minute=minute.astype(np.int32)second=second.astype(np.int32)microsecond=microsecond.astype(np.int32)# check if input was scalar and change return accordinglyisscalar=Falsetry:JD[0]except:isscalar=True# return a 'real' datetime instance if calendar is gregorian.ifcalendarin'proleptic_gregorian'or \
(calendarin['standard','gregorian']andlen(ind_before)==0):ifnotisscalar:returnnp.array([real_datetime(*args)forargsinzip(year,month,day,hour,minute,second,microsecond)])else:returnreal_datetime(year[0],month[0],day[0],hour[0],minute[0],second[0],microsecond[0])else:# or else, return a 'datetime-like' instance.ifnotisscalar:returnnp.array([datetime(*args)forargsinzip(year,month,day,hour,minute,second,microsecond,dayofwk,dayofyr)])else:returndatetime(year[0],month[0],day[0],hour[0],minute[0],second[0],microsecond[0],dayofwk[0],dayofyr[0])def_DateFromNoLeapDay(JD):"""
returns a 'datetime-like' object given Julian Day for a calendar with no leap
days. Julian Day is a fractional day with a resolution of approximately 0.1 seconds.
"""# based on redate.py by David Finlayson.ifJD<0:raiseValueError('Julian Day must be positive')dayofwk=int(math.fmod(int(JD+1.5),7))(F,Z)=math.modf(JD+0.5)Z=int(Z)A=ZB=A+1524C=int((B-122.1)/365.)D=int(365.*C)E=int((B-D)/30.6001)# Convert to dateday=B-D-int(30.6001*E)+Fnday=B-D-123ifnday<=305:dayofyr=nday+60else:dayofyr=nday-305ifE<14:month=E-1else:month=E-13ifmonth>2:year=C-4716else:year=C-4715# Convert fractions of a day to time(dfrac,days)=math.modf(day/1.0)(hfrac,hours)=math.modf(dfrac*24.0)(mfrac,minutes)=math.modf(hfrac*60.0)(sfrac,seconds)=math.modf(mfrac*60.0)microseconds=sfrac*1.e6returndatetime(year,month,int(days),int(hours),int(minutes),int(seconds),int(microseconds),dayofwk,dayofyr)def_DateFromAllLeap(JD):"""
returns a 'datetime-like' object given Julian Day for a calendar where all
years have 366 days.
Julian Day is a fractional day with a resolution of approximately 0.1 seconds.
"""# based on redate.py by David Finlayson.ifJD<0:raiseValueError('Julian Day must be positive')dayofwk=int(math.fmod(int(JD+1.5),7))(F,Z)=math.modf(JD+0.5)Z=int(Z)A=ZB=A+1524C=int((B-122.1)/366.)D=int(366.*C)E=int((B-D)/30.6001)# Convert to dateday=B-D-int(30.6001*E)+Fnday=B-D-123ifnday<=305:dayofyr=nday+60else:dayofyr=nday-305ifE<14:month=E-1else:month=E-13ifmonth>2:dayofyr=dayofyr+1ifmonth>2:year=C-4716else:year=C-4715# Convert fractions of a day to time(dfrac,days)=math.modf(day/1.0)(hfrac,hours)=math.modf(dfrac*24.0)(mfrac,minutes)=math.modf(hfrac*60.0)(sfrac,seconds)=math.modf(mfrac*60.0)microseconds=sfrac*1.e6returndatetime(year,month,int(days),int(hours),int(minutes),int(seconds),int(microseconds),dayofwk,dayofyr)def_DateFrom360Day(JD):"""
returns a 'datetime-like' object given Julian Day for a calendar where all
months have 30 days.
Julian Day is a fractional day with a resolution of approximately 0.1 seconds.
"""ifJD<0:raiseValueError('Julian Day must be positive')#jd = int(360. * (year + 4716)) + int(30. * (month - 1)) + day(F,Z)=math.modf(JD)year=int((Z-0.5)/360.)-4716dayofyr=Z-(year+4716)*360month=int((dayofyr-0.5)/30)+1day=dayofyr-(month-1)*30+F# Convert fractions of a day to time(dfrac,days)=math.modf(day/1.0)(hfrac,hours)=math.modf(dfrac*24.0)(mfrac,minutes)=math.modf(hfrac*60.0)(sfrac,seconds)=math.modf(mfrac*60.0)microseconds=sfrac*1.e6returndatetime(year,month,int(days),int(hours),int(minutes),int(seconds),int(microseconds),-1,dayofyr)def_dateparse(timestr):"""parse a string of the form time-units since yyyy-mm-dd hh:mm:ss
return a tuple (units,utc_offset, datetimeinstance)"""timestr_split=timestr.split()units=timestr_split[0].lower()ifunitsnotin_units:raiseValueError("units must be one of 'seconds', 'minutes', 'hours' or 'days' (or singular version of these), got '%s'"%units)iftimestr_split[1].lower()!='since':raiseValueError("no 'since' in unit_string")# parse the date string.n=timestr.find('since')+6year,month,day,hour,minute,second,utc_offset=_parse_date(timestr[n:].strip())returnunits,utc_offset,datetime(year,month,day,hour,minute,second)classutime:"""
Performs conversions of netCDF time coordinate
data to/from datetime objects.
To initialize: C{t = utime(unit_string,calendar='standard')}
where
B{C{unit_string}} is a string of the form
C{'time-units since <time-origin>'} defining the time units.
Valid time-units are days, hours, minutes and seconds (the singular forms
are also accepted). An example unit_string would be C{'hours
since 0001-01-01 00:00:00'}.
The B{C{calendar}} keyword describes the calendar used in the time calculations.
All the values currently defined in the U{CF metadata convention
<http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.1/cf-conventions.html#time-coordinate>}
are accepted. The default is C{'standard'}, which corresponds to the mixed
Gregorian/Julian calendar used by the C{udunits library}. Valid calendars
are:
C{'gregorian'} or C{'standard'} (default):
Mixed Gregorian/Julian calendar as defined by udunits.
C{'proleptic_gregorian'}:
A Gregorian calendar extended to dates before 1582-10-15. That is, a year
is a leap year if either (i) it is divisible by 4 but not by 100 or (ii)
it is divisible by 400.
C{'noleap'} or C{'365_day'}:
Gregorian calendar without leap years, i.e., all years are 365 days long.
all_leap or 366_day Gregorian calendar with every year being a leap year,
i.e., all years are 366 days long.
C{'360_day'}:
All years are 360 days divided into 30 day months.
C{'julian'}:
Proleptic Julian calendar, extended to dates after 1582-10-5. A year is a
leap year if it is divisible by 4.
The C{L{num2date}} and C{L{date2num}} class methods can used to convert datetime
instances to/from the specified time units using the specified calendar.
The datetime instances returned by C{num2date} are 'real' python datetime
objects if the date falls in the Gregorian calendar (i.e.
C{calendar='proleptic_gregorian', 'standard'} or C{'gregorian'} and
the date is after 1582-10-15). Otherwise, they are 'phony' datetime
objects which are actually instances of C{L{netcdftime.datetime}}. This is
because the python datetime module cannot handle the weird dates in some
calendars (such as C{'360_day'} and C{'all_leap'}) which don't exist in any real
world calendar.
Example usage:
>>> from netcdftime import utime
>>> from datetime import datetime
>>> cdftime = utime('hours since 0001-01-01 00:00:00')
>>> date = datetime.now()
>>> print date
2006-03-17 16:04:02.561678
>>>
>>> t = cdftime.date2num(date)
>>> print t
17577328.0672
>>>
>>> date = cdftime.num2date(t)
>>> print date
2006-03-17 16:04:02
>>>
The resolution of the transformation operation is approximately 0.1 seconds.
Warning: Dates between 1582-10-5 and 1582-10-15 do not exist in the
C{'standard'} or C{'gregorian'} calendars. An exception will be raised if you pass
a 'datetime-like' object in that range to the C{L{date2num}} class method.
Words of Wisdom from the British MetOffice concerning reference dates:
"udunits implements the mixed Gregorian/Julian calendar system, as
followed in England, in which dates prior to 1582-10-15 are assumed to use
the Julian calendar. Other software cannot be relied upon to handle the
change of calendar in the same way, so for robustness it is recommended
that the reference date be later than 1582. If earlier dates must be used,
it should be noted that udunits treats 0 AD as identical to 1 AD."
@ivar origin: datetime instance defining the origin of the netCDF time variable.
@ivar calendar: the calendar used (as specified by the C{calendar} keyword).
@ivar unit_string: a string defining the the netCDF time variable.
@ivar units: the units part of C{unit_string} (i.e. 'days', 'hours', 'seconds').
"""def__init__(self,unit_string,calendar='standard'):"""
@param unit_string: a string of the form
C{'time-units since <time-origin>'} defining the time units.
Valid time-units are days, hours, minutes and seconds (the singular forms
are also accepted). An example unit_string would be C{'hours
since 0001-01-01 00:00:00'}.
@keyword calendar: describes the calendar used in the time calculations.
All the values currently defined in the U{CF metadata convention
<http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.1/cf-conventions.html#time-coordinate>}
are accepted. The default is C{'standard'}, which corresponds to the mixed
Gregorian/Julian calendar used by the C{udunits library}. Valid calendars
are:
- C{'gregorian'} or C{'standard'} (default):
Mixed Gregorian/Julian calendar as defined by udunits.
- C{'proleptic_gregorian'}:
A Gregorian calendar extended to dates before 1582-10-15. That is, a year
is a leap year if either (i) it is divisible by 4 but not by 100 or (ii)
it is divisible by 400.
- C{'noleap'} or C{'365_day'}:
Gregorian calendar without leap years, i.e., all years are 365 days long.
- C{'all_leap'} or C{'366_day'}:
Gregorian calendar with every year being a leap year, i.e.,
all years are 366 days long.
-C{'360_day'}:
All years are 360 days divided into 30 day months.
-C{'julian'}:
Proleptic Julian calendar, extended to dates after 1582-10-5. A year is a
leap year if it is divisible by 4.
@returns: A class instance which may be used for converting times from netCDF
units to datetime objects.
"""calendar=calendar.lower()ifcalendarin_calendars:self.calendar=calendarelse:raiseValueError("calendar must be one of %s, got '%s'"%(str(_calendars),calendar))units,tzoffset,self.origin=_dateparse(unit_string)# real-world calendars limited to positive reference years.ifself.calendarin['julian','standard','gregorian','proleptic_gregorian']:ifself.origin.year==0:msg='zero not allowed as a reference year, does not exist in Julian or Gregorian calendars'raiseValueError(msg)elifself.origin.year<0:msg='negative reference year in time units, must be >= 1'raiseValueError(msg)self.tzoffset=tzoffset# time zone offset in minutesself.units=unitsself.unit_string=unit_stringifself.calendarin['noleap','365_day']andself.origin.month==2andself.origin.day==29:raiseValueError('cannot specify a leap day as the reference time with the noleap calendar')ifself.calendar=='360_day'andself.origin.day>30:raiseValueError('there are only 30 days in every month with the 360_day calendar')ifself.calendarin['noleap','365_day']:self._jd0=_NoLeapDayFromDate(self.origin)elifself.calendarin['all_leap','366_day']:self._jd0=_AllLeapFromDate(self.origin)elifself.calendar=='360_day':self._jd0=_360DayFromDate(self.origin)else:self._jd0=JulianDayFromDate(self.origin,calendar=self.calendar)defdate2num(self,date):"""
Returns C{time_value} in units described by L{unit_string}, using
the specified L{calendar}, given a 'datetime-like' object.
The datetime object must represent UTC with no time-zone offset.
If there is a time-zone offset implied by L{unit_string}, it will
be applied to the returned numeric values.
Resolution is approximately 0.1 seconds.
If C{calendar = 'standard'} or C{'gregorian'} (indicating
that the mixed Julian/Gregorian calendar is to be used), an
exception will be raised if the 'datetime-like' object describes
a date between 1582-10-5 and 1582-10-15.
Works for scalars, sequences and numpy arrays.
Returns a scalar if input is a scalar, else returns a numpy array.
"""isscalar=Falsetry:date[0]except:isscalar=Trueifnotisscalar:date=numpy.array(date)shape=date.shapeifself.calendarin['julian','standard','gregorian','proleptic_gregorian']:ifisscalar:jdelta=JulianDayFromDate(date,self.calendar)-self._jd0else:jdelta=JulianDayFromDate(date.flat,self.calendar)-self._jd0elifself.calendarin['noleap','365_day']:ifisscalar:ifdate.month==2anddate.day==29:raiseValueError('there is no leap day in the noleap calendar')jdelta=_NoLeapDayFromDate(date)-self._jd0else:jdelta=[]fordindate.flat:ifd.month==2andd.day==29:raiseValueError('there is no leap day in the noleap calendar')jdelta.append(_NoLeapDayFromDate(d)-self._jd0)elifself.calendarin['all_leap','366_day']:ifisscalar:jdelta=_AllLeapFromDate(date)-self._jd0else:jdelta=[_AllLeapFromDate(d)-self._jd0fordindate.flat]elifself.calendar=='360_day':ifisscalar:ifdate.day>30:raiseValueError('there are only 30 days in every month with the 360_day calendar')jdelta=_360DayFromDate(date)-self._jd0else:jdelta=[]fordindate.flat:ifd.day>30:raiseValueError('there are only 30 days in every month with the 360_day calendar')jdelta.append(_360DayFromDate(d)-self._jd0)ifnotisscalar:jdelta=numpy.array(jdelta)# convert to desired units, subtract time zone offset.ifself.unitsinmicrosec_units:jdelta=jdelta*86400.*1.e6-self.tzoffset*60.*1.e6elifself.unitsinmillisec_units:jdelta=jdelta*86400.*1.e3-self.tzoffset*60.*1.e3elifself.unitsinsec_units:jdelta=jdelta*86400.-self.tzoffset*60.elifself.unitsinmin_units:jdelta=jdelta*1440.-self.tzoffsetelifself.unitsinhr_units:jdelta=jdelta*24.-self.tzoffset/60.elifself.unitsinday_units:jdelta=jdelta-self.tzoffset/1440elifself.unitsinmonth_units:jdelta=jdelta*12/365.242198781-self.tzoffset/1440.## Totally uncertain!!!!!!!!!!!!!!!!!!!!!!!!!!!!!else:raiseValueError('unsupported time units')ifisscalar:returnjdeltaelse:returnnumpy.reshape(jdelta,shape)defnum2date(self,time_value):"""
Return a 'datetime-like' object given a C{time_value} in units
described by L{unit_string}, using L{calendar}.
dates are in UTC with no offset, even if L{unit_string} contains
a time zone offset from UTC.
Resolution is approximately 0.1 seconds.
Works for scalars, sequences and numpy arrays.
Returns a scalar if input is a scalar, else returns a numpy array.
The datetime instances returned by C{num2date} are 'real' python datetime
objects if the date falls in the Gregorian calendar (i.e.
C{calendar='proleptic_gregorian'}, or C{calendar = 'standard'/'gregorian'} and
the date is after 1582-10-15). Otherwise, they are 'phony' datetime
objects which are actually instances of netcdftime.datetime. This is
because the python datetime module cannot handle the weird dates in some
calendars (such as C{'360_day'} and C{'all_leap'}) which
do not exist in any real world calendar.
"""isscalar=Falsetry:time_value[0]except:isscalar=Trueismasked=Falseifhasattr(time_value,'mask'):mask=time_value.maskismasked=Trueifnotisscalar:time_value=numpy.array(time_value,dtype='d')shape=time_value.shape# convert to desired units, add time zone offset.ifself.unitsinmicrosec_units:jdelta=time_value/86400000000.+self.tzoffset/1440.elifself.unitsinmillisec_units:jdelta=time_value/86400000.+self.tzoffset/1440.elifself.unitsinsec_units:jdelta=time_value/86400.+self.tzoffset/1440.elifself.unitsinmin_units:jdelta=time_value/1440.+self.tzoffset/1440.elifself.unitsinhr_units:jdelta=time_value/24.+self.tzoffset/1440.elifself.unitsinday_units:jdelta=time_value+self.tzoffset/1440.elifself.unitsinmonth_units:jdelta=time_value*365.242198781/12+self.tzoffset/1440.# Months is manually added, please take care on their use. else:raiseValueError('unsupported time units')jd=self._jd0+jdeltaifself.calendarin['julian','standard','gregorian','proleptic_gregorian']:ifnotisscalar:ifismasked:date=[]forj,minzip(jd.flat,mask.flat):ifnotm:date.append(DateFromJulianDay(j,self.calendar))else:date.append(None)else:date=DateFromJulianDay(jd.flat,self.calendar)else:ifismaskedandmask.item():date=Noneelse:date=DateFromJulianDay(jd,self.calendar)elifself.calendarin['noleap','365_day']:ifnotisscalar:date=[_DateFromNoLeapDay(j)forjinjd.flat]else:date=_DateFromNoLeapDay(jd)elifself.calendarin['all_leap','366_day']:ifnotisscalar:date=[_DateFromAllLeap(j)forjinjd.flat]else:date=_DateFromAllLeap(jd)elifself.calendar=='360_day':ifnotisscalar:date=[_DateFrom360Day(j)forjinjd.flat]else:date=_DateFrom360Day(jd)ifisscalar:returndateelse:returnnumpy.reshape(numpy.array(date),shape)def_parse_timezone(tzstring):"""Parses ISO 8601 time zone specs into tzinfo offsets
Adapted from pyiso8601 (http://code.google.com/p/pyiso8601/)
"""iftzstring=="Z":return0# This isn't strictly correct, but it's common to encounter dates without# timezones so I'll assume the default (which defaults to UTC).iftzstringisNone:return0m=TIMEZONE_REGEX.match(tzstring)prefix,hours,minutes=m.groups()hours,minutes=int(hours),int(minutes)ifprefix=="-":hours=-hoursminutes=-minutesreturnminutes+hours*60.def_parse_date(datestring):"""Parses ISO 8601 dates into datetime objects
The timezone is parsed from the date string, assuming UTC
by default.
Adapted from pyiso8601 (http://code.google.com/p/pyiso8601/)
"""ifnotisinstance(datestring,str)andnotisinstance(datestring,unicode):raiseValueError("Expecting a string %r"%datestring)m=ISO8601_REGEX.match(datestring.strip())ifnotm:raiseValueError("Unable to parse date string %r"%datestring)groups=m.groupdict()tzoffset_mins=_parse_timezone(groups["timezone"])ifgroups["hour"]isNone:groups["hour"]=0ifgroups["minute"]isNone:groups["minute"]=0ifgroups["second"]isNone:groups["second"]=0# if groups["fraction"] is None:# groups["fraction"] = 0# else:# groups["fraction"] = int(float("0.%s" % groups["fraction"]) * 1e6)iyear=int(groups["year"])returniyear,int(groups["month"]),int(groups["day"]),\
int(groups["hour"]),int(groups["minute"]),int(groups["second"]),\
tzoffset_minsdef_check_index(indices,times,nctime,calendar,select):"""Return True if the time indices given correspond to the given times,
False otherwise.
Parameters:
indices : sequence of integers
Positive integers indexing the time variable.
times : sequence of times.
Reference times.
nctime : netCDF Variable object
NetCDF time object.
calendar : string
Calendar of nctime.
select : string
Index selection method.
"""N=nctime.shape[0]if(indices<0).any():returnFalseif(indices>=N).any():returnFalsetry:t=nctime[indices]nctime=nctime# WORKAROUND TO CHANGES IN SLICING BEHAVIOUR in 1.1.2# this may be unacceptably slow...# if indices are unsorted, or there are duplicate# values in indices, read entire time variable into numpy# array so numpy slicing rules can be used.exceptIndexError:nctime=nctime[:]t=nctime[indices]# if fancy indexing not available, fall back on this.# t=[]# for ind in indices:# t.append(nctime[ind])ifselect=='exact':returnnumpy.all(t==times)elifselect=='before':ta=nctime[numpy.clip(indices+1,0,N-1)]returnnumpy.all(t<=times)andnumpy.all(ta>times)elifselect=='after':tb=nctime[numpy.clip(indices-1,0,N-1)]returnnumpy.all(t>=times)andnumpy.all(tb<times)elifselect=='nearest':ta=nctime[numpy.clip(indices+1,0,N-1)]tb=nctime[numpy.clip(indices-1,0,N-1)]delta_after=ta-tdelta_before=t-tbdelta_check=numpy.abs(times-t)returnnumpy.all(delta_check<=delta_after)andnumpy.all(delta_check<=delta_before)defdate2index(dates,nctime,calendar=None,select='exact'):"""
date2index(dates, nctime, calendar=None, select='exact')
Return indices of a netCDF time variable corresponding to the given dates.
@param dates: A datetime object or a sequence of datetime objects.
The datetime objects should not include a time-zone offset.
@param nctime: A netCDF time variable object. The nctime object must have a
C{units} attribute. The entries are assumed to be stored in increasing
order.
@param calendar: Describes the calendar used in the time calculation.
Valid calendars C{'standard', 'gregorian', 'proleptic_gregorian'
'noleap', '365_day', '360_day', 'julian', 'all_leap', '366_day'}.
Default is C{'standard'}, which is a mixed Julian/Gregorian calendar
If C{calendar} is None, its value is given by C{nctime.calendar} or
C{standard} if no such attribute exists.
@param select: C{'exact', 'before', 'after', 'nearest'}
The index selection method. C{exact} will return the indices perfectly
matching the dates given. C{before} and C{after} will return the indices
corresponding to the dates just before or just after the given dates if
an exact match cannot be found. C{nearest} will return the indices that
correpond to the closest dates.
"""try:nctime.unitsexceptAttributeError:raiseAttributeError("netcdf time variable is missing a 'units' attribute")# Setting the calendar.ifcalendar==None:calendar=getattr(nctime,'calendar','standard')cdftime=utime(nctime.units,calendar=calendar)times=cdftime.date2num(dates)returntime2index(times,nctime,calendar=calendar,select=select)deftime2index(times,nctime,calendar=None,select='exact'):"""
time2index(times, nctime, calendar=None, select='exact')
Return indices of a netCDF time variable corresponding to the given times.
@param times: A numeric time or a sequence of numeric times.
@param nctime: A netCDF time variable object. The nctime object must have a
C{units} attribute. The entries are assumed to be stored in increasing
order.
@param calendar: Describes the calendar used in the time calculation.
Valid calendars C{'standard', 'gregorian', 'proleptic_gregorian'
'noleap', '365_day', '360_day', 'julian', 'all_leap', '366_day'}.
Default is C{'standard'}, which is a mixed Julian/Gregorian calendar
If C{calendar} is None, its value is given by C{nctime.calendar} or
C{standard} if no such attribute exists.
@param select: C{'exact', 'before', 'after', 'nearest'}
The index selection method. C{exact} will return the indices perfectly
matching the times given. C{before} and C{after} will return the indices
corresponding to the times just before or just after the given times if
an exact match cannot be found. C{nearest} will return the indices that
correpond to the closest times.
"""try:nctime.unitsexceptAttributeError:raiseAttributeError("netcdf time variable is missing a 'units' attribute")# Setting the calendar.ifcalendar==None:calendar=getattr(nctime,'calendar','standard')num=numpy.atleast_1d(times)N=len(nctime)# Trying to infer the correct index from the starting time and the stride.# This assumes that the times are increasing uniformly.iflen(nctime)>=2:t0,t1=nctime[:2]dt=t1-t0else:t0=nctime[0]dt=1.ifselectin['exact','before']:index=numpy.array((num-t0)/dt,int)elifselect=='after':index=numpy.array(numpy.ceil((num-t0)/dt),int)else:index=numpy.array(numpy.around((num-t0)/dt),int)# Checking that the index really corresponds to the given time.# If the times do not correspond, then it means that the times# are not increasing uniformly and we try the bisection method.ifnot_check_index(index,times,nctime,calendar,select):# Use the bisection method. Assumes nctime is ordered.importbisectindex=numpy.array([bisect.bisect_right(nctime,n)forninnum],int)before=index==0index=numpy.array([bisect.bisect_left(nctime,n)forninnum],int)after=index==Nifselectin['before','exact']andnumpy.any(before):raiseValueError('Some of the times given are before the first time in `nctime`.')ifselectin['after','exact']andnumpy.any(after):raiseValueError('Some of the times given are after the last time in `nctime`.')# Find the times for which the match is not perfect.# Use list comprehension instead of the simpler `nctime[index]` since# not all time objects support numpy integer indexing (eg dap).index[after]=N-1ncnum=numpy.squeeze([nctime[i]foriinindex])mismatch=numpy.nonzero(ncnum!=num)[0]ifselect=='exact':iflen(mismatch)>0:raiseValueError('Some of the times specified were not found in the `nctime` variable.')elifselect=='before':index[after]=Nindex[mismatch]-=1elifselect=='after':passelifselect=='nearest':nearest_to_left=num[mismatch]<numpy.array([float(nctime[i-1])+float(nctime[i])foriinindex[mismatch]])/2.index[mismatch]=index[mismatch]-1*nearest_to_leftelse:raiseValueError("%s is not an option for the `select` argument."%select)# Correct for indices equal to -1index[before]=0# convert numpy scalars or single element arrays to python ints.return_toscalar(index)def_toscalar(a):ifa.shapein[(),(1,)]:returna.item()else:returna