How do I do forced-photometry on a set of RA, Dec

Support
, , ,
#1

If I have a set of RA, Dec (e.g., a 2-column text file with RA, Dec in decimal degrees), what are good ways to perform forced-photometry at those RA, Dec positions on an image that has been astrometrically and photometrically calibrated?

I’m interested in the case of both

  1. single image
  2. subtracted image

I understand how forcedPhotCcd.py is supposed to work if I already have a set of sources in a coadd source catalog (CoaddSrcReferencesTask) or multi-band reference catalog (MultiBandReferencesTask). But I’m not clear on how best to easily create a minimal source catalog just based on a simple list of RA, Dec.

#2

I think your best bet is to use lsst.meas.base.ForcedMeasurementTask directly. That isn’t a CmdLineTask, so you’ll have to write a bit of Python to use it, but it has decent reference documentation.

One big catch is that you will - as you’ve guessed - have to pass it a minimal SourceCatalog containing the positions of the objects you want to measure. You’ll also need to set up a dummy centroider to set the output shape slot from the reference catalog. Unfortunately, we usually use ForcedTransformedCentroid to do that, but it actually transforms pixel-coordinate centroids from the reference to the measurement frame (via both WCSs) rather than looking at the ra,dec coord. We probably need to add a similar plugin that uses the coord to enable ra,dec forced photometry.

1 Like
#3

I wrote a task to perform force photometry on DIA sources, which might be helpful as a starting point to insert your own source loading: forcePhotDiaSources.py. Jim’s caveat about having to create your own SourceCatalog still applies, and I think that’s probably the biggest hurdle.

1 Like
#4

@jbosch @ctslater Thanks for confirming that the somewhat annoying part of setting up the catalog is indeed somewhat annoying.

#5

@ctslater Ah, perfect. Thank you.

Would it be appropriate to add your forcePhotDiaSources.py to the pipe_task/examples or ip_diffim/examples?

#6

@jbosch I’ve successfully written a new transform plugin for the centroids, but how do I mock up/set up the footprints?

#7

This is also a use case for the #suit - we believe that users will want to be able to bring in coordinates from other surveys/catalogs and perform forced photometry on LSST single-epoch and coadded images. This has to work both at an API (Python) level and as a GUI-driven capability.

It is therefore definitely a requirement on the SuperTask design that it can accommodate “input: image, list of coordinates; output: catalog”. So at some point - maybe not in the current 6-month cycle, but soon thereafter - we should prove that out and build such a tool.

#8

If you have footprints in some other coordinate system you can transform them. If you have ellipses for the objects, you can create footprints from those. If you have neither, creating circular footprints the size of the PSF is the best I can come up with.

#9

Great. How do I create footprints from a set of RA, Dec points and given circular radius?

I’m having significant difficulty understanding the afw.detection.Footprint interface from the Python side.

#10

There’s a Footprint constructor that takes an integer position and a radius, so I think something like this should work:

coord = lsst.afw.coord.IcrsCoord(lsst.afw.geom.Point2D(ra, dec), lsst.afw.geom.degrees)
fpCenter = lsst.afw.geom.Point2I(wcs.skyToPixel(coord))
footprint = lsst.afw.detection.Footprint(fpCenter, radius)
1 Like
#11

Thank you, that makes sense. How do I add them to a SourceCatalog? record.set('footprint', footprint) doesn’t work.

#12

There’s a special setter for footprints: record.setFootprint(footprint).

1 Like
#13

I implemented forced-position photometry based on an external catalog in the context of Twinkles (obs_lsstsim). See my notes in the following GitHub repo:

*** Note this relies on a specific branch of meas_base: u/wmwv/radeccentroid to implement the ForcedRaDecCentroid plugin.***

Specifically the README:

Which currently (2017-01-30) consists of:

Twinkles single-image subtractions and forced photometry

Explorations by Michael Wood-Vasey wmwv@pitt.edu

source /global/common/cori/contrib/lsst/lsstDM/setupStack-12_1.sh
setup -j -r /global/homes/w/wmwv/lsst/obs_lsstSim
# SHA1 8d21232
setup -j -r /global/homes/w/wmwv/lsst/meas_base
# SHA1 70cb64f
  1. Run a set of r-band subtractions against the first r-band image for Run 1.1
sh run_imageDifference_run1.1.sh
  1. Generate a catalog from a selected diaSrc for a arbitrarily seleted r-band visit
python make_coord_file.py

Currently visit and outfile are hardcoded. Would like to eventually expand make_coord_file.py to accept argument, such as

python make_coord_file.py --out_file coord_ra_dec.csv --id visit=2210949 filter=r
  1. Run forced-photometry on the calexp images
python twinkles_forcedPhotExternalCatalog.py coord_ra_dec.csv --dataset calexp 
python twinkles_forcedPhotExternalCatalog.py coord_ra_dec.csv --dataset deepDiff_differenceExp

This will take a while, so if you’re running interactively, you may wish to use nohup. E.g.,

nohup python twinkles_forcedPhotExternalCatalog.py coord_ra_dec.csv --dataset deepDiff_differenceExp > phot_all.log 2>&1 < /dev/null &
1 Like
#14

I had a similar need recently, so I will share the solution I came to:

def make_sourceCat_from_astropy(astropy_catalog, exposure):
    schema = afwTable.SourceTable.makeMinimalSchema()
    x_key = schema.addField("centroid_x", type="D")
    y_key = schema.addField("centroid_y", type="D")
    alias = schema.getAliasMap()
    alias.set("slot_Centroid", "centroid")

    sourceCat = afwTable.SourceCatalog(schema)
    for source_rec in astropy_catalog:
        rec = sourceCat.addNew()
        coord = afwCoord.IcrsCoord(source_rec['raICRS']*afwGeom.degrees,
                                   source_rec['decICRS']*afwGeom.degrees)
        rec.setCoord(coord)
        point = exposure.getWcs().skyToPixel(coord)
        fp = afwDet.Footprint(afwGeom.Point2I(point), 6.0)
        rec.setFootprint(fp)
        rec[x_key] = point.getX()
        rec[y_key] = point.getY()
    return sourceCat

def do_forcephot(exposure, sourceCat):
    measurement_config = measBase.ForcedMeasurementConfig()
    measurement_config.plugins.names = ["base_TransformedCentroid", "base_PsfFlux"]
    measurement_config.slots.shape = None

    measurement = measBase.ForcedMeasurementTask(schema, config=measurement_config)

    measCat = measurement.generateMeasCat(exposure, sourceCat, exposure.getWcs())

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