Skip to content

Commit

Permalink
Accept references from multiple patches
Browse files Browse the repository at this point in the history 
ForcedPhotCcd assumed each CCD only overlappend one patch, fix.
  • Loading branch information
@natelust
natelust committed Feb 22, 2021
1 parent dbe88b8 commit 01ec831
Show file tree
Hide file tree
Showing 2 changed files with 66 additions and 43 deletions.
106 changes: 64 additions & 42 deletions python/lsst/meas/base/forcedPhotCcd.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>.

import collections
from itertools import chain

import lsst.pex.config
import lsst.pex.exceptions
Expand All @@ -36,6 +35,7 @@
import lsst.pipe.base.connectionTypes as cT

import lsst.pipe.base as pipeBase
from lsst.skymap import BaseSkyMap

from .references import MultiBandReferencesTask
from .forcedMeasurement import ForcedMeasurementTask
Expand Down Expand Up @@ -174,11 +174,11 @@ class ForcedPhotCcdConnections(PipelineTaskConnections,
dimensions=["skymap", "tract", "patch"],
multiple=True
)
refWcs = cT.Input(
doc="Reference world coordinate system.",
name="{inputCoaddName}Coadd.wcs",
storageClass="Wcs",
dimensions=["band", "skymap", "tract", "patch"],
skyMap = cT.Input(
doc="SkyMap dataset that defines the coordinate system of the reference catalog.",
name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
storageClass="SkyMap",
dimensions=["skymap"],
)
measCat = cT.Output(
doc="Output forced photometry catalog.",
Expand Down Expand Up @@ -336,7 +336,14 @@ def __init__(self, butler=None, refSchema=None, initInputs=None, **kwds):

def runQuantum(self, butlerQC, inputRefs, outputRefs):
inputs = butlerQC.get(inputRefs)
inputs['refCat'] = self.filterReferences(inputs['exposure'], inputs['refCat'], inputs['refWcs'])

tract = butlerQC.quantum.dataId['tract']
skyMap = inputs.pop("skyMap")
inputs['refWcs'] = skyMap[tract].getWcs()

inputs['refCat'] = self.mergeAndFilterReferences(inputs['exposure'], inputs['refCat'],
inputs['refWcs'])

inputs['measCat'], inputs['exposureId'] = self.generateMeasCat(inputRefs.exposure.dataId,
inputs['exposure'],
inputs['refCat'], inputs['refWcs'],
Expand All @@ -346,16 +353,16 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
outputs = self.run(**inputs)
butlerQC.put(outputs, outputRefs)

def filterReferences(self, exposure, refCat, refWcs):
def mergeAndFilterReferences(self, exposure, refCats, refWcs):
"""Filter reference catalog so that all sources are within the
boundaries of the exposure.
Parameters
----------
exposure : `lsst.afw.image.exposure.Exposure`
Exposure to generate the catalog for.
refCat : `lsst.afw.table.SourceCatalog`
Catalog of shapes and positions at which to force photometry.
refCats : sequence of `lsst.afw.table.SourceCatalog`
Catalogs of shapes and positions at which to force photometry.
refWcs : `lsst.afw.image.SkyWcs`
Reference world coordinate system.
Expand Down Expand Up @@ -385,37 +392,52 @@ def filterReferences(self, exposure, refCat, refWcs):
expPolygon = lsst.sphgeom.ConvexPolygon(expSkyCorners)

# Step 2: Filter out reference catalog sources that are
# not contained within the exposure boundaries.
# refCat will be a list of at least one element
sources = type(refCat[0])(refCat[0].table)
for record in chain(*refCat):
if expPolygon.contains(record.getCoord().getVector()):
sources.append(record)
refCatIdDict = {ref.getId(): ref.getParent() for ref in sources}

# Step 3: Cull sources that do not have their parent
# source in the filtered catalog. Save two copies of each
# source.
refSources = type(refCat[0])(refCat[0].table)
for record in chain(*refCat):
if expPolygon.contains(record.getCoord().getVector()):
recordId = record.getId()
topId = recordId
while (topId > 0):
if topId in refCatIdDict:
topId = refCatIdDict[topId]
else:
break
if topId == 0:
refSources.append(record)

# Step 4: Transform source footprints from the reference
# coordinates to the exposure coordinates.
for refRecord in refSources:
refRecord.setFootprint(refRecord.getFootprint().transform(refWcs,
expWcs, expRegion))
# Step 5: Replace reference catalog with filtered source list.
return refSources
# not contained within the exposure boundaries, or whose
# parents are not within the exposure boundaries. Note
# that within a single input refCat, the parents always
# appear before the children.
mergedRefCat = lsst.afw.table.SourceCatalog(refCats[0].table)
for refCat in refCats:
containedIds = {0} # zero as a parent ID means "this is a parent"
for record in refCat:
if expPolygon.contains(record.getCoord().getVector()) and record.getParent() in containedIds:
record.setFootprint(record.getFootprint().transform(refWcs, expWcs, expRegion))
mergedRefCat.append(record)
containedIds.add(record.getId())
mergedRefCat.sort(lsst.afw.table.SourceTable.getParentKey())
return mergedRefCat

# # not contained within the exposure boundaries.
# # refCat will be a list of at least one element
# sources = type(refCat[0])(refCat[0].table)
# for record in chain(*refCat):
# if expPolygon.contains(record.getCoord().getVector()):
# sources.append(record)
# refCatIdDict = {ref.getId(): ref.getParent() for ref in sources}

# # Step 3: Cull sources that do not have their parent
# # source in the filtered catalog. Save two copies of each
# # source.
# refSources = type(refCat[0])(refCat[0].table)
# for record in chain(*refCat):
# if expPolygon.contains(record.getCoord().getVector()):
# recordId = record.getId()
# topId = recordId
# while (topId > 0):
# if topId in refCatIdDict:
# topId = refCatIdDict[topId]
# else:
# break
# if topId == 0:
# refSources.append(record)

# # Step 4: Transform source footprints from the reference
# # coordinates to the exposure coordinates.
# for refRecord in refSources:
# refRecord.setFootprint(refRecord.getFootprint().transform(refWcs,
# expWcs, expRegion))
# # Step 5: Replace reference catalog with filtered source list.
# return refSources

def generateMeasCat(self, exposureDataId, exposure, refCat, refWcs, idPackerName):
"""Generate a measurement catalog for Gen3.
Expand Down Expand Up @@ -619,7 +641,7 @@ def getExposure(self, dataRef):
dataRef : `lsst.daf.persistence.ButlerDataRef`
Butler data reference.
"""
exposure = ForcedPhotImageTask.getExposure(self, dataRef)
exposure = dataRef.get(self.dataPrefix + "calexp", immediate=True)

if self.config.doApplyExternalPhotoCalib:
source = f"{self.config.externalPhotoCalibName}_photoCalib"
Expand Down
3 changes: 2 additions & 1 deletion python/lsst/meas/base/forcedPhotCoadd.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ class ForcedPhotCoaddConnections(pipeBase.PipelineTaskConnections,
name="{inputCoaddName}Coadd.wcs",
storageClass="Wcs",
dimensions=["band", "skymap", "tract", "patch"],
)
) # used in place of a skymap wcs because of DM-28880
measCat = pipeBase.connectionTypes.Output(
doc="Output forced photometry catalog.",
name="{outputCoaddName}Coadd_forced_src",
Expand Down Expand Up @@ -196,6 +196,7 @@ def __init__(self, butler=None, refSchema=None, initInputs=None, **kwds):

def runQuantum(self, butlerQC, inputRefs, outputRefs):
inputs = butlerQC.get(inputRefs)

refCatInBand = inputs.pop('refCatInBand')
inputs['measCat'], inputs['exposureId'] = self.generateMeasCat(inputRefs.exposure.dataId,
inputs['exposure'],
Expand Down

0 comments on commit 01ec831

Please sign in to comment.