{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## DC2 Retrieve Visit-Level Forced Src Photometry and Postage Stamps from Object Catalog.\n",
"
Owner: **Michael Wood-Vasey** ([@wmwv](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@wmwv))\n",
"
Last Verified to Run: **2019-05-08**\n",
"\n",
"### Learning Objectives\n",
"After studying this Notebook you should be able to\n",
"1. Retrieve the catalog of Object-position forced photometry from a given tract, patch, filter, and visit.\n",
"2. Compute the tract, patch for a given RA, Dec and skymap.\n",
"2. Retrieve a postage stamp from the per-visit calibrated image (calexp) for a given object RA, Dec.\n",
"3. Display a full coadd image and full calexp image with the location of an object highlighted.\n",
"\n",
"### See Also\n",
"1. The [Coadd Postage Stamp Notebook](dm_butler_postage_stamps.ipynb) 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.\n",
"\n",
"### Logistics\n",
"1. 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 `repo` below.\n",
"2. If you use `%matplotlib notebook` you 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 see `frame=2; plt.figure(frame)` commands that are effectively noops for `%matplotlib inline` but ensure that the proper frame is being referenced for the interactive windows of `%matplotlib notebook`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"import lsst.daf.persistence as dafPersist\n",
"import lsst.afw.geom as afwGeom\n",
"import lsst.afw.coord as afwCoord\n",
"import lsst.afw.display as afwDisplay\n",
"import lsst.afw.image as afwImage"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# %matplotlib notebook\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"afw_backend = 'matplotlib'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use a simple repo `LSSTDESC/desc-dc2-dm-data` to track where the data repositories are actually stored."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from desc_dc2_dm_data import REPOS\n",
"\n",
"run = '1.2p'\n",
"repo = REPOS[run]\n",
"print(repo)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"butler = dafPersist.Butler(repo)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Reading a Forced-Position Photometry Catalog"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tract = 4850\n",
"patch = '4,5'\n",
"\n",
"filt = 'r'\n",
"partial_data_id = {'tract': tract, 'patch': patch, 'filter': filt}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dataset_type = 'forced_src'\n",
"data_refs = butler.subset(datasetType=dataset_type, dataId=partial_data_id)\n",
"data_ids = [dr.dataId for dr in data_refs\n",
" if butler.datasetExists(datasetType=dataset_type,\n",
" dataId=dr.dataId)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dId = {'filter': 'r', 'visit': 181898, 'raftName': 'R31', 'detectorName': 'S00'}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use this same data_id to get both the calibrated exposure (calexp) and the forced-position photometry ('forced_src')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"forced_src = butler.get(datasetType='forced_src', dataId=dId)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I think in AstroPy rather than AFW Tables. We can easily convert here with `asAstropy()`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cat = forced_src.asAstropy()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(len(forced_src))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# We also could have directly loaded the calibration information using\n",
"calib = butler.get(datasetType='calexp_calib', dataId=dId)\n",
"calib.setThrowOnNegativeFlux(False)\n",
"\n",
"mag, mag_err = calib.getMagnitude(cat['base_PsfFlux_instFlux'],\n",
" cat['base_PsfFlux_instFluxErr'])\n",
"cat['mag'] = mag\n",
"cat['mag_err'] = mag_err\n",
"cat['snr'] = np.abs(cat['base_PsfFlux_instFlux'])/cat['base_PsfFlux_instFluxErr']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(cat)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is per tract, so ~20k sources seems potentially reasonble.\n",
"\n",
"The `id` is the same `id` as in the Object catalog."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def get_photometry_for_id(butler, obj_id, data_id,\n",
" calibration_dataset_type='calexp_calib',\n",
" catalog_dataset_type='forced_src'):\n",
" \"\"\"Return the photometry from a data_id and datasteType for a specified object id.\n",
"\n",
" The motiviation is for reading out photometry from a known Object ID,\n",
" such as that from an Object catalog and reading a forced-photometry catalog.\n",
" But in implementation it will work fine as long as you know the Object ID\n",
" in the catalog datasetType requested.\n",
" \"\"\"\n",
" calib = butler.get(datasetType=calibration_dataset_type, dataId=data_id)\n",
" forced_src = butler.get(datasetType=catalog_dataset_type, dataId=data_id)\n",
"\n",
" matching_id_idx, = np.where(forced_src['id'] == obj_id)\n",
" if len(matching_id_idx) < 1:\n",
" return None\n",
" print(\"Could not find obj_id: %d\" % obj_id)\n",
" \n",
" cat = forced_src.asAstropy()\n",
" calib.setThrowOnNegativeFlux(False)\n",
"\n",
" mag, mag_err = calib.getMagnitude(forced_src['base_PsfFlux_instFlux'],\n",
" forced_src['base_PsfFlux_instFluxErr'])\n",
" cat['mag'] = mag\n",
" cat['mag_err'] = mag_err\n",
" cat['snr'] = np.abs(cat['base_PsfFlux_instFlux'])/cat['base_PsfFlux_instFluxErr']\n",
" \n",
" return cat[matching_id_idx]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"frame = 1\n",
"plt.figure(frame)\n",
"plt.scatter(np.rad2deg(cat['coord_ra']),\n",
" np.rad2deg(cat['coord_dec']))\n",
"plt.xlabel('RA [deg]')\n",
"plt.ylabel('Dec [deg]')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"frame = 2\n",
"plt.figure(frame)\n",
"\n",
"fig, ax = plt.subplots(1, 2, figsize=(12,6))\n",
"ax[0].scatter(cat['mag'], cat['snr'])\n",
"ax[0].set_xlabel('%s mag' % dId['filter'])\n",
"ax[0].set_ylabel('S/N')\n",
"\n",
"ax[1].scatter(cat['mag'], cat['snr'])\n",
"ax[1].set_xlabel('%s mag' % dId['filter'])\n",
"ax[1].set_ylabel('S/N')\n",
"ax[1].set_ylim(0, 10)\n",
"ax[1].axhline(5, color='green', ls='--')\n",
"ax[1].axvline(24.25, color='green', ls='--')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Looks like a reasonable set of magnitudes and uncertainties.\n",
"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.\n",
"\n",
"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)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finding the Visits that Contain a Given Position\n",
"Let's pick one of the Objects and go through all of the data IDs in that filter for the tract, patch."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"good_snr = cat[cat['snr'] > 25]\n",
"obj = good_snr[100]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def get_all_images_for_object_id_from_data_ids(object_id, data_refs, dataset_type='forced_src'):\n",
" \"\"\"Take an Object ID and data_ref list and return matching data_refs that contain that object and exist\n",
" \n",
" Note: This perhaps isn't quite the function you wanted.\n",
" You may have wanted a function that just took an object_id and filter.\n",
" But we'll get there.\n",
" \n",
" Note that we're being a little inefficient and loading each raft, sensor catalog\n",
" and checking to see if it contains the given ID.\n",
" It would likely be more efficient to just load the WCS for the Visit\n",
" and look up the raft, sensor.\n",
" \"\"\"\n",
" matching_data_refs = []\n",
" for data_ref in data_refs:\n",
" # Should add check to make sure it exists\n",
" if not data_ref.datasetExists(datasetType=dataset_type):\n",
" continue\n",
" cat = data_ref.get(datasetType=dataset_type)\n",
" if object_id in cat['id']:\n",
" matching_data_refs.append(data_ref)\n",
"\n",
" return matching_data_refs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"matching_data_refs = get_all_images_for_object_id_from_data_ids(obj['id'], data_refs, dataset_type='forced_src')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data_ids = [dr.dataId for dr in matching_data_refs]\n",
"print(data_ids)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Find the Tract(s), Patch(es) for a given RA, Dec"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We're still missing one step.\n",
"\n",
"How do we get the tract, patch for a given ObjectId or RA, Dec?\n",
"\n",
"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.\n",
"\n",
"To do this we\n",
"1. Get the skymap for the coadd dataset (by default, `deepCoadd`)\n",
"2. Use the skymap object to look up the tract\n",
"3. Use the tract object to look up the patch\n",
"4. Create a partial data Id dict and query the butler for `forced_src` catalogs that match this partial data Id.\n",
"5. Go through each forced_src catalog in that tract, patch and save the ones that match the given Id."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"skymap = butler.get(datasetType='deepCoadd_skyMap')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"help(skymap.findTractPatchList)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ra, dec = obj['coord_ra'], obj['coord_dec']\n",
"# Note the catalog returns coord_ra, coord_dec in RADIANS\n",
"radec = afwGeom.SpherePoint(ra, dec, afwGeom.radians)\n",
"\n",
"tracts_and_patches = skymap.findTractPatchList([radec])\n",
"\n",
"partial_data_ids = [{'tract': tractInfo.getId(), 'patch': '%d,%d' % patch.getIndex()} \\\n",
" for tractInfo, patchList in tracts_and_patches\n",
" for patch in patchList]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"filt = 'r'\n",
"dataset_type = 'forced_src'\n",
"\n",
"data_refs = []\n",
"for partial in partial_data_ids:\n",
" this_data_id = partial.copy()\n",
" this_data_id['filter'] = 'r'\n",
" print(this_data_id)\n",
" \n",
" these_data_refs = butler.subset(datasetType=dataset_type, dataId=this_data_id)\n",
" data_refs.extend(these_data_refs)\n",
"\n",
"data_ids = [dr.dataId for dr in data_refs\n",
" if butler.datasetExists(datasetType=dataset_type,\n",
" dataId=dr.dataId)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"matching_data_refs = get_all_images_for_object_id_from_data_ids(obj['id'], data_refs, dataset_type='forced_src')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data_ids = [dr.dataId for dr in matching_data_refs]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(data_ids)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Here's Some Code to Generate Postage Stamps"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def cutout_ra_dec(butler, ra, dec, data_id, datasetType='calexp',\n",
" cutoutSideLength=51, verbose=False,\n",
" **kwargs):\n",
" \"\"\"\n",
" Produce a cutout from the given image at the given afw SpherePoint radec position.\n",
" \n",
" Parameters\n",
" ----------\n",
" butler: lsst.daf.persistence.Butler\n",
" Servant providing access to a data repository\n",
" ra, dec: Right Ascension, Declination in decimal degrees\n",
" Coordinates of the center of the cutout.\n",
" data_id: Data Id\n",
" datasetType: string ['calexp'] \n",
" cutoutSideLength: float [optional] \n",
" Side of the cutout region in pixels.\n",
" \n",
" Returns\n",
" -------\n",
" MaskedImage\n",
" \"\"\"\n",
" cutoutSize = afwGeom.ExtentI(cutoutSideLength, cutoutSideLength)\n",
"\n",
" radec = afwGeom.SpherePoint(ra, dec, afwGeom.degrees)\n",
"\n",
" calexp = butler.get(datasetType, dataId=data_id)\n",
" xy = afwGeom.PointI(calexp.getWcs().skyToPixel(radec))\n",
" if verbose:\n",
" print(\"Making cutout at (x, y) {xy:} of size ({cutoutSize:}, {cutoutSize:})\".format({'xy': xy, 'cutoutSize': cutoutSize}))\n",
" print(xy, cutoutSize)\n",
"\n",
" bbox = afwGeom.BoxI(xy - cutoutSize//2, cutoutSize)\n",
" \n",
" cutout_image = butler.get(datasetType+'_sub', bbox=bbox, immediate=True, dataId=data_id)\n",
" \n",
" return cutout_image"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def display_cutout_image(butler, ra, dec, data_id,\n",
" vmin=None, vmax=None, label=None,\n",
" frame=None, display=None, backend='matplotlib',\n",
" show=True, saveplot=False, savefits=False,\n",
" datasetType='calexp'):\n",
" \"\"\"\n",
" Display a postage stamp for a given RA, Dec using LSST lsst.afw.display.\n",
" \n",
" Parameters\n",
" ----------\n",
" ra: float [degrees]\n",
" dec: float [degrees]\n",
" backend: string\n",
" Backend can be anything that lsst.afw.display and your configuration supports: \n",
" e.g. matplotlib, ds9, ginga, firefly.\n",
" \n",
" Returns\n",
" -------\n",
" MaskedImage\n",
" \n",
" Notes\n",
" -----\n",
" Parameters are the same as for make_cutout_image, except for the backend.\n",
" You can rely on definitely having the matplotlib backend.\n",
" ds9, ginga, and firefly can be set up but are non-trivial on the scale of a simple Notebook.\n",
" \"\"\"\n",
" cutout_image = cutout_ra_dec(butler, ra, dec, data_id, datasetType='deepCoadd')\n",
" if savefits:\n",
" if isinstance(savefits, str):\n",
" filename = savefits\n",
" else:\n",
" filename = 'postage-stamp.fits'\n",
" cutout_image.writeFits(filename)\n",
" \n",
" if display is None:\n",
" display = afwDisplay.Display(frame=frame, backend=backend)\n",
"\n",
" radec = afwGeom.SpherePoint(ra, dec, afwGeom.degrees)\n",
" xy = cutout_image.getWcs().skyToPixel(radec)\n",
" \n",
" display.mtv(cutout_image)\n",
" display.scale(\"asinh\", \"zscale\")\n",
" display.dot('o', xy.getX(), xy.getY(), ctype='red')\n",
" display.show_colorbar()\n",
" \n",
" return cutout_image"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Finally, Let's Make Some Stamps\n",
"We can now use these data_ids to \n",
"1. Generate postage stamps from the calexps\n",
"2. Extract the photometry from the forced-source catalog by matching to Object Id."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"frame = 3\n",
"plt.figure(frame)\n",
"\n",
"ra, deg = obj['coord_ra'], obj['coord_dec']\n",
"ra_deg, dec_deg = np.rad2deg(ra), np.rad2deg(dec)\n",
"\n",
"for did in data_ids:\n",
" cutout = display_cutout_image(butler, ra_deg, dec_deg, did, datasetType='calexp',\n",
" frame=frame)\n",
" phot = get_photometry_for_id(butler, obj['id'], did,\n",
" catalog_dataset_type='forced_src')\n",
" print(phot['id', 'coord_ra', 'coord_dec', 'mag', 'mag_err'])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for maskName, maskBit in cutout.mask.getMaskPlaneDict().items():\n",
" print('{}: {}'.format(maskName, display.getMaskPlaneColor(maskName)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### For a Bit of Perspective, let's end by looking at the full image"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"did = data_ids[0]\n",
"coadd = butler.get('deepCoadd', dataId=did)\n",
"calexp = butler.get('calexp', dataId=did)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"frame = 4\n",
"plt.figure(frame)\n",
"\n",
"xy = coadd.getWcs().skyToPixel(radec)\n",
"\n",
"display = afwDisplay.Display(frame=frame, backend=afw_backend)\n",
"display.scale(\"asinh\", min=0, max=5)\n",
"\n",
"display.mtv(coadd)\n",
"display.dot('o', xy.getX(), xy.getY(), ctype='red')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"frame = 5 \n",
"plt.figure(frame)\n",
"\n",
"xy = calexp.getWcs().skyToPixel(radec)\n",
"\n",
"display = afwDisplay.Display(frame=frame, backend=afw_backend)\n",
"display.scale(\"asinh\", min=-0.5, max=50)\n",
"display.mtv(calexp)\n",
"display.dot('o', xy.getX(), xy.getY(), ctype='red')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "desc-stack",
"language": "python",
"name": "desc-stack"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}