''' Implementation of Julian day.
.. codeauthor:: Yoshifumi Futaana
This :mod:`irfpy.util.julday` package provides implementation of Julian day
and related functionality.
It consists of mainly three classes:
- :class:`Julday` expresses a Julian day.
- :class:`JdObject` expresses a pair of (:class:`Julday` and any object).
- :class:`JdSeries` expresses a series of :class:`JdObject`.
Some helper functions can bridge the :class:`Julday` object and
``datetime.datetime`` object.
.. note::
Currently :class:`JdObject` and :class:`JdSeries` are not used frequently.
Similar functionalities are found at
:class:`irfpy.util.timeseries.ScalarSeries` or
:class:`irfpy.util.timeseries.Vector3dSeries`
These are limited to the floating point values, but much faster and intuitive.
Prioritize to use :mod:`irfpy.util.timeseries` functionality.
.. note::
It is recommended to use `datetime.datetime` for nominal use.
Using :class:`irfpy.util.julday.Julday` is not very recommended;
only for converting to `datetime.datetime`.
See also :mod:`irfpy.util.utc`.
'''
import datetime
import logging
_logger = logging.getLogger(__name__)
import numpy
import bisect
[docs]class Julday:
""" Julian day class.
:param yr_or_jd: Float, :class:`Julday` or ``datetime.datetime``
instance specifying the time.
>>> from irfpy.util.julday import Julday
>>> jd = Julday(2008, 11, 10, 13, 51, 25.371)
>>> print(jd) # doctest: +NORMALIZE_WHITESPACE
2454781.07738(2008-11-10T13:51:25.371)
>>> jd2 = Julday(jd)
>>> print(jd2)
2454781.07738(2008-11-10T13:51:25.371)
>>> dt = jd.getDatetime()
>>> print(dt) # Note: double precision!!
2008-11-10 13:51:25.370980
>>> jd3 = Julday(dt)
>>> print(jd3) # doctest: +NORMALIZE_WHITESPACE
2454781.07738(2008-11-10T13:51:25.371)
Deprecateed method:
Previously , one can get julday in double format by ``jd.jd``
However, this call is not supported anymore. Use :meth:`juld`.
>>> jd = Julday(1, 1, 1, 0, 0, 0) # corresponding to 1721423.5
>>> print('%.5f' % jd.juld())
1721423.50000
>>> jd = Julday(2000000.5) # JD=2000000.5 is 0763-09-15 0:0:0
>>> print(jd.to_string())
0763-09-15T00:00:00.000
.. todo::
It is planned to implement comparison with a natural python way.
"""
HOUR = 1 / 24.0
MINUTE = 1 / 24.0 / 60.0
SECOND = 1 / 24.0 / 60.0 / 60.0
def __init__(self, yr_or_jd, mo=None, dy=None, hr=None, mn=None, se=None):
"""
Instance the Julday object.
"""
if mo is not None:
self.__jd = Julday.cal2jd(yr_or_jd, mo, dy, hr, mn, se)
else:
if isinstance(yr_or_jd, Julday):
self.__jd = yr_or_jd.__jd
elif isinstance(yr_or_jd, datetime.datetime):
self.__jd = Julday.cal2jd(yr_or_jd.year, yr_or_jd.month,
yr_or_jd.day, yr_or_jd.hour, yr_or_jd.minute,
yr_or_jd.second)
self.__jd += (yr_or_jd.microsecond / 86400e6)
else:
self.__jd = yr_or_jd
[docs] def juld(self):
'''Return double precision expression of julian day.
'''
return self.__jd
[docs] @classmethod
def cal2jd(self, year, month, day, hour, minute, second):
''' A class method to calculate the julian day.
'''
if(month < 3):
year = year - 1
month = month + 12
y100 = year // 100
greg = 0
if(year > 1582 or (year == 1582 and month > 10) or
(year == 1582 and month == 10 and day > 4)):
greg = 2 - y100 + (y100 // 4)
else:
greg = 0
ddy = day + hour / 24.0 + minute / 1440. + second / 86400.
jd = int(365.25 * (year + 4716)) + (int)(30.6001 *
(month + 1)) + ddy + greg - 1524.5
return jd
[docs] @classmethod
def jd2cal(self, jd):
''' A class method to return a tuple of calendar.
'''
jd = jd + 0.5
jd_int = (int)(jd)
jd_dbl = jd - jd_int
if jd_int < 2299161:
pass
else:
a = int((jd_int - 1867216.25) / 36524.25)
jd_int = jd_int + 1 + a - (int)(a // 4)
dy = jd_int + 1524
yr = (int)((dy - 122.1) / 365.25)
dy1 = (int)(365.25 * yr)
mo = (int)((dy - dy1) / 30.6001)
day = int(dy - dy1 - (int)(30.6001 * mo))
if mo < 14:
month = mo - 1
else:
month = mo - 13
if month > 2:
year = yr - 4716
else:
year = yr - 4715
hour = (int)(jd_dbl * 24)
jd_dbl = jd_dbl - hour * (1 / 24.)
minute = (int)(jd_dbl * 1440)
jd_dbl = jd_dbl - minute * (1 / 1440.)
second = jd_dbl * 86400
return [year, month, day, hour, minute, second, -1, -1, -1]
[docs] def getDatetime(self):
''' Return the corresponding time in ``datetime.datetime`` format.
'''
yy, mm, dd, HH, MM, SS, t0, t1, t2 = self.jd2cal(self.__jd)
ss = int(SS) # Here SS is double
usec = int((SS - ss) * 1e6)
return datetime.datetime(yy, mm, dd, HH, MM, ss, usec)
[docs] def dayFrom(self, julday):
''' Return the difference in time, i.e. self.juld() - julday.juld().
>>> jd0 = Julday(2009, 1, 10, 12, 0, 0)
>>> jd1 = Julday(2009, 1, 12, 6, 0, 0)
>>> print(jd0.dayFrom(jd1))
-1.75
'''
return self.__jd - julday.__jd
[docs] def dayAfter(self, days):
''' Return the Julday instance after the specified day.
>>> jd0 = Julday(2010, 5, 3, 0, 0, 0)
>>> jd1 = jd0.dayAfter(31.5)
>>> print(jd1.getDatetime())
2010-06-03 12:00:00
'''
return Julday(self.__jd + days)
[docs] def isEarlier(self, julday):
''' Return True if the self is earlier than the specified julday
>>> jd0 = Julday(2005, 1, 1, 0, 0, 0)
>>> jd1 = Julday(2006, 1, 1, 0, 0, 0)
>>> print(jd0.isEarlier(jd1))
True
>>> print(jd0.isEarlier(jd0))
False
>>> print(jd1.isEarlier(jd0))
False
'''
return (self.__jd < julday.__jd)
[docs] def isLater(self, julday):
''' Return True if the self is later than the specified julday
>>> jd0 = Julday(2005, 1, 1, 0, 0, 0)
>>> jd1 = Julday(2006, 1, 1, 0, 0, 0)
>>> print(jd0.isLater(jd1))
False
>>> print(jd0.isLater(jd0))
False
>>> print(jd1.isLater(jd0))
True
'''
return (self.__jd > julday.__jd)
[docs] def isSame(self, julday):
''' Return ``True`` if both the time is the same.
Comparison by values, not object.
>>> jd0 = Julday(2009, 1, 10, 12, 0, 0)
>>> jd1 = Julday(2012, 1, 10, 12, 0, 0)
>>> jd2 = Julday(2009, 1, 10, 12, 0, 0)
>>> jd3 = jd0
>>> print(jd0.isSame(jd0))
True
>>> print(jd0.isSame(jd1))
False
>>> print(jd0.isSame(jd2))
True
>>> print(jd0.isSame(jd3))
True
>>> print(jd0 is jd2)
False
>>> print(jd0 is jd3)
True
'''
return (self.__jd == julday.__jd)
[docs] def to_string(self):
date = self.jd2cal(self.__jd)
return "%04d-%02d-%02dT%02d:%02d:%06.3f" % (
date[0], date[1], date[2], date[3], date[4], date[5])
def __str__(self):
date = self.jd2cal(self.__jd)
return "%14.5f" % self.__jd + "(%04d-%02d-%02dT%02d:%02d:%06.3f)" % (
date[0], date[1], date[2], date[3], date[4], date[5])
def __repr__(self):
date = self.jd2cal(self.__jd)
return 'Julday(%4d, %2d, %2d, %2d, %2d, %06.3f)' % (
date[0], date[1], date[2], date[3], date[4], date[5])
def __hash__(self):
''' Return hash.
>>> jd0 = Julday(2012, 1, 10, 15, 30, 0)
>>> print(jd0)
2455937.14583(2012-01-10T15:30:00.000)
>>> print(hash(jd0))
336268772537366913
'''
return hash(self.__jd)
def __cmp__(self, other):
if other is None:
return -1
return self.__jd - other.__jd
def __eq__(self, other):
if other is None:
return False
# print "COMPARING %f __eq__ %f"%(self.__jd, other.__jd)
return self.__jd == other.__jd
def __lt__(self, other):
if other is None:
return False
# print "COMPARING %f __lt__ %f"%(self.__jd, other.__jd)
return self.__jd < other.__jd
def __le__(self, other):
if other is None:
return False
# print "COMPARING %f __le__ %f"%(self.__jd, other.__jd)
return self.__jd <= other.__jd
def __ne__(self, other):
# print "COMPARING %f __ne__ %f"%(self.__jd, other.__jd)
return not self.__eq__(other)
def __gt__(self, other):
# print "COMPARING %f __gt__ %f"%(self.__jd, other.__jd)
return not self.__le__(other)
def __ge__(self, other):
# print "COMPARING %f __ge__ %f"%(self.__jd, other.__jd)
return not self.__lt__(other)
[docs] def getTimedelta(self, other):
'''Return ``datetime.timedelta`` object for self-other
'''
return self.getDatetime() - other.getDatetime()
[docs]class JdObject:
'''
A pair of julian day and any object.
The object is NOT deep-copied. It is stored by reference (just by '=').
Therefore, if you changed the object after making the JdObject instance,
the contents inside the JdObject would be changed.
>>> jdo = JdObject(Julday(2007, 5, 10, 20, 15, 38), 2734.8)
>>> jd = jdo.getJd()
>>> print(jd)
2454231.34419(2007-05-10T20:15:38.000)
>>> dat = jdo.getData()
>>> print(dat)
2734.8
'''
def __init__(self, julday, object):
if not isinstance(julday, Julday):
raise TypeError(
'Specified julday (%s) is not an instance of Julday'
% julday)
self._jd = julday
self._obj = object
[docs] def getData(self):
''' Return the object.
'''
return self._obj
[docs] def getJd(self):
''' Return the time as :class:`Julday` instance.
'''
return self._jd
[docs] def getDatetime(self):
''' Return the time as ``datetime.datetime`` instance.
'''
return self._jd.getDatetime()
[docs] def jd(self):
''' Return the time as double-precision float.
'''
return self._jd.juld()
def __str__(self):
return self._obj.__str__() + " @ " + self._jd.__str__()
def __hash__(self):
self._jd.__hash__()
def __cmp__(self, other):
if other is None:
return -1
return self._jd.__cmp__(other._jd)
def __eq__(self, other):
if other is None:
return False
return self._jd == other._jd
def __lt__(self, other):
if other is None:
return False
return self._jd < other._jd
def __le__(self, other):
if other is None:
return False
return self._jd <= other._jd
def __ne__(self, other):
return not self.__eq__(other)
def __gt__(self, other):
return not self.__le__(other)
def __ge__(self, other):
return not self.__lt__(other)
[docs]class JdSeriesKeyAlreadyExistException(Exception):
""" Exception thrown when the JdSeries has data on the specified julian day
"""
def __init__(self, value):
self.value = value
def __str__(self):
return repr(self.value)
[docs]class JdSeries:
r"""
Map of julian day and any data (object).
The ``JdSeries`` behaves as a kind of container constituting of
:class:`Julday` and any type of object.
You can instance the container as
>>> ser = JdSeries()
You can add items as
>>> ser.add(Julday(2008, 12, 3, 10, 0, 0), 350.55)
>>> ser.add(Julday(2008, 12, 3, 11, 0, 0), 348.73)
>>> ser.add(Julday(2008, 12, 3, 12, 0, 0), 342.91)
You can know the size of the container.
>>> print(ser.size())
3
If you add same Julian day data,
:class:`JdSeriesKeyAlreadyExistException` is thrown.
>>> try:
... ser.add(Julday(2008, 12, 3, 10, 0, 0), 328.99)
... print('No exeption caught. Failed')
... except JdSeriesKeyAlreadyExistException:
... print('Exception correctly caught. Successful.')
Exception correctly caught. Successful.
You may get data by :meth:`getNeighbor` method with specifying Julian day.
You get :class:`JdObject` instance.
>>> data = ser.getNeighbor(Julday(2008, 12, 3, 10, 0, 0))
>>> print(data.getJd().juld()) # doctest: +ELLIPSIS
2454803.91666...
>>> print(data.getData())
350.55
You can also use :meth:`getNeighbor` method with specifying
non-existing Julian day, but in this case,
you get a data closest to the specified Julday,
and execution time will be 10-100 times slower.
>>> data = ser.getNeighbor(Julday(2008, 12, 3, 10, 10, 0))
>>> print(data)
350.55 @ 2454803.91667(2008-12-03T09:59:60.000)
Especially, just after adding data, :meth:`getNeighbor`
will be slowened because :meth:`getNeighbor` must refresh
an internal cache list. Therefore, data adding is recommended
to be done before getting data if it allows.
To avoid this slower situation, it is recommended to use exact
Julian day you specified as far as it may work.
The registered Julian day list can be obtained by
:class:`getJuldayList` method.
>>> jdlist = ser.getJuldayList()
>>> print(len(jdlist))
3
>>> print('%r\n%r\n%r' % (jdlist[:]))
Julday(2008, 12, 3, 9, 59, 60.000)
Julday(2008, 12, 3, 11, 0, 00.000)
Julday(2008, 12, 3, 12, 0, 00.000)
>>> print(ser.getNeighbor(jdlist[2]).getData())
342.91
To obtain a set of data that includes the specified Julian day,
use :meth:`getBetween` method, which it is also slow.
>>> jd = ser.getNeighbor(Julday(2008, 12, 3, 11, 50, 0))
>>> print(jd)
342.91 @ 2454804.00000(2008-12-03T12:00:00.000)
>>> jd0, jd1 = ser.getBetween(Julday(2008, 12, 3, 11, 50, 0))
>>> print(jd0)
348.73 @ 2454803.95833(2008-12-03T11:00:00.000)
>>> print(jd1)
342.91 @ 2454804.00000(2008-12-03T12:00:00.000)
You can use the iteration. The :meth:`JdObject` instance is iterated.
>>> n = 0
>>> for jdo in ser:
... n = n + 1
>>> print(n)
3
Iteration can of course be repeated.
>>> n = 0
>>> for jdo in ser:
... print(jdo.getData())
350.55
348.73
342.91
"""
def __init__(self):
# _v is a dictionary. Key is Julday and the data is Object
self._v = {}
# Once sorted, _sk is an sorted array of julday in double.
# If you change _v, _sk should be set to None.
self._sk = None
self.iterjd = None # for iteration.
def __getitem__(self, idx):
''' Get the ``JdObject`` for the corresponding index.
'''
jdlist = self.getJdList()
return self.getNeighbor(Julday(jdlist[idx]))
[docs] def add(self, jdobj_or_jd, none_or_obj=None):
"""
Add the data to the series.
If argument size is 1, the jdobj_or_jd is considered as JdObject.
If argument size is 2, the jdobj_or_jd is considered as Julday
and none_or_obj is considered as Object.
"""
# Argument size is 2
jd = jdobj_or_jd
obj = none_or_obj
# Argument size is 1
if none_or_obj is None:
jd = jdobj_or_jd.getJd()
obj = jdobj_or_jd.getData()
if jd in self._v:
raise JdSeriesKeyAlreadyExistException(
"Key already exist. JD=" + jd.__str__())
self._v[jd] = obj
self._sk = None
[docs] def clear(self):
"""Clear the saved data.
>>> x = JdSeries()
>>> x.isEmpty()
True
>>> x.add(Julday(2008, 1, 2, 3, 4, 5), 300)
>>> x.isEmpty()
False
>>> x.clear()
>>> x.isEmpty()
True
"""
self._v = {}
self._sk = None
[docs] def jdsort(self):
''' Make internal cache.
'''
if self._sk is None:
_logger.info('Sorting JdSeries.')
self._sk = []
keys = list(self._v.keys())
for jd in keys:
self._sk.append(jd.juld())
self._sk.sort()
[docs] def size(self):
"""Returns the size of the saved data"""
return len(self._v)
[docs] def isEmpty(self):
""" Return ture if the container size is 0"""
return (len(self._v) == 0)
[docs] def first(self):
""" Return the first :class:`JdObject` """
firstjd = self.firstJd()
firstobj = self._v[firstjd]
return JdObject(firstjd, firstobj)
[docs] def firstJd(self):
""" Return the first :class:`Julday` """
jdlist = self.getJdList()
return Julday(jdlist[0])
[docs] def last(self):
""" Return the last :class:`JdObject` """
lastjd = self.lastJd()
lastobj = self._v[lastjd]
return JdObject(lastjd, lastobj)
[docs] def lastJd(self):
""" Return the last :class:`Julday` """
jdlist = self.getJdList()
return Julday(jdlist[-1])
[docs] def getJdList(self):
""" Return the sorted list of julian day in double precision float """
self.jdsort()
dbllist = [jd for jd in self._sk]
dbllist.sort()
return dbllist
[docs] def getJuldayList(self):
""" Return the sorted list of julian day (:class:`Julday` instances).
This is rather expensive calculations.
Using :meth:`getJdList` that returns the float array is faster.
"""
jdlist = list(self._v.keys())
return tuple(sorted(jdlist))
[docs] def getDatetimeList(self):
''' Return the sorted list of ``datatime.datetime`` instance.
>>> t0 = Julday(2000, 1, 1, 0, 0, 0)
>>> t1 = Julday(2000, 3, 10, 0, 0, 0)
>>> t2 = Julday(1900, 5, 10, 0, 0, 0)
>>> ser = JdSeries()
>>> ser.add(t0, 30)
>>> ser.add(t1, 40)
>>> ser.add(t2, 50)
>>> dtlist = ser.getDatetimeList()
>>> print(dtlist[0])
1900-05-10 00:00:00
'''
dtlist = [jd.getDatetime() for jd in list(self._v.keys())]
return tuple(sorted(dtlist))
def __iter__(self):
"""
Iteration is enabled for JdSeries class.
Iteration is done the order of the julian day.
:class:`JdObject` instance is iterated.
"""
# print >> sys.stderr, " ## Just to know... __iter__ is called."
return self
def __next__(self):
"""
Support iteration.
"""
if self.size() == 0:
self.__stopIteration()
raise StopIteration
if self.iterjd is None:
self.iterjd = 0
else:
self.iterjd = self.iterjd + 1
try:
list = self.getJdList()
jd = list[self.iterjd]
except IndexError as e:
self.__stopIteration()
raise StopIteration
returnobj = self.getNeighbor(Julday(jd))
return returnobj
def __stopIteration(self):
'''Internal function when iteration want to be stopped.
'''
self.iterjd = None
[docs] def hasElement(self, julday):
"""Return that the specified julian day is contained in the dataset """
return (julday in self._v)
[docs] def isIncluded(self, julday):
""" Check the given julday is in the time period.
Return True if the specified julian day is included in the time period
"""
jd0 = self.firstJd()
jd1 = self.lastJd()
return (not jd0.isLater(julday)) and (not jd1.isEarlier(julday))
def __getTimedeltaList(self, julday):
## Obtain keys (JD)
jdkey = list(self._v.keys())
jddiff = []
for jd in jdkey:
dt = jd.getTimedelta(julday)
dts = abs(dt.days * 86400. + dt.seconds)
jddiff.append(dts)
jddiff = numpy.array(jddiff)
return jddiff
[docs] def getNeighbor(self, julday):
"""Get the neighboring data for the specified :class:`Julday`.
Return the :class:`JdObject` that is closest data to the
specified :class:`Julday`.
This method has been extremely slow if the specified julday is
not exactly in the key. However, by using cache (``_sk``),
the performance of multiple calls of this method are now quite quick.
However, using the :class:`Julday` instances returned by
:meth:`getJuldayList` is faster.
"""
if self.hasElement(julday):
# _logger.debug('Exact key specified (%s)'%julday)
return JdObject(julday, self._v[julday])
julday_f = julday.juld()
_logger.debug('### Jd specified=%.5f' % julday_f)
self.jdsort() # sort julday to store _sk if needed.
# _logger.debug('### Key sorted=%s' % jdkey)
idx = bisect.bisect(self._sk, julday_f)
_logger.debug('### Index=%d' % idx)
if idx == 0:
jd = self._sk[0]
elif idx == len(self._sk):
jd = self._sk[-1]
else:
jd0 = self._sk[idx - 1]
jd1 = self._sk[idx]
dt0 = abs(julday_f - jd0)
dt1 = abs(julday_f - jd1)
if dt0 <= dt1:
jd = jd0
else:
jd = jd1
jd = Julday(jd)
_logger.debug('### Jd=%s' % jd)
# jddiff = self.__getTimedeltaList(julday)
# minjd = jddiff.argmin()
# jd = jdkey[minjd]
# _logger.debug('Specified JD=%s/Getting JD=%s' % (julday, jd))
if jd in self._v:
return JdObject(jd, self._v[jd])
else:
raise RuntimeError(
'Fatal. Should not happen in single thread.\n' +
'Not found the key %s.' % str(jd))
[docs] def getBetween(self, julday):
""" Return two JdObject between which the specified julday is located.
Return the 2-element array of JdObject between which the specified
julday is located.
@retval JdObject0 A JdObject of the closest data before julday
@retval JdObject1 A JdObject of the closest data after julday
Note that this method is extremely slow.
"""
if self.size() < 2:
raise RuntimeError(
'This JdSeries only has %d element, but must have >=2.' %
self.size())
jd0 = None
jd1 = None
self.jdsort() # sort julday to store _sk if needed.
if not self.isIncluded(julday):
raise RuntimeError(
'Specified JD=%s is not included in the data range.' % julday)
idx = bisect.bisect(self._sk, julday.juld())
if idx == 0:
jd0 = Julday(self._sk[0])
jd1 = Julday(self._sk[1])
elif idx == len(self._sk):
jd0 = Julday(self._sk[-2])
jd1 = Julday(self._sk[-1])
else:
jd0 = Julday(self._sk[idx - 1])
jd1 = Julday(self._sk[idx])
return (JdObject(jd0, self._v[jd0]), JdObject(jd1, self._v[jd1]))
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')