"""
Helpers for time parsing and conversion.
"""
#c Copyright 2008-2023, the GAVO project <gavo@ari.uni-heidelberg.de>
#c
#c This program is free software, covered by the GNU GPL. See the
#c COPYING file in the source distribution.
import bisect
import datetime
import math
from gavo import utils
from gavo.stc import common
JD_MJD = 2400000.5
[docs]def parseISODT(value):
try:
return utils.parseISODT(value)
except ValueError as ex:
raise common.STCLiteralError(str(ex), value)
[docs]@utils.document
def jdnToDateTime(jd):
"""returns a ``datetime.datetime`` instance for a julian day number.
"""
return jYearToDateTime((jd-2451545.0)/365.25+2000.0)
[docs]@utils.document
def mjdToDateTime(mjd):
"""returns a ``datetime.datetime`` instance for a modified julian day number.
Beware: This loses a couple of significant digits due to transformation
to jd.
"""
return jdnToDateTime(mjd+JD_MJD)
[docs]@utils.document
def bYearToDateTime(bYear):
"""returns a datetime.datetime instance for a fractional Besselian year.
This uses the formula given by Lieske, J.H., A&A 73, 282 (1979).
"""
jdn = (bYear-1900.0)*common.tropicalYear+2415020.31352
return jdnToDateTime(jdn)
[docs]@utils.document
def jYearToDateTime(jYear):
"""returns a datetime.datetime instance for a fractional (julian) year.
This refers to time specifications like J2001.32.
"""
return datetime.datetime(2000, 1, 1, 12)+datetime.timedelta(
days=(jYear-2000.0)*365.25)
dtJ2000 = jYearToDateTime(2000.0)
dtB1950 = bYearToDateTime(1950.0)
[docs]@utils.document
def dateTimeToJdn(dt):
"""returns a julian day number (including fractionals) from a datetime
instance.
"""
a = (14-dt.month)//12
y = dt.year+4800-a
m = dt.month+12*a-3
jdn = dt.day+(153*m+2)//5+365*y+y//4-y//100+y//400-32045
try:
secsOnDay = dt.hour*3600+dt.minute*60+dt.second+dt.microsecond/1e6
except AttributeError:
secsOnDay = 0
return jdn+(secsOnDay-43200)/86400.
[docs]@utils.document
def dateTimeToMJD(dt):
"""returns a modified julian date for a datetime instance.
"""
return dateTimeToJdn(dt)-JD_MJD
[docs]def dateTimeToBYear(dt):
return (dateTimeToJdn(dt)-2415020.31352)/common.tropicalYear+1900.0
[docs]@utils.document
def dateTimeToJYear(dt):
"""returns a fractional (julian) year for a datetime.datetime instance.
"""
return (dateTimeToJdn(dt)-2451545)/365.25+2000
[docs]def getSeconds(td):
"""returns the number of seconds corresponding to a timedelta object.
"""
return td.days*86400+td.seconds+td.microseconds*1e-6
############ Time scale conversions
# All these convert to/from TT.
_TDTminusTAI = datetime.timedelta(seconds=32.184)
[docs]@utils.document
def TTtoTAI(tdt):
"""returns TAI for a (datetime.datetime) TDT.
"""
return tdt+_TDTminusTAI
[docs]@utils.document
def TAItoTT(tai):
"""returns TDT for a (datetime.datetime) TAI.
"""
return tai-_TDTminusTAI
def _getTDBOffset(tdb):
"""returns the TDB-TDT according to [EXS] 2.222-1.
"""
g = (357.53+0.9856003*(dateTimeToJdn(tdb)-2451545.0))/180*math.pi
return datetime.timedelta(0.001658*math.sin(g)+0.000014*math.sin(2*g))
[docs]def TDBtoTT(tdb):
"""returns an approximate TT from a TDB.
The simplified formula 2.222-1 from [EXS] is used.
"""
return tdb-_getTDBOffset(tdb)
[docs]def TTtoTDB(tt):
"""returns approximate TDB from TT.
The simplified formula 2.222-1 from [EXS] is used.
"""
return tt+_getTDBOffset(tt)
_L_G = 6.969291e-10 # [EXS], p. 47
def _getTCGminusTT(dt): # [EXS], 2.223-5
return datetime.timedelta(seconds=
_L_G*(dateTimeToJdn(dt)-2443144.5)*86400)
[docs]def TTtoTCG(tt):
"""returns TT from TCG.
This uses 2.223-5 from [EXS].
"""
return tt+_getTCGminusTT(tt)
[docs]def TCGtoTT(tcg):
"""returns TT from TCG.
This uses 2.223-5 from [EXS].
"""
return tcg+_getTCGminusTT(tcg)
_L_B = 1.550505e-8
def _getTCBminusTDB(dt): # [EXS], 2.223-2
return datetime.timedelta(
seconds=_L_B*(dateTimeToJdn(dt)-2443144.5)*86400)
[docs]def TCBtoTT(tcb):
"""returns an approximate TCB from a TT.
This uses [EXS] 2.223-2 and the approximate conversion from TDB to TT.
"""
return TDBtoTT(tcb+_getTCBminusTDB(tcb))
[docs]def TTtoTCB(tt):
"""returns an approximate TT from a TCB.
This uses [EXS] 2.223-2 and the approximate conversion from TT to TDB.
"""
return TTtoTDB(tt)-_getTCBminusTDB(tt)
def _makeLeapSecondTable():
lsTable = []
for lsCount, lsMoment in enumerate([ # from Lenny tzinfo
datetime.datetime(1971, 12, 31, 23, 59, 59),
datetime.datetime(1972, 0o6, 30, 23, 59, 59),
datetime.datetime(1972, 12, 31, 23, 59, 59),
datetime.datetime(1973, 12, 31, 23, 59, 59),
datetime.datetime(1974, 12, 31, 23, 59, 59),
datetime.datetime(1975, 12, 31, 23, 59, 59),
datetime.datetime(1976, 12, 31, 23, 59, 59),
datetime.datetime(1977, 12, 31, 23, 59, 59),
datetime.datetime(1978, 12, 31, 23, 59, 59),
datetime.datetime(1979, 12, 31, 23, 59, 59),
datetime.datetime(1981, 0o6, 30, 23, 59, 59),
datetime.datetime(1982, 0o6, 30, 23, 59, 59),
datetime.datetime(1983, 0o6, 30, 23, 59, 59),
datetime.datetime(1985, 0o6, 30, 23, 59, 59),
datetime.datetime(1987, 12, 31, 23, 59, 59),
datetime.datetime(1989, 12, 31, 23, 59, 59),
datetime.datetime(1990, 12, 31, 23, 59, 59),
datetime.datetime(1992, 0o6, 30, 23, 59, 59),
datetime.datetime(1993, 0o6, 30, 23, 59, 59),
datetime.datetime(1994, 0o6, 30, 23, 59, 59),
datetime.datetime(1995, 12, 31, 23, 59, 59),
datetime.datetime(1997, 0o6, 30, 23, 59, 59),
datetime.datetime(1998, 12, 31, 23, 59, 59),
datetime.datetime(2005, 12, 31, 23, 59, 59),
datetime.datetime(2008, 12, 31, 23, 59, 59),
]):
lsTable.append((lsMoment, datetime.timedelta(seconds=lsCount+10)))
return lsTable
# A table of TAI-UTC, indexed by UTC
leapSecondTable = _makeLeapSecondTable()
del _makeLeapSecondTable
_sentinelTD = datetime.timedelta(seconds=0)
[docs]def getLeapSeconds(dt, table=leapSecondTable):
"""returns TAI-UTC for the datetime dt.
"""
ind = bisect.bisect_left(leapSecondTable, (dt, _sentinelTD))
if ind==0:
return datetime.timedelta(seconds=9.)
return table[ind-1][1]
[docs]def UTCtoTT(utc):
"""returns TT from UTC.
The leap second table is complete through 2009-5.
>>> getLeapSeconds(datetime.datetime(1998,12,31,23,59,58))
datetime.timedelta(seconds=31)
>>> TTtoTAI(UTCtoTT(datetime.datetime(1998,12,31,23,59,59)))
datetime.datetime(1999, 1, 1, 0, 0, 30)
>>> TTtoTAI(UTCtoTT(datetime.datetime(1999,1,1,0,0,0)))
datetime.datetime(1999, 1, 1, 0, 0, 32)
"""
return TAItoTT(utc+getLeapSeconds(utc))
# A table of TAI-UTC, indexed by TT
ttLeapSecondTable = [(UTCtoTT(t), dt)
for t, dt in leapSecondTable]
[docs]def TTtoUTC(tt):
"""returns UTC from TT.
The leap second table is complete through 2009-5.
>>> TTtoUTC(UTCtoTT(datetime.datetime(1998,12,31,23,59,59)))
datetime.datetime(1998, 12, 31, 23, 59, 59)
>>> TTtoUTC(UTCtoTT(datetime.datetime(1999,1,1,0,0,0)))
datetime.datetime(1999, 1, 1, 0, 0)
"""
# XXX TODO: leap seconds need to be computed from UTC, so this will
# be one second off in the immediate vicinity of a leap second.
return TTtoTAI(tt)-getLeapSeconds(tt, ttLeapSecondTable)
# A dict mapping timescales to conversions to/from TT.
timeConversions = {
"UTC": (UTCtoTT, TTtoUTC),
"TCB": (TCBtoTT, TTtoTCB),
"TCG": (TCGtoTT, TTtoTCG),
"TDB": (TDBtoTT, TTtoTDB),
"TAI": (TAItoTT, TTtoTAI),
"TT": (utils.identity, utils.identity),
}
[docs]@ utils.document
def isMJD(col):
"""returns True if the rscdef.Column instance col likely contains MJD values.
This has a long and winding history in DaCHS, and so this is a disaster
of heuristics.
"""
if col.unit!="d":
return False
# classic DaCHS: have a name with mjd in it
if "mjd" in col.name:
return True
# SIAP: a stupid uype
if col.ucd and "MJD" in col.ucd.upper():
return True
# xtype abuse: just here for backward compatibility; this is pain because
# something else will have to kill the xtype.
if col.xtype=="mjd":
return True
# the current recommended way to declare MJD: use an appropriate time0
# in a votable:Coords annotation
for role in col.dmRoles:
# resolve weakref
ann = role()
if ann.name!="location":
# it's not the annotation of a time part of votable:Coords
continue
try:
time0 = ann.instance()["time"]["frame"]["time0"]
if time0=="MJD-origin":
return True
elif time0=="JD-origin":
return False
else:
return float(time0)==JD_MJD
except KeyError:
# no proper frame annotation in this role; keep on trying
pass
# I could inspect column statistics (JD 50000 would be extremely unlikely),
# but there's an overdose of heuristics here as is.
return False
[docs]def datetimeMapperFactory(colDesc):
import time
# TODO: take a long, hard look at this. I *think* essentially all of
# this should be replaced by displayHint-s, if only because neither
# timestamp nor date should have a unit. However, this stuff is executed.
# We need to work out where...
if (colDesc["dbtype"]=="timestamp"
or colDesc["dbtype"]=="date"
# must look in original, as the one in colDesc comes from typesystems
or colDesc.original.xtype=="adql:TIMESTAMP" # legacy, delete ~2020
or colDesc.original.xtype=="timestamp"):
unit = colDesc["unit"]
if (colDesc["displayHint"].get("format")=="humanDate"
or colDesc.original.xtype=="adql:TIMESTAMP" # legacy, delete ~2020
or colDesc.original.xtype=="timestamp"):
fun = lambda val: (val and val.isoformat()) or None
destType = ("char", "*")
colDesc["nullvalue"] = ""
elif (colDesc["ucd"] and "MJD" in colDesc["ucd"].upper()
or colDesc["xtype"]=="mjd"
or "mjd" in colDesc["name"]):
colDesc["unit"] = "d"
fun = lambda val: (val and dateTimeToMJD(val))
destType = ("double", '1')
colDesc["nullvalue"] = "NaN"
colDesc["xtype"] = None
elif unit=="yr" or unit=="a":
fun = lambda val: (val and dateTimeToJYear(val))
def fun(val):
return (val and dateTimeToJYear(val))
return str(val)
destType = ("double", '1')
colDesc["nullvalue"] = "NaN"
colDesc["xtype"] = None
elif unit=="d":
fun = lambda val: (val and dateTimeToJdn(val))
destType = ("double", '1')
colDesc["nullvalue"] = "NaN"
colDesc["xtype"] = None
elif unit=="s":
fun = lambda val: (val and time.mktime(val.timetuple()))
destType = ("double", '1')
colDesc["nullvalue"] = "NaN"
colDesc["xtype"] = None
else:
# default for datetime column: serialise to timestamp
fun = lambda val: (val and val.isoformat()) or None
destType = ("char", "*")
colDesc["nullvalue"] = ""
colDesc["xtype"] = "timestamp"
colDesc["datatype"], colDesc["arraysize"] = destType
return fun
utils.registerDefaultMF(datetimeMapperFactory)
def _test():
import doctest
doctest.testmod()
if __name__=="__main__":
_test()