-
Notifications
You must be signed in to change notification settings - Fork 18
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-14720: Implement forced photometry on PVIs in AP pipe. #220
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -30,6 +30,7 @@ | |
import lsst.afw.math as afwMath | ||
import lsst.afw.table as afwTable | ||
from lsst.meas.astrom import AstrometryConfig, AstrometryTask | ||
from lsst.meas.base import ForcedMeasurementTask | ||
from lsst.meas.extensions.astrometryNet import LoadAstrometryNetObjectsTask | ||
from lsst.pipe.tasks.registerImage import RegisterTask | ||
from lsst.meas.algorithms import SourceDetectionTask, SingleGaussianPsf, \ | ||
|
@@ -81,6 +82,8 @@ class ImageDifferenceConfig(pexConfig.Config): | |
doc="Match diaSources with input calexp sources and ref catalog sources") | ||
doMeasurement = pexConfig.Field(dtype=bool, default=True, doc="Measure diaSources?") | ||
doDipoleFitting = pexConfig.Field(dtype=bool, default=True, doc="Measure dipoles using new algorithm?") | ||
doForcedMeasurement = pexConfig.Field(dtype=bool, default=True, | ||
doc="Force photometer diaSource locations on PVI?") | ||
doWriteSubtractedExp = pexConfig.Field(dtype=bool, default=True, doc="Write difference exposure?") | ||
doWriteMatchedExp = pexConfig.Field(dtype=bool, default=False, | ||
doc="Write warped and PSF-matched template coadd exposure?") | ||
|
@@ -131,6 +134,10 @@ class ImageDifferenceConfig(pexConfig.Config): | |
target=DipoleFitTask, | ||
doc="Enable updated dipole fitting method", | ||
) | ||
forcedMeasurement = pexConfig.ConfigurableField( | ||
target=ForcedMeasurementTask, | ||
doc="Subtask to force photometer PVI at diaSource location.", | ||
) | ||
getTemplate = pexConfig.ConfigurableField( | ||
target=GetCoaddAsTemplateTask, | ||
doc="Subtask to retrieve template exposure and sources", | ||
|
@@ -189,6 +196,12 @@ def setDefaults(self): | |
# after the user has set doPreConvolve. | ||
self.measurement.algorithms.names.add('base_PeakLikelihoodFlux') | ||
|
||
self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"] | ||
self.forcedMeasurement.copyColumns = { | ||
"id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"} | ||
self.forcedMeasurement.slots.centroid = "base_TransformedCentroid" | ||
self.forcedMeasurement.slots.shape = None | ||
|
||
# For shuffling the control sample | ||
random.seed(self.controlRandomSeed) | ||
|
||
|
@@ -249,6 +262,12 @@ def __init__(self, butler=None, **kwargs): | |
if self.config.doMeasurement: | ||
self.makeSubtask("measurement", schema=self.schema, | ||
algMetadata=self.algMetadata) | ||
if self.config.doForcedMeasurement: | ||
self.totFluxKey = self.schema.addField( | ||
"totFlux", "D", "Forced flux measured on the PVI") | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. the name of these keys is not very descriptive, and does not seem to match the documentation exactly. I assume tot means total? How does this "total" relate to the forced flux measured? Does the forced flux give a total in some way that is clearer if you know about PVI (it may just be my ignorance). Consider changing the name to something more descriptive (or just like pviFlux etc), or changing the documentation string to clarify. I would say add a comment here, which may not be a bad idea, but that will not be documentation to someone who later just sees the schema and not what produced it. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. |
||
self.totFluxErrKey = self.schema.addField( | ||
"totFluxErr", "D", "Forced flux error measured on the PVI") | ||
self.makeSubtask("forcedMeasurement", refSchema=self.schema) | ||
if self.config.doMatchSources: | ||
self.schema.addField("refMatchId", "L", "unique id of reference catalog match") | ||
self.schema.addField("srcMatchId", "L", "unique id of source match") | ||
|
@@ -603,6 +622,18 @@ def runDataRef(self, sensorRef, templateIdList=None): | |
else: | ||
self.measurement.run(diaSources, subtractedExposure, exposure) | ||
|
||
if self.config.doForcedMeasurement: | ||
# Run forced psf photometry on the PVI at the diaSource locations. | ||
# Copy the measured flux and error into the diaSource. | ||
forcedSources = self.forcedMeasurement.generateMeasCat( | ||
exposure, diaSources, subtractedExposure.getWcs()) | ||
self.forcedMeasurement.run(forcedSources, exposure, diaSources, subtractedExposure.getWcs()) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. If this code is run in a performance critical section of code, consider doing the copying of records with a schema mapper instead of a loop. For this review it does not really matter to me either way, but wanted to bring it up in case it made a difference. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I've tried working with schema mappers in the past and found them obtuse and unintuitive. If you can provide an example of this type of operation using a schema mapper that would be great. This may be asking too much but if you could provide a code snippet to use a schema mapper in this case I would be all for it. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Hey @morriscb, does this help? import random
import lsst.afw.table as afwTable
import lsst.afw.geom as afwGeom
N_SOURCES = 1
def make_catalogs(num_sources):
# Create some placeholder catalogs with random data.
diaSourceSchema = afwTable.SourceTable.makeMinimalSchema()
totFluxKey = diaSourceSchema.addField("totFlux", doc="xx", type="F")
totFluxErrKey = diaSourceSchema.addField("totFluxErr", doc="xx", type="F")
diaSources = afwTable.SourceCatalog(diaSourceSchema)
diaSources.table.preallocate(num_sources) # required for columnar operations
for i in range(num_sources):
rec = diaSources.addNew()
rec.set("coord_ra", random.random() * afwGeom.degrees)
rec.set("coord_dec", random.random() * afwGeom.degrees)
forcedSourceSchema = afwTable.SourceTable.makeMinimalSchema()
base_PsfFlux_fluxKey = forcedSourceSchema.addField("base_PsfFlux_flux", doc="xx", type="F")
base_PsfFlux_fluxErrKey = forcedSourceSchema.addField("base_PsfFlux_fluxErr", doc="xx", type="F")
forcedSources = afwTable.SourceCatalog(forcedSourceSchema)
forcedSources.table.preallocate(num_sources) # required for columnar operations
for i in range(num_sources):
rec = forcedSources.addNew()
rec.set(base_PsfFlux_fluxKey, random.random())
rec.set(base_PsfFlux_fluxErrKey, random.random())
return diaSources, forcedSources
def copy_mapper(diaSources, forcedSources):
# Copy PsfFlux from forceSources to diaSources with a mapper
mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flux")[0], "totFlux", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_fluxErr")[0], "totFluxErr", True)
for forced, dia in zip(forcedSources, diaSources):
dia.assign(forced, mapper)
return diaSources
def copy_brute(diaSources, forcedSources):
# Copy PsfFlux from forceSources to diaSources with "brute force" assignment
for forced, dia in zip(forcedSources, diaSources):
dia.set("totFlux", forced["base_PsfFlux_flux"])
dia.set("totFluxErr", forced["base_PsfFlux_fluxErr"])
return diaSources
def copy_column(diaSources, forcedSources):
# Copy PsfFlux from forceSources to diaSources as a column
# NB requires that the tables be contiguous in memory (ie, the preallocate step, above).
diaSources["totFlux"] = forcedSources["base_PsfFlux_flux"]
diaSources["totFluxErr"] = forcedSources["base_PsfFlux_fluxErr"]
return diaSources
if __name__ == "__main__":
# Just to check the copies worked: all entries should have both positions & fluxes.
print(copy_mapper(*make_catalogs(N_SOURCES)))
print(copy_brute(*make_catalogs(N_SOURCES)))
print(copy_column(*make_catalogs(N_SOURCES))) For what it's worth, using the mapper seems to be a bit faster than doing things the “brute force” way, although columnar operations are faster still... but require your tables to be contiguous in memory. I don't know if that's guaranteed in your situation.
I also give no guarantees that a loop with There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Awesome, thanks John. From my previous experience the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Committed and tested the new schema mapper code. @natelust if you could please take one last look and let me know if there's anything else needed. Thanks! |
||
mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema) | ||
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flux")[0], "totFlux", True) | ||
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_fluxErr")[0], "totFluxErr", True) | ||
for diaSource, forcedSource in zip(diaSources, forcedSources): | ||
diaSource.assign(forcedSource, mapper) | ||
|
||
# Match with the calexp sources if possible | ||
if self.config.doMatchSources: | ||
if sensorRef.datasetExists("src"): | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Just to verify these are the only plugins you want right? I assume it is doing what you want since you must have tested it, but sometimes people get bitten by =, |=, +=, etc for this datatype so I wanted to double check.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Correct, I only want these pluggins to be run for simplicity.