Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-37943: Add starCatalog output #13

Merged
merged 2 commits into from
May 1, 2023
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
26 changes: 18 additions & 8 deletions python/lsst/drp/tasks/gbdesAstrometricFit.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@
import yaml
import wcsfit
import astshim
import pyarrow as pa

import lsst.geom
import lsst.pex.config as pexConfig
Expand Down Expand Up @@ -261,9 +260,15 @@ class GbdesAstrometricFitConnections(pipeBase.PipelineTaskConnections,
doc=("Source table with stars used in fit, along with residuals in pixel coordinates and tangent "
"plane coordinates and chisq values."),
name='gbdesAstrometricFit_fitStars',
storageClass='ArrowTable',
storageClass='ArrowNumpyDict',
dimensions=('instrument', 'skymap', 'tract', 'physical_filter'),
)
starCatalog = pipeBase.connectionTypes.Output(
doc="",
name='gbdesAstrometricFit_starCatalog',
storageClass='ArrowNumpyDict',
dimensions=('instrument', 'skymap', 'tract', 'physical_filter')
)


class GbdesAstrometricFitConfig(pipeBase.PipelineTaskConfig,
Expand Down Expand Up @@ -441,6 +446,7 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
visit = outputRef.dataId['visit']
butlerQC.put(output.outputWCSs[visit], outputRef)
butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
butlerQC.put(output.starCatalog, outputRefs.starCatalog)

def run(self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch=None,
refObjectLoader=None):
Expand Down Expand Up @@ -478,7 +484,7 @@ def run(self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch

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

# Get information about the extent of the input visits
fields, fieldCenter, fieldRadius = self._prep_sky(inputVisitSummaries, exposureInfo.medianEpoch)
Expand All @@ -491,8 +497,9 @@ def run(self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch
(self.config.matchRadius * u.arcsec).to(u.degree).value)

# Add the reference catalog to the associator
medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format='decimalyear').mjd
refObjects, refCovariance = self._load_refcat(associations, refObjectLoader, fieldCenter, fieldRadius,
extensionInfo, epoch=refEpoch)
extensionInfo, epoch=medianEpoch)

# Add the science catalogs and associate new sources as they are added
sourceIndices, usedColumns = self._load_catalogs_and_associate(associations, inputCatalogRefs,
Expand Down Expand Up @@ -533,11 +540,13 @@ def run(self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch
self.log.info("WCS fitting done")

outputWCSs = self._make_outputs(wcsf, inputVisitSummaries, exposureInfo)
outputCatalog = pa.Table.from_pydict(wcsf.getOutputCatalog())
outputCatalog = wcsf.getOutputCatalog()
starCatalog = wcsf.getStarCatalog()

return pipeBase.Struct(outputWCSs=outputWCSs,
fitModel=wcsf,
outputCatalog=outputCatalog)
outputCatalog=outputCatalog,
starCatalog=starCatalog)

def _prep_sky(self, inputVisitSummaries, epoch, fieldName='Field'):
"""Get center and radius of the input tract. This assumes that all
Expand Down Expand Up @@ -692,7 +701,8 @@ def _get_exposure_info(self, inputVisitSummaries, instrument, fieldNumber=0, ins

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

# Add information for the reference catalog. Most of the values are
# not used.
Expand All @@ -707,7 +717,7 @@ def _get_exposure_info(self, inputVisitSummaries, instrument, fieldNumber=0, ins
decs.append(0.0)
airmasses.append(0.0)
exposureTimes.append(0)
mjds.append((refEpoch if (refEpoch is not None) else medianEpoch))
mjds.append((refEpoch if (refEpoch is not None) else medianMJD))
observatories.append(np.array([0, 0, 0]))
identity = wcsfit.IdentityMap()
icrs = wcsfit.SphericalICRS()
Expand Down