DC2 Retrieve Visit-Level Forced Src Photometry and Postage Stamps from Object Catalog.#
Owner: Michael Wood-Vasey (@wmwv)
Last Verified to Run: 2019-05-08
Learning Objectives#
After studying this Notebook you should be able to
Retrieve the catalog of Object-position forced photometry from a given tract, patch, filter, and visit.
Compute the tract, patch for a given RA, Dec and skymap.
Retrieve a postage stamp from the per-visit calibrated image (calexp) for a given object RA, Dec.
Display a full coadd image and full calexp image with the location of an object highlighted.
See Also#
The Coadd Postage Stamp Notebook provides a basic introduction for getting postage stamps from the coadd image. The actually generation of postage stamps in this present Notebook is essentially the same. The key new things are figuring out what images to look at.
Logistics#
Meant to be run on NERSC, where these data and environment are available. But beyond having an environment with the stack, the only NERSC-specific aspect is the location of the data. If you have your own set of data elsewhere, you just need to change the
repobelow.If you use
%matplotlib notebookyou can get interactive windows. This works in Jupyter, but not in the full JupyterLab environment, which disallows Javascript. To ensure something that runs in the default JupyterLab environment that we present for DESC, this Notebook uses%matplotlib inline. But you will seeframe=2; plt.figure(frame)commands that are effectively noops for%matplotlib inlinebut ensure that the proper frame is being referenced for the interactive windows of%matplotlib notebook.
import numpy as np
import lsst.daf.persistence as dafPersist
import lsst.afw.geom as afwGeom
import lsst.afw.coord as afwCoord
import lsst.afw.display as afwDisplay
import lsst.afw.image as afwImage
# %matplotlib notebook
%matplotlib inline
import matplotlib.pyplot as plt
afw_backend = 'matplotlib'
We use a simple repo LSSTDESC/desc-dc2-dm-data to track where the data repositories are actually stored.
from desc_dc2_dm_data import REPOS
run = '1.2p'
repo = REPOS[run]
print(repo)
butler = dafPersist.Butler(repo)
Reading a Forced-Position Photometry Catalog#
tract = 4850
patch = '4,5'
filt = 'r'
partial_data_id = {'tract': tract, 'patch': patch, 'filter': filt}
dataset_type = 'forced_src'
data_refs = butler.subset(datasetType=dataset_type, dataId=partial_data_id)
data_ids = [dr.dataId for dr in data_refs
if butler.datasetExists(datasetType=dataset_type,
dataId=dr.dataId)]
The following data_id should be in the list returned by the above query. We explicitly define it here to show how to construct a data_id for the forced source photometry which requires the visit-level information (visit, filter, raftName, detectorName).
dId = {'filter': 'r', 'visit': 181898, 'raftName': 'R31', 'detectorName': 'S00'}
We can use this same data_id to get both the calibrated exposure (calexp) and the forced-position photometry (‘forced_src’)
forced_src = butler.get(datasetType='forced_src', dataId=dId)
I think in AstroPy rather than AFW Tables. We can easily convert here with asAstropy():
cat = forced_src.asAstropy()
print(len(forced_src))
The photometric calibration information is stored in the calibrated exposure. We can extract that and use it to produce calibrated magnitudes for sources with counts > 0. Because these are forced-position photometry based on the location of obejects detected in a deeper coadd, a significant fraction of sources may have estimated counts consistent with zero and thus sometimes negative.
# We also could have directly loaded the calibration information using
calib = butler.get(datasetType='calexp_calib', dataId=dId)
calib.setThrowOnNegativeFlux(False)
mag, mag_err = calib.getMagnitude(cat['base_PsfFlux_instFlux'],
cat['base_PsfFlux_instFluxErr'])
cat['mag'] = mag
cat['mag_err'] = mag_err
cat['snr'] = np.abs(cat['base_PsfFlux_instFlux'])/cat['base_PsfFlux_instFluxErr']
print(cat)
This is per tract, so ~20k sources seems potentially reasonble.
The id is the same id as in the Object catalog.
Let’s write all of this up in a function to extract the photometry for a given object Id. This seems trivial, but it’s worth putting in a function to put the calibration stuff all together.
def get_photometry_for_id(butler, obj_id, data_id,
calibration_dataset_type='calexp_calib',
catalog_dataset_type='forced_src'):
"""Return the photometry from a data_id and datasteType for a specified object id.
The motiviation is for reading out photometry from a known Object ID,
such as that from an Object catalog and reading a forced-photometry catalog.
But in implementation it will work fine as long as you know the Object ID
in the catalog datasetType requested.
"""
calib = butler.get(datasetType=calibration_dataset_type, dataId=data_id)
forced_src = butler.get(datasetType=catalog_dataset_type, dataId=data_id)
matching_id_idx, = np.where(forced_src['id'] == obj_id)
if len(matching_id_idx) < 1:
return None
print("Could not find obj_id: %d" % obj_id)
cat = forced_src.asAstropy()
calib.setThrowOnNegativeFlux(False)
mag, mag_err = calib.getMagnitude(forced_src['base_PsfFlux_instFlux'],
forced_src['base_PsfFlux_instFluxErr'])
cat['mag'] = mag
cat['mag_err'] = mag_err
cat['snr'] = np.abs(cat['base_PsfFlux_instFlux'])/cat['base_PsfFlux_instFluxErr']
return cat[matching_id_idx]
frame = 1
plt.figure(frame)
plt.scatter(np.rad2deg(cat['coord_ra']),
np.rad2deg(cat['coord_dec']))
plt.xlabel('RA [deg]')
plt.ylabel('Dec [deg]')
frame = 2
plt.figure(frame)
fig, ax = plt.subplots(1, 2, figsize=(12,6))
ax[0].scatter(cat['mag'], cat['snr'])
ax[0].set_xlabel('%s mag' % dId['filter'])
ax[0].set_ylabel('S/N')
ax[1].scatter(cat['mag'], cat['snr'])
ax[1].set_xlabel('%s mag' % dId['filter'])
ax[1].set_ylabel('S/N')
ax[1].set_ylim(0, 10)
ax[1].axhline(5, color='green', ls='--')
ax[1].axvline(24.25, color='green', ls='--')
Looks like a reasonable set of magnitudes and uncertainties. The left is the full distribution, while the right plot shows what happens as we go down into the noise floor. Recall that these are the forced photometry on an individual visit based on the detections in the stacked coadd. So we get objects that are below the detection limit in an individual image. The noise is a property of the background, so the smooth continuation of the line even beyond a reasonable single-image limit is due to a constant noise source for a decreasing measured flux.
We get a 5-sigma at r=24.25 mag. So overall on the visit we have a reasonable span of 7 magnitudes (from 17.25 to 24.25 mag).
Finding the Visits that Contain a Given Position#
Let’s pick one of the Objects and go through all of the data IDs in that filter for the tract, patch.
good_snr = cat[cat['snr'] > 25]
obj = good_snr[100]
def get_all_images_for_object_id_from_data_ids(object_id, data_refs, dataset_type='forced_src'):
"""Take an Object ID and data_ref list and return matching data_refs that contain that object and exist
Note: This perhaps isn't quite the function you wanted.
You may have wanted a function that just took an object_id and filter.
But we'll get there.
Note that we're being a little inefficient and loading each raft, sensor catalog
and checking to see if it contains the given ID.
It would likely be more efficient to just load the WCS for the Visit
and look up the raft, sensor.
"""
matching_data_refs = []
for data_ref in data_refs:
# Should add check to make sure it exists
if not data_ref.datasetExists(datasetType=dataset_type):
continue
cat = data_ref.get(datasetType=dataset_type)
if object_id in cat['id']:
matching_data_refs.append(data_ref)
return matching_data_refs
Warning: Actually running the following get_all_images_for_object_id_from_data_ids can take a long time - somehwere between 30 minutes and never.
matching_data_refs = get_all_images_for_object_id_from_data_ids(obj['id'], data_refs, dataset_type='forced_src')
data_ids = [dr.dataId for dr in matching_data_refs]
print(data_ids)
Find the Tract(s), Patch(es) for a given RA, Dec#
We’re still missing one step.
How do we get the tract, patch for a given ObjectId or RA, Dec?
You probably wanted to start with a given ObjectID from the Object (coadd) catalog and then get all of the images that include that object.
To do this we
Get the skymap for the coadd dataset (by default,
deepCoadd)Use the skymap object to look up the tract
Use the tract object to look up the patch
Create a partial data Id dict and query the butler for
forced_srccatalogs that match this partial data Id.Go through each forced_src catalog in that tract, patch and save the ones that match the given Id.
skymap = butler.get(datasetType='deepCoadd_skyMap')
help(skymap.findTractPatchList)
ra, dec = obj['coord_ra'], obj['coord_dec']
# Note the catalog returns coord_ra, coord_dec in RADIANS
radec = afwGeom.SpherePoint(ra, dec, afwGeom.radians)
tracts_and_patches = skymap.findTractPatchList([radec])
partial_data_ids = [{'tract': tractInfo.getId(), 'patch': '%d,%d' % patch.getIndex()} \
for tractInfo, patchList in tracts_and_patches
for patch in patchList]
filt = 'r'
dataset_type = 'forced_src'
data_refs = []
for partial in partial_data_ids:
this_data_id = partial.copy()
this_data_id['filter'] = 'r'
print(this_data_id)
these_data_refs = butler.subset(datasetType=dataset_type, dataId=this_data_id)
data_refs.extend(these_data_refs)
data_ids = [dr.dataId for dr in data_refs
if butler.datasetExists(datasetType=dataset_type,
dataId=dr.dataId)]
matching_data_refs = get_all_images_for_object_id_from_data_ids(obj['id'], data_refs, dataset_type='forced_src')
data_ids = [dr.dataId for dr in matching_data_refs]
print(data_ids)
Here’s Some Code to Generate Postage Stamps#
def cutout_ra_dec(butler, ra, dec, data_id, datasetType='calexp',
cutoutSideLength=51, verbose=False,
**kwargs):
"""
Produce a cutout from the given image at the given afw SpherePoint radec position.
Parameters
----------
butler: lsst.daf.persistence.Butler
Servant providing access to a data repository
ra, dec: Right Ascension, Declination in decimal degrees
Coordinates of the center of the cutout.
data_id: Data Id
datasetType: string ['calexp']
cutoutSideLength: float [optional]
Side of the cutout region in pixels.
Returns
-------
MaskedImage
"""
cutoutSize = afwGeom.ExtentI(cutoutSideLength, cutoutSideLength)
radec = afwGeom.SpherePoint(ra, dec, afwGeom.degrees)
calexp = butler.get(datasetType, dataId=data_id)
xy = afwGeom.PointI(calexp.getWcs().skyToPixel(radec))
if verbose:
print("Making cutout at (x, y) {xy:} of size ({cutoutSize:}, {cutoutSize:})".format({'xy': xy, 'cutoutSize': cutoutSize}))
print(xy, cutoutSize)
bbox = afwGeom.BoxI(xy - cutoutSize//2, cutoutSize)
cutout_image = butler.get(datasetType+'_sub', bbox=bbox, immediate=True, dataId=data_id)
return cutout_image
def display_cutout_image(butler, ra, dec, data_id,
vmin=None, vmax=None, label=None,
frame=None, display=None, backend='matplotlib',
show=True, saveplot=False, savefits=False,
datasetType='calexp'):
"""
Display a postage stamp for a given RA, Dec using LSST lsst.afw.display.
Parameters
----------
ra: float [degrees]
dec: float [degrees]
backend: string
Backend can be anything that lsst.afw.display and your configuration supports:
e.g. matplotlib, ds9, ginga, firefly.
Returns
-------
MaskedImage
Notes
-----
Parameters are the same as for make_cutout_image, except for the backend.
You can rely on definitely having the matplotlib backend.
ds9, ginga, and firefly can be set up but are non-trivial on the scale of a simple Notebook.
"""
cutout_image = cutout_ra_dec(butler, ra, dec, data_id, datasetType='deepCoadd')
if savefits:
if isinstance(savefits, str):
filename = savefits
else:
filename = 'postage-stamp.fits'
cutout_image.writeFits(filename)
if display is None:
display = afwDisplay.Display(frame=frame, backend=backend)
radec = afwGeom.SpherePoint(ra, dec, afwGeom.degrees)
xy = cutout_image.getWcs().skyToPixel(radec)
display.mtv(cutout_image)
display.scale("asinh", "zscale")
display.dot('o', xy.getX(), xy.getY(), ctype='red')
display.show_colorbar()
return cutout_image
Finally, Let’s Make Some Stamps#
We can now use these data_ids to
Generate postage stamps from the calexps
Extract the photometry from the forced-source catalog by matching to Object Id.
frame = 3
plt.figure(frame)
ra, deg = obj['coord_ra'], obj['coord_dec']
ra_deg, dec_deg = np.rad2deg(ra), np.rad2deg(dec)
for did in data_ids:
cutout = display_cutout_image(butler, ra_deg, dec_deg, did, datasetType='calexp',
frame=frame)
phot = get_photometry_for_id(butler, obj['id'], did,
catalog_dataset_type='forced_src')
print(phot['id', 'coord_ra', 'coord_dec', 'mag', 'mag_err'])
for maskName, maskBit in cutout.mask.getMaskPlaneDict().items():
print('{}: {}'.format(maskName, display.getMaskPlaneColor(maskName)))
For a Bit of Perspective, let’s end by looking at the full image#
did = data_ids[0]
coadd = butler.get('deepCoadd', dataId=did)
calexp = butler.get('calexp', dataId=did)
frame = 4
plt.figure(frame)
xy = coadd.getWcs().skyToPixel(radec)
display = afwDisplay.Display(frame=frame, backend=afw_backend)
display.scale("asinh", min=0, max=5)
display.mtv(coadd)
display.dot('o', xy.getX(), xy.getY(), ctype='red')
Note the outlines above showing the pattern of the chip gaps in the visit images that went into the coadd (green), saturated spikes from bright stars (also green), and the identified objects (blue).
frame = 5
plt.figure(frame)
xy = calexp.getWcs().skyToPixel(radec)
display = afwDisplay.Display(frame=frame, backend=afw_backend)
display.scale("asinh", min=-0.5, max=50)
display.mtv(calexp)
display.dot('o', xy.getX(), xy.getY(), ctype='red')