Owner(s): Bryce Kalmbach (@jbkalmbach)
Last Verified to Run: 2021-04-04
Verified Stack Release: v21.0.0
This notebook shows
After working through this lesson you should be able to:
This notebook is intended to be runnable on lsst-lsp-stable.ncsa.illinois.edu
from a local git clone of https://github.com/LSSTScienceCollaborations/StackClubCourse.
This notebook uses methods from these other Stack Club notebooks:
Low-Surface Brightness Source Detection
as well as previous notebooks in the Stack Club Course.
The image data in this notebook is DECam data from the HiTS survey processed by Meredith Rawls (@mrawls) (original dataset location: /project/mrawls/hits2015/rerun/cw_2020_04
).
You can find the Stack version that this notebook is running by using eups list -s
on the terminal command line:
# What version of the Stack am I using?
! echo $HOSTNAME
! eups list lsst_distrib -s
nb-kadrlica-r21-0-0 21.0.0+973e4c9e85 current v21_0_0 setup
We will need the following packages
import lsst.daf.persistence as dafPersist
import lsst.afw.image as afwImage
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lsst.geom
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.time import Time
# Enable interactive, in-notebook plotting with matplotlib
%matplotlib ipympl
# Some other options if the above doesn't work
#%matplotlib widget
#%matplotlib inline
import lsst.afw.display as afw_display
afw_display.setDefaultBackend('matplotlib')
Our first task is (not surprisingly) to point the butler at the location of the repository we want to use.
data_dir = '/project/stack-club/decam_hits_2015_subset/'
butler = dafPersist.Butler(data_dir)
We happen to know a priori that there are difference images in this dataset, but what are the actual datasets called when we access them through the butler? We can remind ourselves with getDatasetTypes
.
data_types = butler.getDatasetTypes()
diff_data_types = [x for x in data_types if x.startswith('deepDiff') ]
list(np.sort(list(diff_data_types)))
['deepDiff_config', 'deepDiff_config_filename', 'deepDiff_diaSrc', 'deepDiff_diaSrc_filename', 'deepDiff_diaSrc_len', 'deepDiff_diaSrc_md', 'deepDiff_diaSrc_schema', 'deepDiff_diaSrc_schema_filename', 'deepDiff_diaSrc_schema_len', 'deepDiff_diaSrc_schema_md', 'deepDiff_differenceExp', 'deepDiff_differenceExp_bbox', 'deepDiff_differenceExp_detector', 'deepDiff_differenceExp_filename', 'deepDiff_differenceExp_filter', 'deepDiff_differenceExp_header_wcs', 'deepDiff_differenceExp_md', 'deepDiff_differenceExp_photoCalib', 'deepDiff_differenceExp_sub', 'deepDiff_differenceExp_visitInfo', 'deepDiff_differenceExp_wcs', 'deepDiff_kernelSrc', 'deepDiff_kernelSrc_filename', 'deepDiff_kernelSrc_len', 'deepDiff_kernelSrc_md', 'deepDiff_kernelSrc_schema', 'deepDiff_matchedExp', 'deepDiff_matchedExp_bbox', 'deepDiff_matchedExp_detector', 'deepDiff_matchedExp_filename', 'deepDiff_matchedExp_filter', 'deepDiff_matchedExp_header_wcs', 'deepDiff_matchedExp_md', 'deepDiff_matchedExp_photoCalib', 'deepDiff_matchedExp_sub', 'deepDiff_matchedExp_visitInfo', 'deepDiff_matchedExp_wcs', 'deepDiff_metadata', 'deepDiff_metadata_filename', 'deepDiff_warpedExp', 'deepDiff_warpedExp_bbox', 'deepDiff_warpedExp_detector', 'deepDiff_warpedExp_filename', 'deepDiff_warpedExp_filter', 'deepDiff_warpedExp_header_wcs', 'deepDiff_warpedExp_md', 'deepDiff_warpedExp_photoCalib', 'deepDiff_warpedExp_sub', 'deepDiff_warpedExp_visitInfo', 'deepDiff_warpedExp_wcs']
We can also use queryMetadata
to see what visits are available.
butler.queryMetadata('deepDiff_differenceExp', ('visit'), dataId={'filter': 'g'})[:10]
[410790, 410791, 410792, 410793, 410794, 410795, 410796, 410797, 410798, 410799]
Now we are ready to load some data.
visit_num = 410929
ccd_num = 9
First let's load the calexp and get the maskedImage
.
calexp = butler.get('calexp', {'visit': visit_num, 'ccdnum': ccd_num, 'filter': 'g'})
Two ways to get the maskedImage
out of the Exposure
object.
calexp_im = calexp.getMaskedImage()
# Alternate Way
calexp_im = calexp.maskedImage
Get the source catalog (src
)
calexp_src_cat = butler.get('src', {'visit': visit_num, 'ccdnum': ccd_num, 'filter': 'g'})
Get the difference image itself
diffexp = butler.get('deepDiff_differenceExp', {'visit': visit_num, 'ccdnum': ccd_num, 'filter':'g'})
diffexp_src_cat = butler.get('deepDiff_diaSrc', {'visit': visit_num, 'ccdnum': ccd_num, 'filter': 'g'})
diffexp_im = diffexp.getMaskedImage()
Since we want to get light curves in the end we'll need to learn how to get the time of the visits. Let's use getInfo
found in the Exposure
object.
exp_visit_info = calexp.getInfo().getVisitInfo()
visit_date = exp_visit_info.getDate()
print(visit_date)
DateTime("2015-02-17T03:51:13.595192500", TAI)
visit_date_python = exp_visit_info.getDate().toPython()
print(visit_date_python)
2015-02-17 03:51:13.595192
visit_date_astropy = Time(visit_date_python)
print(visit_date_astropy.mjd)
57070.160574018424
calexp
and diffExp
side-by-side¶fig = plt.figure()#figsize=(16, 14))
display = []
fig.add_subplot(1,2,1)
display.append(afw_display.Display(frame=fig))
display[0].scale('linear', 'zscale')
display[0].mtv(calexp_im)
fig.add_subplot(1,2,2)
display.append(afw_display.Display(frame=fig))
display[1].scale('linear', 'zscale')
#display[1].setMaskTransparency(10)
display[1].mtv(diffexp_im)
plt.tight_layout()
Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …
First, we get the mask plane from the difference image. Then we ask what items it contains, and what colors are used to display them.
mask = diffexp_im.mask
mask.getMaskPlaneDict().items()
dict_items([('BAD', 0), ('CLIPPED', 10), ('CR', 3), ('DETECTED', 5), ('DETECTED_NEGATIVE', 6), ('EDGE', 4), ('INEXACT_PSF', 11), ('INTRP', 2), ('NOT_DEBLENDED', 12), ('NO_DATA', 8), ('REJECTED', 13), ('SAT', 1), ('SENSOR_EDGE', 14), ('SUSPECT', 7), ('UNMASKEDNAN', 9)])
mask = diffexp.getMask()
for mask_name, mask_bit in mask.getMaskPlaneDict().items():
print('{:20}: {}'.format(mask_name, display[1].getMaskPlaneColor(mask_name)))
BAD : red CLIPPED : None CR : magenta DETECTED : blue DETECTED_NEGATIVE : cyan EDGE : yellow INEXACT_PSF : None INTRP : green NOT_DEBLENDED : None NO_DATA : orange REJECTED : None SAT : green SENSOR_EDGE : None SUSPECT : yellow UNMASKEDNAN : None
You can also mouse over the difference image mask and matplotlib will display the mask plane bit.
Using an asteroid finding algorithm called KBMOD we found Kuiper Belt Objects in the 2015 HiTS data that we are using here. In this folder, we've provided astrometry for the two objects we found in this field. Let's use 2015 DQ249
as an example and build its light curve.
The astrometry for these objects is provided in text files easily loaded in a pandas.DataFrame
.
hits_object_df = pd.read_csv('NotebookData/hits_kbmod_2015_DQ249_coords.dat', delimiter=' ')
hits_object_df.head()
visit | year | month | day | ra_hour | ra_min | ra_sec | dec_deg | dec_min | dec_sec | |
---|---|---|---|---|---|---|---|---|---|---|
0 | 410929 | 2015 | 2 | 17.15939 | 10 | 19 | 38.442 | -5 | 56 | 40.41 |
1 | 410985 | 2015 | 2 | 17.24529 | 10 | 19 | 37.616 | -5 | 56 | 33.08 |
2 | 411035 | 2015 | 2 | 17.31416 | 10 | 19 | 36.969 | -5 | 56 | 27.42 |
3 | 411069 | 2015 | 2 | 17.36113 | 10 | 19 | 36.501 | -5 | 56 | 23.76 |
4 | 411269 | 2015 | 2 | 18.08845 | 10 | 19 | 29.672 | -5 | 55 | 23.12 |
Let's use the information for visit number 410985 and find the asteroid in our difference image.
hits_visit = hits_object_df.loc[1]
Load the visit data. From our KBMOD search, I already know it is in CCD #9 so we can start with the data id already complete.
data_id = {'visit': hits_visit['visit'], 'ccdnum': 9, 'filter':'g'}
calexp = butler.get('calexp', data_id)
calexp_src_cat = butler.get('src', data_id)
diffexp = butler.get('deepDiff_differenceExp', data_id)
diffexp_src_cat = butler.get('deepDiff_diaSrc', data_id)
diffexp_src_df = diffexp_src_cat.asAstropy().to_pandas()
Brief Aside: There are lots of columns in the source catalog including quality flags. Here we show some potentially useful quality flags.
[x for x in diffexp_src_df.columns if x.startswith('base_PixelFlags')]
['base_PixelFlags_flag', 'base_PixelFlags_flag_offimage', 'base_PixelFlags_flag_edge', 'base_PixelFlags_flag_interpolated', 'base_PixelFlags_flag_saturated', 'base_PixelFlags_flag_cr', 'base_PixelFlags_flag_bad', 'base_PixelFlags_flag_suspect', 'base_PixelFlags_flag_interpolatedCenter', 'base_PixelFlags_flag_saturatedCenter', 'base_PixelFlags_flag_crCenter', 'base_PixelFlags_flag_suspectCenter']
To find the object we will load the image's WCS and use it to convert the ra, dec to pixel location.
wcs = diffexp.getWcs()
wcs
FITS standard SkyWcs: Sky Origin: (154.787109, -5.951107) Pixel Origin: (1171.94, 1762.08) Pixel Scale: 0.262966 arcsec/pixel
Use astropy SkyCoords to translate ra, dec into radians.
obj_pos = SkyCoord('%i %i %f %i %i %f' % (hits_visit['ra_hour'],
hits_visit['ra_min'],
hits_visit['ra_sec'],
hits_visit['dec_deg'],
hits_visit['dec_min'],
hits_visit['dec_sec']),
unit=(u.hourangle, u.degree))
obj_pos.ra.deg, obj_pos.dec.deg
(154.9067333333333, -5.942522222222222)
Use lsst.geom
package to create SpherePoint
that describes a position on the sky.
lsst.geom
also has units that we can provide.
obj_pos_lsst = lsst.geom.SpherePoint(obj_pos.ra.deg, obj_pos.dec.deg, lsst.geom.degrees)
Finally, we can use skyToPixel
to convert our ra, dec coordinates to pixel coordinates in the image.
x_pix, y_pix = wcs.skyToPixel(obj_pos_lsst)
We can also double check that our object is on this CCD by using getDimensions
and comparing to the pixel location the WCS gives us.
print(x_pix, y_pix)
print(diffexp.getDimensions())
1053.1148631559506 3388.9460528498885 (2048, 4096)
Factory
method to create cutouts¶One way to create a postage stamp is using the Factory
method. To use this we need to create a bounding box with lsst.geom
.
The origin
argument in the call to Factory
specifies that image pixel origin for our bounding box will be local to the cutout. (For more info on afwImage.LOCAL
vs afwImage.PARENT
see here.)
If deep
is set then the cutout will copy the data rather than using a reference.
x_half_width = 40
y_half_width = 40
# Define bounding box for cutout
bbox = lsst.geom.Box2I()
bbox.include(lsst.geom.Point2I(x_pix - x_half_width, y_pix - y_half_width))
bbox.include(lsst.geom.Point2I(x_pix + x_half_width, y_pix + y_half_width))
# Generate cutouts with Factory
calexp_cutout = calexp.Factory(calexp, bbox, origin=afwImage.LOCAL, deep=False)
diffexp_cutout = diffexp.Factory(diffexp, bbox, origin=afwImage.LOCAL, deep=False)
There is a new and easier way to get a cutout!
calexp_cutout = calexp.getCutout(obj_pos_lsst, size=lsst.geom.Extent2I(80, 80))
diffexp_cutout = diffexp.getCutout(obj_pos_lsst, size=lsst.geom.Extent2I(80, 80))
fig = plt.figure(figsize=(10, 5))
stamp_display = []
fig.add_subplot(1,2,1)
stamp_display.append(afw_display.Display(frame=fig))
stamp_display[0].scale('linear', 'zscale')
stamp_display[0].mtv(calexp_cutout.maskedImage)
#stamp_display[0].dot('o', x_pix, y_pix, size=4)
for src in calexp_src_cat:
stamp_display[0].dot('o', src.getX(), src.getY(), ctype='cyan', size=4)
plt.title('Calexp Image and Source Catalog')
fig.add_subplot(1,2,2)
stamp_display.append(afw_display.Display(frame=fig))
stamp_display[1].scale('linear', 'zscale')
stamp_display[1].mtv(diffexp_cutout.maskedImage)
#stamp_display[1].dot('o', x_pix, y_pix, size=4)
for src in diffexp_src_cat:
stamp_display[1].dot('o', src.getX(), src.getY(), ctype='cyan', size=4)
plt.title('Diffexp Image and Source Catalog')
plt.tight_layout()
Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …
In this folder there is astrometry for another asteroid that was found in the same field on the first night of HiTS observations. Load hits_kbmod_2014_XW40_coords.dat
into a dataframe and try to recreate the plot above for one of the visits.
Hint: If you run into an error trying to create a cutout it may be because while the asteroid is in the same field it may not fall on the same CCD. To find the correct CCD use the wcs.skyToPixel
function and this map of the DECam focal plane.
We are going to use Astropy and the match_to_catalog_sky
method to match our asteroid to the closest source in difference exposure source catalog. We already loaded our asteroid ra, dec into an Astropy SkyCoord
object above and called it obj_pos
. A SkyCoord
object can hold more than one set of coordinates though. So, we will load the coordinates from the source catalog into a SkyCoord
object called visit_coords
.
visit_coords = SkyCoord(diffexp_src_cat['coord_ra']*u.rad, diffexp_src_cat['coord_dec']*u.rad)
# obj_pos is only one row: The location of our asteroid in this visit.
obj_pos
<SkyCoord (ICRS): (ra, dec) in deg (154.90673333, -5.94252222)>
# visit_coords is many rows: The location of all detected sources in the source catalog
visit_coords[:10]
<SkyCoord (ICRS): (ra, dec) in deg [(154.65984136, -5.87743136), (154.65947515, -6.00393691), (154.66276295, -5.9437437 ), (154.66326343, -5.94459549), (154.66322514, -5.99797342), (154.66715343, -5.87290338), (154.66768744, -6.00104431), (154.66791396, -6.00582827), (154.66821111, -5.99871225), (154.66894971, -5.99688253)]>
Now we want to find the closest match to obj_pos
in visit_coords
so we use matchToCatalogSky
with visit_coords
as the argument. What we get back are the index in visit_coords
of the closest match and the 2-dimensional separation on the sky to that match. (If we had distance information we could use this to also get a closest 3-dimensional match.)
idx, sep2d, sep3d = obj_pos.match_to_catalog_sky(visit_coords)
print(idx, sep2d.arcsec)
84 [0.39408103]
We can use the source catalog directly to get instrumental fluxes and errors.
obj_instFlux = diffexp_src_cat.getPsfInstFlux()[idx]
print(obj_instFlux)
355.4253564211041
obj_instFlux_err = diffexp_src_cat.getPsfInstFluxErr()[idx]
print(obj_instFlux_err)
56.061742251235195
But what if we want calibrated fluxes and magnitudes along with the errors?
Use photoCalib
product.
deepDiff_photoCalib = diffexp.getPhotoCalib()
obj_g_flux = deepDiff_photoCalib.instFluxToNanojansky(obj_instFlux, obj_instFlux_err)
print(obj_g_flux)
# Access flux and error separately
print(obj_g_flux.value, obj_g_flux.error)
value=1280.067439057850379, error=201.9069922125679568 1280.0674390578504 201.90699221256796
obj_g_mag = deepDiff_photoCalib.instFluxToMagnitude(obj_instFlux, obj_instFlux_err)
print(obj_g_mag)
# Access flux and error separately
print(obj_g_mag.value, obj_g_mag.error)
value=23.63191787346008965, error=0.1712548298239273958 23.63191787346009 0.1712548298239274
def return_obj_skycoord(visit_data):
obj_pos = SkyCoord('%i %i %f %i %i %f' % (visit_data['ra_hour'],
visit_data['ra_min'],
visit_data['ra_sec'],
visit_data['dec_deg'],
visit_data['dec_min'],
visit_data['dec_sec']),
unit=(u.hourangle, u.degree))
return obj_pos
To do this we need to load in the differenceExp
and get the time of the observation. Then we get the source catalog and find the closest match to the known position of our asteroid coordinates.
To make sure we only keep good matches in our light curve we set a threshold of 1 arcsec on our matches. If the closest detected object is more than 1 arcsec we will move on to the next visit without a flux measurement.
If we do have a good match then we will use the photoCalib
to get a calibrated flux measurement for our light curve.
visit_time = []
visit_flux = []
visit_flux_err = []
visit_mag = []
for obs_idx in range(len(hits_object_df)):
# Load data
hits_visit = hits_object_df.iloc[obs_idx]
data_id = {'visit': hits_visit['visit'], 'ccdnum': 9, 'filter':'g'}
diffexp = butler.get('deepDiff_differenceExp', data_id)
diffexp_src_cat = butler.get('deepDiff_diaSrc', data_id)
exp_visit_info = diffexp.getInfo().getVisitInfo()
# Get Times
visit_date_python = exp_visit_info.getDate().toPython()
visit_date_astropy = Time(visit_date_python)
# Match to Difference Image Source Catalog
obj_pos = return_obj_skycoord(hits_visit)
visit_coords = SkyCoord(diffexp_src_cat['coord_ra']*u.rad,
diffexp_src_cat['coord_dec']*u.rad)
match_idx, match_sep2d, _ = obj_pos.match_to_catalog_sky(visit_coords)
# Only keep matches with 1 arcsecond. Otherwise skip this visit.
if match_sep2d.arcsec > 1.0:
print('No close matches for visit %i. Distance to closest match: %.2f arcsec' % (hits_visit['visit'], match_sep2d.arcsec))
continue
else:
print('Match within %.2f arcsec for visit %i' % (match_sep2d.arcsec, hits_visit['visit']))
# Load Flux for matched object
visit_time.append(visit_date_astropy.mjd)
inst_flux = diffexp_src_cat.getPsfInstFlux()[match_idx]
inst_flux_err = diffexp_src_cat.getPsfInstFluxErr()[match_idx]
deepDiff_photoCalib = diffexp.getPhotoCalib()
obj_flux = deepDiff_photoCalib.instFluxToNanojansky(inst_flux, inst_flux_err)
visit_flux.append(obj_flux.value)
visit_flux_err.append(obj_flux.error)
# For Exercise 2
# Load Magnitude for matched object or interesting columns from source catalog or visit information here.
visit_time = np.array(visit_time)
No close matches for visit 410929. Distance to closest match: 16.27 arcsec Match within 0.39 arcsec for visit 410985 No close matches for visit 411035. Distance to closest match: 33.96 arcsec No close matches for visit 411069. Distance to closest match: 35.25 arcsec No close matches for visit 411269. Distance to closest match: 85.68 arcsec Match within 0.53 arcsec for visit 411319 Match within 0.33 arcsec for visit 411369 Match within 0.64 arcsec for visit 411420 No close matches for visit 411470. Distance to closest match: 36.11 arcsec No close matches for visit 411671. Distance to closest match: 50.35 arcsec Match within 0.48 arcsec for visit 411772 Match within 0.61 arcsec for visit 411822 No close matches for visit 411872. Distance to closest match: 34.88 arcsec
fig = plt.figure()
plt.errorbar(visit_time - visit_time[0], visit_flux, yerr=visit_flux_err, marker='o', lw=1, elinewidth=2)
plt.xlabel('Time from First Observation (Days)')
plt.ylabel('Flux (nanojansky)')
plt.title('2015 DQ249 Light Curve in HiTS')
Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …
Text(0.5, 1.0, '2015 DQ249 Light Curve in HiTS')
And there is our light curve!
It seems that a few of our measurements did not have matching sources in the source catalog. KBMOD is designed to find objects that are not likely to be above the standard 5-sigma detection threshold in a single measurement so this is not surprising. An advanced exercise would be to rerun source detection with a lower threshold to try and get measurements for those visits and fill in the light curve!
Now that we have the basic infrastructure in place go back and make some other plots over time. An easy first step would be to make this plot with the magnitude of the source instead of the flux. Other ideas are to plot values from other columns in the source catalog or plot properties of the exposures from the exp_visit_info
object. Or remake this plot with the astrometry from 2014_XW40
.