Skip to content

Commit

Permalink
Fix epochs and add observatory information
Browse files Browse the repository at this point in the history
  • Loading branch information
@cmsaunders
cmsaunders committed Oct 5, 2022
1 parent dcba639 commit 02ce01f
Show file tree
Hide file tree
Showing 2 changed files with 33 additions and 14 deletions.
42 changes: 30 additions & 12 deletions python/lsst/drp/tasks/wcsfit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
import numpy as np
import astropy.time
import astropy.units as u
import astropy.coordinates
import yaml
import wcsfit
import pandas as pd
Expand Down Expand Up @@ -278,12 +279,12 @@ class WCSFitConfig(pipeBase.PipelineTaskConfig,
)
systematicError = pexConfig.Field(
dtype=float,
doc="Systematic error padding for the science images (arcsec)",
doc="Systematic error padding for the science images (marcsec)",
default=0.0
)
referenceSystematicError = pexConfig.Field(
dtype=float,
doc="Systematic error padding for the reference catalog (arcsec)",
doc="Systematic error padding for the reference catalog (marcsec)",
default=0.0
)
modelComponents = pexConfig.ListField(
Expand Down Expand Up @@ -378,6 +379,9 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):

instrumentName = butlerQC.quantum.dataId['instrument']

sampleRefCat = inputs['referenceCatalog'][0].get()
refEpoch = sampleRefCat[0]['epoch']

self.refObjectLoader = ReferenceObjectLoader(
dataIds=[ref.datasetRef.dataId
for ref in inputRefs.referenceCatalog],
Expand All @@ -386,14 +390,14 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
self.refObjectLoader.config.requireProperMotion = self.config.requireProperMotion
self.refObjectLoader.config.anyFilterMapsToThis = self.config.anyFilterMapsToThis

output = self.run(**inputs, instrumentName=instrumentName)
output = self.run(**inputs, instrumentName=instrumentName, refEpoch=refEpoch)

for outputRef in outputRefs.outputWcs:
visit = outputRef.dataId['visit']
butlerQC.put(output.outputWCSs[visit], outputRef)
butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)

def run(self, inputCatalogRefs, inputVisitSummary, instrumentName=None):
def run(self, inputCatalogRefs, inputVisitSummary, instrumentName=None, refEpoch=None):
"""Run the WCS fit for a given set of visits
Parameters
Expand Down Expand Up @@ -422,7 +426,8 @@ def run(self, inputCatalogRefs, inputVisitSummary, instrumentName=None):
instrument = wcsfit.Instrument(instrumentName)

# Get RA, Dec, MJD, etc., for the input visits
exposuresHelper, medianEpoch = self._get_exposure_info(inputVisitSummary, instrument)
exposuresHelper, medianEpoch = self._get_exposure_info(inputVisitSummary, instrument,
refEpoch=refEpoch)

# Get information about the extent of the input visits
fields, fieldCenter, fieldRadius = self._prep_sky(inputVisitSummary, medianEpoch)
Expand All @@ -435,7 +440,7 @@ def run(self, inputCatalogRefs, inputVisitSummary, instrumentName=None):
(self.config.matchRadius * u.arcsec).to(u.degree).value)

# Add the reference catalog to the associator
self._load_refCat(associations, fieldCenter, fieldRadius, epoch=medianEpoch)
self._load_refCat(associations, fieldCenter, fieldRadius, epoch=refEpoch)

# Add the science catalogs and associate new sources as they are added
self._load_catalogs_and_associate(associations, inputCatalogRefs)
Expand Down Expand Up @@ -508,7 +513,8 @@ def _prep_sky(self, inputVisitSummaries, epoch, fieldName='Field'):

return fields, center, radius

def _get_exposure_info(self, visitSummaryTables, instrument, fieldNumber=0, instrumentNumber=0):
def _get_exposure_info(self, visitSummaryTables, instrument, fieldNumber=0, instrumentNumber=0,
refEpoch=None):
"""Get various information about the input visits to feed to the
fitting routines.
Expand Down Expand Up @@ -538,6 +544,7 @@ def _get_exposure_info(self, visitSummaryTables, instrument, fieldNumber=0, inst
airmasses = []
exposureTimes = []
mjds = []
observatories = []
wcss = []

extensionType = []
Expand All @@ -558,7 +565,17 @@ def _get_exposure_info(self, visitSummaryTables, instrument, fieldNumber=0, inst
airmasses.append(visitInfo.getBoresightAirmass())
exposureTimes.append(visitInfo.getExposureTime())
obsDate = visitInfo.getDate()
mjds.append(obsDate.get(obsDate.MJD))
obsMJD = obsDate.get(obsDate.MJD)
mjds.append(obsMJD)
# Get the observatory ICRS position for use in fitting parallax
obsLon = visitInfo.observatory.getLongitude().asDegrees()
obsLat = visitInfo.observatory.getLatitude().asDegrees()
obsElev = visitInfo.observatory.getElevation()
earthLocation = astropy.coordinates.EarthLocation.from_geodetic(obsLon, obsLat, obsElev)
observatory_gcrs = earthLocation.get_gcrs(astropy.time.Time(obsMJD, format='mjd'))
observatory_icrs = observatory_gcrs.transform_to(astropy.coordinates.ICRS())
# We want the position in AU in Cartesian coordinates
observatories.append(observatory_icrs.cartesian.xyz.to(u.AU).value)

for row in visitSummary:
detector = row['id']
Expand All @@ -583,7 +600,7 @@ def _get_exposure_info(self, visitSummaryTables, instrument, fieldNumber=0, inst

# Set the reference epoch to be the median of the science visits.
# The reference catalog will be shifted to this date.
medianEpoch = np.median(mjds)
medianEpoch = astropy.time.Time(np.median(mjds), format='mjd').decimalyear

# Add information for the reference catalog. Most of the values are
# not used.
Expand All @@ -598,7 +615,8 @@ def _get_exposure_info(self, visitSummaryTables, instrument, fieldNumber=0, inst
self.decs.append(0.0)
airmasses.append(0.0)
exposureTimes.append(0)
mjds.append(medianEpoch)
mjds.append((refEpoch if (refEpoch is not None) else medianEpoch))
observatories.append(np.array([0, 0, 0]))
identity = wcsfit.IdentityMap()
icrs = wcsfit.SphericalICRS()
refWcs = wcsfit.Wcs(identity, icrs, "Identity", np.pi / 180.)
Expand All @@ -624,7 +642,8 @@ def _get_exposure_info(self, visitSummaryTables, instrument, fieldNumber=0, inst
self.decs,
airmasses,
exposureTimes,
mjds)
mjds,
observatories)

return exposuresHelper, medianEpoch

Expand Down Expand Up @@ -684,7 +703,6 @@ def _load_refCat(self, associations, center, radius, epoch=None, fieldIndex=0, i
detectorIndex = self.extensionInfo.iloc[extensionIndex]['detectorIndex']
refWcs = self.extensionInfo.iloc[extensionIndex]['wcs']

# TODO: Any way to include the proper motion information here?
associations.addCatalog(refWcs, "STELLAR", visitIndex, fieldIndex, instrumentIndex, detectorIndex,
extensionIndex, np.ones(len(refCat), dtype=bool),
ra, dec, np.arange(len(ra)))
Expand Down
5 changes: 3 additions & 2 deletions tests/test_wcsfit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,7 @@ def setUpClass(cls):
cls.fieldNumber = 0
cls.instrumentName = 'HSC'
cls.instrument = wcsfit.Instrument(cls.instrumentName)
cls.refEpoch = 57205.5

# Make test inputVisitSummary. VisitSummaryTables are taken from
# collection HSC/runs/RC2/w_2022_20/DM-34794
Expand Down Expand Up @@ -473,7 +474,7 @@ def test_make_outputs(self):
WCSFitTask.refObjectLoader = self.refObjectLoader

outputs = WCSFitTask.run(self.inputCatalogRefs, self.inputVisitSummary,
instrumentName=self.instrumentName)
instrumentName=self.instrumentName, refEpoch=self.refEpoch)

for v, visit in enumerate(self.testVisits):
visitSummary = self.inputVisitSummary[v]
Expand All @@ -499,7 +500,7 @@ def test_run(self):
WCSFitTask.refObjectLoader = self.refObjectLoader

outputs = WCSFitTask.run(self.inputCatalogRefs, self.inputVisitSummary,
instrumentName=self.instrumentName)
instrumentName=self.instrumentName, refEpoch=self.refEpoch)

outputMaps = outputs.fitModel.mapCollection.getParamDict()

Expand Down

0 comments on commit 02ce01f

Please sign in to comment.