Forced photometry on DP0.2

Support Data Preview 0
#1

Hi all,

I’ve been trying to run the LSST forced photometry pipelines on specific ra/dec in images from DP0.2 on the RSP, but no luck so far. It would be great if anyone can point me towards any existing code and/or tutorials on how to do this. I went through the tutorials on the RSP but they do not cover doing forced photometry yourself.

This is what I have so far:

I’ve come to the conclusion that I likely need to use ForcedMeasurementTask (ForcedMeasurementTask — LSST Science Pipelines). The run method needs:
exposure: the DP0.2 image
refWcs: Coordinate system, I can get that from the simulated DP0.2 image
measCat: a SourceCatalog that will collect the measurement results
refCat: a SourceCatalog that contains the information of where I want to do forced photometry.

I’m having trouble with generating measCat and refCat correctly.

For the measCat variable, I currently do:

schema = afwTable.SourceTable.makeMinimalSchema()
measurement_config = measBase.ForcedMeasurementConfig()
measurement = measBase.ForcedMeasurementTask(schema, config=measurement_config)
measCat = measurement.generateMeasCat(exposure, sourceCat, exposure.getWcs())

which results in this error:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[40], line 6
      4 schema = afwTable.SourceTable.makeMinimalSchema()
      5 measurement_config = measBase.ForcedMeasurementConfig()
----> 6 measurement = measBase.ForcedMeasurementTask(schema, config=measurement_config)
      7 measCat = measurement.generateMeasCat(exposure, sourceCat, exposure.getWcs())

File [/opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe](https://data.lsst.cloud/opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe)+a4f2ea4efa[/python/lsst/meas/base/forcedMeasurement.py:254](https://data.lsst.cloud/python/lsst/meas/base/forcedMeasurement.py#line=253), in ForcedMeasurementTask.__init__(self, refSchema, algMetadata, **kwds)
    252 self.config.slots.setupSchema(self.mapper.editOutputSchema())
    253 for refName, targetName in self.config.copyColumns.items():
--> 254     refItem = refSchema.find(refName)
    255     self.mapper.addMapping(refItem.key, targetName)
    256 self.config.slots.setupSchema(self.mapper.editOutputSchema())

KeyError: "Field with name 'deblend_nChild' not found"

And I’m not sure how to proceed from this.

For the refCat variable, I currently create a SourceCatalog and add my ra/dec:

source = pd.Series({'ra': 56.5498060077, 'dec':-36.4269781414}) #my ra/dec where I want to do forced photometry

schema = afwTable.SourceTable.makeMinimalSchema()
sourceCat = afwTable.SourceCatalog(schema)
coord = geom.SpherePoint(source['ra']*geom.degrees, source['dec']*geom.degrees)
sourceRec.setCoord(coord)

The next step seems to be that I need to add a Footprint to this, but I can’t seem to get that to work (the documentation on this is not super clear). I currently do:

coord = geom.SpherePoint(source['ra']*geom.degrees, 
                             source['dec']*geom.degrees)
point = exposure.getWcs().skyToPixel(coord)
fpCenter = geom.Point2I(point)
radius = 6
footprint = afwDetection.Footprint(fpCenter, radius)

which results in the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[17], line 6
      4 fpCenter = geom.Point2I(point)
      5 radius = 6
----> 6 footprint = afwDetection.Footprint(fpCenter, radius)

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. lsst.afw.detection.Footprint(inputSpans: lsst.afw.geom.SpanSet, region: lsst.geom.Box2I = Box2I(corner=Point2I(0, 0), dimensions=Extent2I(0, 0)))
    2. lsst.afw.detection.Footprint(inputSpans: lsst.afw.geom.SpanSet, peakSchema: lsst.afw.table.Schema, region: lsst.geom.Box2I = Box2I(corner=Point2I(0, 0), dimensions=Extent2I(0, 0)))
    3. lsst.afw.detection.Footprint(arg0: lsst.afw.detection.Footprint)
    4. lsst.afw.detection.Footprint()

Invoked with: Point2I(2039, 2028), 6

I can get it to run with a Box2I instead, but then I’m confused on how to define the aperture (radius), and I’m not sure if this is the recommended way to do it.

Any advice on how to correctly set up the refCat and measCat variables is greatly appreciated. Also please let me know if I’m heading towards the wrong solution entirely.

FYI I have found this thread (How do I do forced-photometry on a set of RA, Dec - #13 by mwv), which covers the same question, but the last answer is from 7 years ago and the code shown there doesn’t work anymore (presumably since LSST pipelines seems to have changed since then).

Thanks!

Cheers,
Tobias

2 Likes
#2

Hi @TobiasGeron, just a quick note to say thanks for posting this. We’re following-up internally, but in the meantime everyone should please feel free to chime in on this thread.

#3

I’m afraid ForcedMeasurementTask is both the best thing we’ve got for this important problem and also not quite suited for it; it makes some assumptions (especially in its default configuration) about the context in which its being run that are probably only valid in the handful of ways we run it in our own pipelines.

That said, I think you can get a lot farther - and maybe all the way to something working - if you modify the config object passed to the task at construction, setting:

# deblend_nChild is being requested because it's in copyColumns.
config.copyColumns  = {}
# Run only a very minimal set of measurements to avoid complications.
config.plugins.names = [
    "base_PixelFlags",
    "base_TransformedCentroid",
    "base_PsfFlux",
]
# "replace with noise" is how we use deblender outputs; disable that
# (measurements will be affected by neighbors).  This should mean
# that we don't need Footprints, which are for deblending, not
# aperture sizes (those are set the configuration for the aperture flux
# plugin, which we've disabled for now anyway).
config.doReplaceWithNoise = False
#4

Hi Jim, thanks for your reply, I appreciate it!

When I add the changes you recommend to create measCat, I get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[21], line 23
     15 # "replace with noise" is how we use deblender outputs; disable that
     16 # (measurements will be affected by neighbors).  This should mean
     17 # that we don't need Footprints, which are for deblending, not
     18 # aperture sizes (those are set the configuration for the aperture flux
     19 # plugin, which we've disabled for now anyway).
     20 measurement_config.doReplaceWithNoise = False
---> 23 measurement = measBase.ForcedMeasurementTask(schema, config=measurement_config)
     24 measCat = measurement.generateMeasCat(exposure, sourceCat, exposure.getWcs())

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/forcedMeasurement.py:249, in ForcedMeasurementTask.__init__(self, refSchema, algMetadata, **kwds)
    248 def __init__(self, refSchema, algMetadata=None, **kwds):
--> 249     super(ForcedMeasurementTask, self).__init__(algMetadata=algMetadata, **kwds)
    250     self.mapper = lsst.afw.table.SchemaMapper(refSchema)
    251     self.mapper.addMinimalSchema(lsst.afw.table.SourceTable.makeMinimalSchema(), False)

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/baseMeasurement.py:253, in BaseMeasurementTask.__init__(self, algMetadata, **kwds)
    252 def __init__(self, algMetadata=None, **kwds):
--> 253     super(BaseMeasurementTask, self).__init__(**kwds)
    254     self.plugins = PluginMap()
    255     self.undeblendedPlugins = PluginMap()

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/pipe_base/g6fe576e57d+11c5857c39/python/lsst/pipe/base/task.py:197, in Task.__init__(self, config, name, parentTask, log)
    195 self.log: lsst.utils.logging.LsstLogAdapter = lsst.utils.logging.getLogger(loggerName)
    196 self.config: Config = config
--> 197 self.config.validate()
    198 if lsstDebug:
    199     self._display = lsstDebug.Info(self.__module__).display

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/baseMeasurement.py:205, in BaseMeasurementConfig.validate(self)
    203     raise ValueError("source centroid slot algorithm is not being run.")
    204 if self.slots.shape is not None and self.slots.shape not in self.plugins.names:
--> 205     raise ValueError("source shape slot algorithm '%s' is not being run." % self.slots.shape)
    206 for slot in (self.slots.psfFlux, self.slots.apFlux, self.slots.modelFlux,
    207              self.slots.gaussianFlux, self.slots.calibFlux):
    208     if slot is not None:

ValueError: source shape slot algorithm 'base_TransformedShape' is not being run.

That error goes away after I add:

measurement_config.slots.shape = None

but I’m not sure if that is the correct way of dealing with that error?

Either way, if I then not add any footprint to refCat (so refCat is generated with the code I shared in the previous post) and go ahead and run the Forced Photometry Task:

measurement.run(measCat, exposure, sourceCat, exposure.getWcs()),

I get the following error:

Traceback (most recent call last):
  File "/opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/baseMeasurement.py", line 387, in doMeasurement
    plugin.measure(measRecord, *args, **kwds)
  File "/opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/plugins.py", line 842, in measure
    measRecord.set(self.centroidKey, refRecord.getCentroid())
                                     ^^^^^^^^^^^^^^^^^^^^^^^
lsst.pex.exceptions.wrappers.LogicError: 
  File "include/lsst/afw/table/BaseRecord.h", line 113, in const typename lsst::afw::table::Field<T>::Element* lsst::afw::table::BaseRecord::getElement(const lsst::afw::table::Key<U>&) const [with T = double; typename lsst::afw::table::Field<T>::Element = double]
    Key is not valid (if this is a SourceRecord, make sure slot aliases have been set up). {0}
lsst::pex::exceptions::LogicError: 'Key is not valid (if this is a SourceRecord, make sure slot aliases have been set up).'

---------------------------------------------------------------------------
LogicError                                Traceback (most recent call last)
File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/baseMeasurement.py:387, in BaseMeasurementTask.doMeasurement(self, plugin, measRecord, *args, **kwds)
    386 try:
--> 387     plugin.measure(measRecord, *args, **kwds)
    388 except FATAL_EXCEPTIONS:

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/plugins.py:842, in ForcedTransformedCentroidPlugin.measure(self, measRecord, exposure, refRecord, refWcs)
    841 else:
--> 842     measRecord.set(self.centroidKey, refRecord.getCentroid())
    843 if self.flagKey is not None:

LogicError: 
  File "include/lsst/afw/table/BaseRecord.h", line 113, in const typename lsst::afw::table::Field<T>::Element* lsst::afw::table::BaseRecord::getElement(const lsst::afw::table::Key<U>&) const [with T = double; typename lsst::afw::table::Field<T>::Element = double]
    Key is not valid (if this is a SourceRecord, make sure slot aliases have been set up). {0}
lsst::pex::exceptions::LogicError: 'Key is not valid (if this is a SourceRecord, make sure slot aliases have been set up).'


During handling of the above exception, another exception occurred:

NotImplementedError                       Traceback (most recent call last)
Cell In[38], line 1
----> 1 measurement.run(measCat, exposure, sourceCat, exposure.getWcs())

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/forcedMeasurement.py:372, in ForcedMeasurementTask.run(self, measCat, exposure, refCat, refWcs, exposureId, beginOrder, endOrder)
    370 # Then process the parent record
    371 noiseReplacer.insertSource(refParentRecord.getId())
--> 372 self.callMeasure(measParentRecord, exposure, refParentRecord, refWcs,
    373                  beginOrder=beginOrder, endOrder=endOrder)
    374 self.callMeasureN(measParentCat[parentIdx:parentIdx+1], exposure,
    375                   refParentCat[parentIdx:parentIdx+1],
    376                   beginOrder=beginOrder, endOrder=endOrder)
    377 # Measure all the children simultaneously

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/baseMeasurement.py:357, in BaseMeasurementTask.callMeasure(self, measRecord, *args, **kwds)
    355 if endOrder is not None and plugin.getExecutionOrder() >= endOrder:
    356     break
--> 357 self.doMeasurement(plugin, measRecord, *args, **kwds)

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/baseMeasurement.py:399, in BaseMeasurementTask.doMeasurement(self, plugin, measRecord, *args, **kwds)
    395 except Exception as error:
    396     self.log.getChild(plugin.name).debug(
    397         "Exception in %s.measure on record %s: %s",
    398         plugin.name, measRecord.getId(), error)
--> 399     plugin.fail(measRecord)

File /opt/lsst/software/stack/stack/miniconda3-py38_4.9.2-8.0.0/Linux64/meas_base/g212a7c68fe+a4f2ea4efa/python/lsst/meas/base/pluginsBase.py:166, in BasePlugin.fail(self, measRecord, error)
    162 traceback.print_exc()
    163 message = ("The algorithm '%s' thinks it cannot fail, but it did; "
    164            "please report this as a bug (the full traceback is above)."
    165            % (self.__class__.__name__,))
--> 166 raise NotImplementedError(message)

NotImplementedError: The algorithm 'ForcedTransformedCentroidPlugin' thinks it cannot fail, but it did; please report this as a bug (the full traceback is above).

And I’m not sure how to solve this. I tried to set up slot aliases but no luck so far. Any advice is appreciated :slight_smile:

#5

Melissa and I worked on this issue this week and wanted to share our progress. We ended up doing avoiding the previous error by adding the base_TransformedShape plugin. I also had some errors related to slot aliases on the refCat variable, but was able to resolve that by adding shape columns. I put everything I had into the function below, which now takes a DP0.2 exposure and ra/dec pair and performs forced photometry:

def forced_photometry(exposure, coords):
    '''
    Does forced photometry on a ra/dec on an exposure.
    Exposure currently should be a DP0.2 ExposureF. 
    Coords should be a tuple of ra,dec, or list of tuples with ra,decs.  Ra and dec in degrees.

    We're using the ForcedMeasurementTask (https://pipelines.lsst.io/py-api/lsst.meas.base.ForcedMeasurementTask.html#lsst.meas.base.ForcedMeasurementTask)

    The run method needs the following arguments:  
    measCat: seems to be an empty catalogue that the function will populate? Not sure.  
    exposure: the actual image  
    refCat: information of the sources. RA/Dec etc.  
    refWcs:  Wcs information, obtained from exposure.
    '''

    if type(coords) == tuple: 
        coords = [coords]

    
    ### Step 1: Create and configure catalogue with all the sources to do forced photometry on (refCat) ###
 
    schema = afwTable.SourceTable.makeMinimalSchema() # Create a SourceTable Schema, which can be transformed to a Source Catalog

    alias = schema.getAliasMap() 
    x_key = schema.addField("centroid_x", type="D") # Add centroids columns to schema, ForcedMeasurementTask needs these columns to run
    y_key = schema.addField("centroid_y", type="D")
    alias.set("slot_Centroid", "centroid")

    xx_key = schema.addField("shape_xx", type="D") # Add shape columns schema, ForcedMeasurementTask needs these columns to run
    yy_key = schema.addField("shape_yy", type="D")
    xy_key = schema.addField("shape_xy", type="D")
    alias.set("slot_Shape", "shape")

    sourceCat = afwTable.SourceCatalog(schema) #Change schema into sourceCatalog
    
    for i in range(len(coords)): # Add a new source record for every entry
        ra, dec = coords[i][0], coords[i][1]
        sourceRec = sourceCat.addNew() 
        coord = geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)
        sourceRec.setCoord(coord) 
        
        point = exposure.getWcs().skyToPixel(coord)
        sourceRec[x_key] = point.getX()
        sourceRec[y_key] = point.getY()
    


    ### Step 2: Configure forced photometry and create measCat ###

    #schema = afwTable.SourceTable.makeMinimalSchema()
    measurement_config = measBase.ForcedMeasurementConfig()
    measurement_config.copyColumns  = {}  # deblend_nChild is being requested because it's in copyColumns.
    measurement_config.plugins.names = [ # Run only a very minimal set of measurements to avoid complications. Note to self: can also add 'base_CircularApertureFlux' as plugin, if we want circular aperture rather than PSF. 
        "base_PixelFlags",
        "base_TransformedCentroid",
        "base_PsfFlux",
        "base_TransformedShape" 
    ] 
    measurement_config.doReplaceWithNoise = False #Disable so that we don't need footprints
    #measurement_config.slots.shape = None #To avoid ValueError: "source shape slot algorithm 'base_TransformedShape' is not being run.", but can also just run that plugin instead
    measurement = measBase.ForcedMeasurementTask(schema, config=measurement_config)
    measCat = measurement.generateMeasCat(exposure, sourceCat, exposure.getWcs()) #Create measCat

    
    ### Step 3: Run forced photometry ###
    measurement.run(measCat, exposure, sourceCat, exposure.getWcs())

    return measCat.asAstropy().to_pandas()

The flux, its error and area of the PSF are in the following columns: base_PsfFlux_instFlux, base_PsfFlux_instFluxErr, base_PsfFlux_area.

Any additional thoughts are obviously welcome, but the code above seems to work as expected. I ran SourceDetectionTask and SingleFrameMeasurementTask on the same exposure and then used the forced photometry function above on the same coordinates, and the psfFlux output (as well as its error and the PSF area) is virtually the same in both cases.

Thanks for the help!

1 Like
#6

Thank you Tobias!

And now we have updated the DP0.2 tutorial notebook 05, “Source Detection and Measurement” to include a demonstration of forced photometry on a calexp in Section 5, just like you have figured out how to do.

This tutorial notebook can be found:

1 Like