Source code for lsst.obs.base.makeRawVisitInfo

#
# 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
[docs] def defaultMetadata(value, defaultValue, minimum=None, maximum=None): """Return the value if it is not NaN and within min/max, otherwise return defaultValue. Parameters ---------- value : `float` metadata value returned by popItem, popFloat, or popAngle defaultValue : `float`` default value to use if the metadata value is invalid minimum : `float` Minimum possible valid value, optional maximum : `float` Maximum possible valid value, optional Returns ------- `float` The "validated" value. """ if np.isnan(value): retVal = defaultValue else: if minimum is not None and value < minimum: retVal = defaultValue elif maximum is not None and value > maximum: retVal = defaultValue else: retVal = value return retVal