#
# LSST Data Management System
# Copyright 2016 LSST Corporation.
#
# This product includes software developed by the
# LSST Project (http://www.lsst.org/).
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the LSST License Statement and
# the GNU General Public License along with this program. If not,
# see <http://www.lsstcorp.org/LegalNotices/>.
#
import math
import numpy as np
import astropy.coordinates
import astropy.time
import astropy.units
from lsst.log import Log
from lsst.daf.base import DateTime
from lsst.afw.geom import degrees
from lsst.afw.image import VisitInfo
__all__ = ["MakeRawVisitInfo"]
PascalPerMillibar = 100.0
PascalPerMmHg = 133.322387415 # from Wikipedia; exact
PascalPerTorr = 101325.0/760.0 # from Wikipedia; exact
KelvinMinusCentigrade = 273.15 # from Wikipedia; exact
# have these read at need, to avoid unexpected errors later
NaN = float("nan")
BadDate = DateTime()
[docs]class MakeRawVisitInfo(object):
"""Base class functor to make a VisitInfo from the FITS header of a raw image.
A subclass will be wanted for each camera. Subclasses should override:
- `setArgDict`, The override can call the base implementation,
which simply sets exposure time and date of observation
- `getDateAvg`
The design philosophy is to make a best effort and log warnings of problems,
rather than raising exceptions, in order to extract as much VisitInfo information as possible
from a messy FITS header without the user needing to add a lot of error handling.
However, the methods that transform units are less forgiving; they assume
the user provides proper data types, since type errors in arguments to those
are almost certainly due to coding mistakes.
Parameters
----------
log : `lsst.log.Log` or None
Logger to use for messages.
(None to use ``Log.getLogger("MakeRawVisitInfo")``).
"""
def __init__(self, log=None):
if log is None:
log = Log.getLogger("MakeRawVisitInfo")
self.log = log
[docs] def __call__(self, md, exposureId):
"""Construct a VisitInfo and strip associated data from the metadata.
Parameters
----------
md : `lsst.daf.base.PropertyList` or `lsst.daf.base.PropertySet`
Metadata to pull from.
Items that are used are stripped from the metadata (except TIMESYS,
because it may apply to other keywords).
exposureId : `int`
exposure ID
Notes
-----
The basic implementation sets `date` and `exposureTime` using typical values
found in FITS files and logs a warning if neither can be set.
"""
argDict = dict(exposureId=exposureId)
self.setArgDict(md, argDict)
for key in list(argDict.keys()): # use a copy because we may delete items
if argDict[key] is None:
self.log.warn("argDict[{}] is None; stripping".format(key, argDict[key]))
del argDict[key]
return VisitInfo(**argDict)
[docs] def setArgDict(self, md, argDict):
"""Fill an argument dict with arguments for VisitInfo and pop associated metadata
Subclasses are expected to override this method, though the override
may wish to call this default implementation, which:
- sets exposureTime from "EXPTIME"
- sets date by calling getDateAvg
Parameters
----------
md : `lsst.daf.base.PropertyList` or `PropertySet`
Metadata to pull from.
Items that are used are stripped from the metadata (except TIMESYS,
because it may apply to other keywords).
argdict : `dict`
dict of arguments
Notes
-----
Subclasses should expand this or replace it.
"""
argDict["exposureTime"] = self.popFloat(md, "EXPTIME")
argDict["date"] = self.getDateAvg(md=md, exposureTime=argDict["exposureTime"])
[docs] def getDateAvg(self, md, exposureTime):
"""Return date at the middle of the exposure.
Parameters
----------
md : `lsst.daf.base.PropertyList` or `PropertySet`
Metadata to pull from.
Items that are used are stripped from the metadata (except TIMESYS,
because it may apply to other keywords).
exposureTime : `float`
Exposure time (sec)
Notes
-----
Subclasses must override. Here is a typical implementation::
dateObs = self.popIsoDate(md, "DATE-OBS")
return self.offsetDate(dateObs, 0.5*exposureTime)
"""
raise NotImplementedError()
[docs] def getDarkTime(self, argDict):
"""Get the darkTime from the DARKTIME keyword, else expTime, else NaN,
If dark time is available then subclasses should call this method by
putting the following in their `__init__` method::
argDict['darkTime'] = self.getDarkTime(argDict)
Parameters
----------
argdict : `dict`
Dict of arguments.
Returns
-------
`float`
Dark time, as inferred from the metadata.
"""
darkTime = argDict.get("darkTime", NaN)
if np.isfinite(darkTime):
return darkTime
self.log.info("darkTime is NaN/Inf; using exposureTime")
exposureTime = argDict.get("exposureTime", NaN)
if not np.isfinite(exposureTime):
raise RuntimeError("Tried to substitute exposureTime for darkTime but it is not available")
return exposureTime
[docs] def offsetDate(self, date, offsetSec):
"""Return a date offset by a specified number of seconds.
date : `lsst.daf.base.DateTime`
Date baseline to offset from.
offsetSec : `float`
Offset, in seconds.
Returns
-------
`lsst.daf.base.DateTime`
The offset date.
"""
if not date.isValid():
self.log.warn("date is invalid; cannot offset it")
return date
if math.isnan(offsetSec):
self.log.warn("offsetSec is invalid; cannot offset date")
return date
dateNSec = date.nsecs(DateTime.TAI)
return DateTime(dateNSec + int(offsetSec*1.0e9), DateTime.TAI)
[docs] def popItem(self, md, key, default=None):
"""Remove an item of metadata and return the value.
Log a warning if the key is not found.
Parameters
----------
md : `lsst.daf.base.PropertyList` or `PropertySet`
Metadata to pull `key` from and remove.
key : `str`
Metadata key to extract.
default : `object`
Value to return if key not found.
Returns
-------
`object`
The value of the specified key, using whatever type md.getScalar(key)
returns.
"""
try:
if not md.exists(key):
self.log.warn("Key=\"{}\" not in metadata".format(key))
return default
val = md.getScalar(key)
md.remove(key)
return val
except Exception as e:
# this should never happen, but is a last ditch attempt to avoid exceptions
self.log.warn("Could not read key=\"{}\" in metadata: {}" % (key, e))
return default
[docs] def popFloat(self, md, key):
"""Pop a float with a default of NaN.
Parameters
----------
md : `lsst.daf.base.PropertyList` or `PropertySet`
Metadata to pull `key` from and remove.
key : `str`
Key to read and remove from md.
Returns
-------
`float`
Value of the requested key as a float; float("nan") if the key is
not found.
"""
val = self.popItem(md, key, default=NaN)
try:
return float(val)
except Exception as e:
self.log.warn("Could not interpret {} value {} as a float: {}".format(key, repr(val), e))
return NaN
[docs] def popAngle(self, md, key, units=astropy.units.deg):
"""Pop an lsst.afw.geom.Angle, whose metadata is in the specified units, with a default of Nan
The angle may be specified as a float or sexagesimal string with 1-3 fields.
Parameters
----------
md : `lsst.daf.base.PropertyList` or `PropertySet`
Metadata to pull `key` from and remove.
key : `str`
Key to read and remove from md.
Returns
-------
`lsst.afw.geom.Angle`
Value of the requested key as an angle; Angle(NaN) if the key is
not found.
"""
angleStr = self.popItem(md, key, default=None)
if angleStr is not None:
try:
return (astropy.coordinates.Angle(angleStr, unit=units).deg)*degrees
except Exception as e:
self.log.warn("Could not intepret {} value {} as an angle: {}".format(key, repr(angleStr), e))
return NaN*degrees
[docs] def popIsoDate(self, md, key, timesys=None):
"""Pop a FITS ISO date as an lsst.daf.base.DateTime
Parameters
----------
md : `lsst.daf.base.PropertyList` or `PropertySet`
Metadata to pull `key` from and remove.
key : `str`
Date key to read and remove from md.
timesys : `str`
Time system as a string (not case sensitive), e.g. "UTC" or None;
if None then look for TIMESYS (but do NOT pop it, since it may be
used for more than one date) and if not found, use UTC.
Returns
-------
`lsst.daf.base.DateTime`
Value of the requested date; `DateTime()` if the key is not found.
"""
isoDateStr = self.popItem(md=md, key=key)
if isoDateStr is not None:
try:
if timesys is None:
timesys = md.getScalar("TIMESYS") if md.exists("TIMESYS") else "UTC"
if isoDateStr.endswith("Z"): # illegal in FITS
isoDateStr = isoDateStr[0:-1]
astropyTime = astropy.time.Time(isoDateStr, scale=timesys.lower(), format="fits")
# DateTime uses nanosecond resolution, regardless of the resolution of the original date
astropyTime.precision = 9
# isot is ISO8601 format with "T" separating date and time and no time zone
return DateTime(astropyTime.tai.isot, DateTime.TAI)
except Exception as e:
self.log.warn("Could not parse {} = {} as an ISO date: {}".format(key, isoDateStr, e))
return BadDate
[docs] def popMjdDate(self, md, key, timesys=None):
"""Get a FITS MJD date as an ``lsst.daf.base.DateTime``.
Parameters
----------
md : `lsst.daf.base.PropertyList` or `PropertySet`
Metadata to pull `key` from and remove.
key : `str`
Date key to read and remove from md.
timesys : `str`
Time system as a string (not case sensitive), e.g. "UTC" or None;
if None then look for TIMESYS (but do NOT pop it, since it may be
used for more than one date) and if not found, use UTC.
Returns
-------
`lsst.daf.base.DateTime`
Value of the requested date; `DateTime()` if the key is not found.
"""
mjdDate = self.popFloat(md, key)
try:
if timesys is None:
timesys = md.getScalar("TIMESYS") if md.exists("TIMESYS") else "UTC"
astropyTime = astropy.time.Time(mjdDate, format="mjd", scale=timesys.lower())
# DateTime uses nanosecond resolution, regardless of the resolution of the original date
astropyTime.precision = 9
# isot is ISO8601 format with "T" separating date and time and no time zone
return DateTime(astropyTime.tai.isot, DateTime.TAI)
except Exception as e:
self.log.warn("Could not parse {} = {} as an MJD date: {}".format(key, mjdDate, e))
return BadDate
@staticmethod
[docs] def eraFromLstAndLongitude(lst, longitude):
"""
Return an approximate Earth Rotation Angle (afw:Angle) computed from
local sidereal time and longitude (both as afw:Angle; Longitude shares
the afw:Observatory covention: positive values are E of Greenwich).
NOTE: if we properly compute ERA via UT1 a la DM-8053, we should remove
this method.
"""
return lst - longitude
@staticmethod
[docs] def altitudeFromZenithDistance(zd):
"""Convert zenith distance to altitude (lsst.afw.geom.Angle)"""
return 90*degrees - zd
@staticmethod
[docs] def centigradeFromKelvin(tempK):
"""Convert temperature from Kelvin to Centigrade"""
return tempK - KelvinMinusCentigrade
@staticmethod
[docs] def pascalFromMBar(mbar):
"""Convert pressure from millibars to Pascals
"""
return mbar*PascalPerMillibar
@staticmethod
[docs] def pascalFromMmHg(mmHg):
"""Convert pressure from mm Hg to Pascals
Notes
-----
Could use the following, but astropy.units.cds is not fully compatible with Python 2
as of astropy 1.2.1 (see https://github.com/astropy/astropy/issues/5350#issuecomment-248612824):
astropy.units.cds.mmHg.to(astropy.units.pascal, mmHg)
"""
return mmHg*PascalPerMmHg
@staticmethod
[docs] def pascalFromTorr(torr):
"""Convert pressure from torr to Pascals
"""
return torr*PascalPerTorr
@staticmethod