DIA Analysis: Supernovae from the Run 1.2p Test#
Michael Wood-Vasey Last Verified to Run: 2019-07-17
After completing this Notebook, the user will be able to
Get a list of simulated SN from the Run 1.2p truth catalog
Select SNe within the test patch in the Run 1.2p DIA Test
Plot lightcurves from the DIA analysis
Plot lightcurves from the truth catalog
Compare the input and recovered lightcurves.
See the Truth GCR Variables for a basic overview of the Truth information and how to access it for variables. LSSTDESC/DC2-analysis
# Inject gcr-catalogs that supports DIA source into path.
import os
import math
import sys
import numpy as np
import pandas as pd
from astropy.coordinates import SkyCoord
import astropy.units as u
import lsst.afw.display as afwDisplay
import lsst.afw.geom as afwGeom
from lsst.daf.persistence import Butler
from lsst.geom import SpherePoint
import lsst.geom
import GCRCatalogs
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
repo = '/global/cscratch1/sd/rearmstr/new_templates/diffim_template'
butler = Butler(repo)
diaSrc = GCRCatalogs.load_catalog('dc2_dia_source_run1.2p_test')
diaObject = GCRCatalogs.load_catalog('dc2_dia_object_run1.2p_test')
truth_cat = GCRCatalogs.load_catalog('dc2_truth_run1.2_variable_summary')
truth_lc = GCRCatalogs.load_catalog('dc2_truth_run1.2_variable_lightcurve')
(We presently will get a warning from the catalog reader in the initalization above because there is no u-band in the subtractions.)
Select SNe from the truth catalog#
columns = ['ra', 'dec', 'redshift', 'uniqueId', 'galaxy_id', 'sn']
truth_all_sne = pd.DataFrame(truth_cat.get_quantities(columns, filters=[f'sn == 1']))
You’ll get a dataset.value deprecation warning. Don’t worry about this. The Data Access Team will fix this someday.
print(f'We found {len(truth_all_sne)} simulated SNe.')
Select SNe from the truth catalog in the DIA test region#
tract = 4849
patch = (6, 6)
skymap = butler.get('deepCoadd_skyMap')
tract_info = skymap[tract]
foo = tract_info.getPatchInfo(patch)
bar = foo.getOuterSkyPolygon(tract_info.getWcs())
tract_box = afwGeom.Box2D(tract_info.getBBox())
tract_pos_list = tract_box.getCorners()
wcs = tract_info.getWcs()
corners = wcs.pixelToSky(tract_pos_list)
corners = np.array([[c.getRa().asDegrees(), c.getDec().asDegrees()] for c in corners])
print(corners)
ra = corners[:, 0]
dec = corners[:, 1]
min_ra, max_ra = np.min(ra), np.max(ra)
min_dec, max_dec = np.min(dec), np.max(dec)
area_cut = [f'ra > {min_ra}', f'ra < {max_ra}', f'dec > {min_dec}', f'dec < {max_dec}']
sn_cut = ['sn == 1']
all_cuts = area_cut + sn_cut
columns = ['ra', 'dec', 'redshift', 'uniqueId', 'galaxy_id', 'sn']
truth_sne = pd.DataFrame(truth_cat.get_quantities(columns, filters=all_cuts))
print(f'We found {len(truth_sne)} simulated SNe.')
avg_dec = (min_dec + max_dec)/2
size = 8
dec_size = size
ra_size = dec_size * np.cos(np.deg2rad(avg_dec))
aspect_ratio = dec_size / ra_size
# fig = plt.figure(figsize=(size, size))
ax = plt.gca()
ax.set_aspect(aspect_ratio)
patch_region = Polygon(corners, color='red', fill=False)
ax.scatter(truth_sne['ra'], truth_sne['dec'])
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
ax.set_xlim(plt.xlim()[::-1])
ax.add_patch(patch_region);
ax.set_title(f'tract: {tract}, patch: {patch}');
A Lightcurve#
Search for diaObjects that match input simulated SNe.
i = 0
sn = truth_sne.iloc[0]
print(sn)
Oops, we need to fix up the dtypes somewhere. Those shouldn’t all be floats.
Get the lightcurve
columns = ['obshistid', 'mjd', 'mag', 'filter']
sn_lc = pd.DataFrame(truth_lc.get_quantities(columns,
native_filters=[f'uniqueId == {sn["uniqueId"]}']))
sn_lc.rename(columns={'filter': 'filter_code'}, inplace=True)
# Translate filter codes to filter names
filter_names = ['u', 'g', 'r', 'i', 'z', 'y']
sn_lc['filter'] = [filter_names[f] for f in sn_lc['filter_code']]
sn_lc = sn_lc.sort_values('mjd')
def plot_lightcurve(df, plot='mag', flux_col_names=None,
title=None, marker='o', linestyle='none',
colors=None, label_prefix='',
**kwargs):
"""Plot a lightcurve from a DataFrame.
"""
# At lexigraphical order, if not wavelength order.
# Assume fixed list of filters.
filter_order = ['u', 'g', 'r', 'i', 'z', 'y']
if colors is None:
colors = {'u': 'violet', 'g': 'indigo', 'r': 'blue', 'i': 'green', 'z': 'orange', 'y': 'red'}
if flux_col_names is not None:
flux_col, flux_err_col = flux_col_names
else:
if plot == 'flux':
flux_col = 'psFlux'
flux_err_col = 'psFluxErr'
else:
flux_col = 'mag'
flux_err_col = 'mag_err'
for filt in filter_order:
this_filter = df.query(f'filter == "{filt}"')
if this_filter.empty:
continue
# This if sequence is a little silly.
plot_kwargs = {'linestyle': linestyle, 'marker': marker, 'color': colors[filt],
'label': f'{label_prefix} {filt}'}
plot_kwargs.update(kwargs)
if flux_err_col in this_filter.columns:
plt.errorbar(this_filter['mjd'], this_filter[flux_col], this_filter[flux_err_col],
**plot_kwargs)
else:
if marker is None:
plt.plot(this_filter['mjd'], this_filter[flux_col], **plot_kwargs)
else:
plot_kwargs.pop('linestyle')
plt.scatter(this_filter['mjd'], this_filter[flux_col], **plot_kwargs)
plt.xlabel('MJD')
if plot == 'flux':
plt.ylabel('psFlux [nJy]')
else:
# Ensure that y-axis decreases as one goes up
# Because plot_lightcurve could be called several times on the same axis,
# simply inverting is not correct. We have to reverse a sorted list.
plt.ylim(sorted(plt.ylim(), reverse=True))
plt.ylabel('mag [AB]')
if title is not None:
plt.title(title)
plt.legend()
plt.figure(figsize=(12, 8))
plot_lightcurve(sn_lc, plot='mag', linestyle='-', marker=None)
Match to DIAObject Catalog#
ra, dec = sn['ra'], sn['dec']
print(ra, dec)
# Match on RA, Dec
sn_position = SkyCoord(sn['ra'], sn['dec'], unit='deg')
diaObject_cat = diaObject.get_quantities(['ra', 'dec', 'diaObjectId'])
diaObject_positions = SkyCoord(diaObject_cat['ra'], diaObject_cat['dec'], unit='deg')
idx, sep2d, _ = sn_position.match_to_catalog_sky(diaObject_positions)
print(f'Index: {idx} is {sep2d.to(u.arcsec)[0]:0.6f} away')
diaObjectId = diaObject_cat['diaObjectId'][idx]
How did we do with the lightcurve?
sn_lc
# We can't use a direct filters = match in the GCR wrapper for the diaSrc table.
# So we have to use a lambda function here to match the ID
dia_lc = pd.DataFrame(diaSrc.get_quantities(['visit', 'mjd', 'psFlux', 'psFluxErr', 'mag', 'mag_err', 'filter'],
filters=[(lambda x: x == diaObjectId, 'diaObjectId')]))
dia_lc = dia_lc.sort_values('mjd')
plt.figure(figsize=(12, 8))
plot_lightcurve(sn_lc, plot='mag', linestyle='-', marker=None, label_prefix='Sim')
plot_lightcurve(dia_lc, plot='mag', label_prefix='DIA')
sim_date_range = [min(sn_lc['mjd']), max(sn_lc['mjd'])]
sim_date_delta = sim_date_range[1] - sim_date_range[0]
buffer_fraction = 0.05
plot_date_range = sim_date_range
plot_date_range[0] -= sim_date_delta * buffer_fraction
plot_date_range[1] += sim_date_delta * buffer_fraction
plt.xlim(plot_date_range)
plt.ylim(25, 19);
The shape seems good. But the calibration is magnitudes off. It does look like it’s a constant magnitude offset.
sn_lc.query('(60567 < mjd) & (mjd < 60568)')
dia_lc.query('(60567 < mjd) & (mjd < 60568)')
The MJD for the same visits (recall visit == obshistid) are slightly different between the truth catalog and the DIA lightcurve:
(60567.219180 - 60567.219782) * 24 * 3600
Why the 52 second difference?
Let’s match on obshistid == visit to get a matched set of dataframes that we can compare.
joint_lc = pd.merge(sn_lc, dia_lc, left_on='obshistid', right_on='visit', suffixes=('_sim', '_dia'))
joint_lc['mjd'] = joint_lc['mjd_dia']
joint_lc['filter'] = joint_lc['filter_sim']
joint_lc['delta_mag'] = joint_lc['mag_dia'] - joint_lc['mag_sim']
joint_lc
plot_lightcurve(joint_lc, flux_col_names=('delta_mag', 'mag_err'))
Hmmm…. so not really a constant offset in magnitude. Nevertheless, I can only suspect there’s some calibration or flux interpretation error somewhere.
If one wants to look at the postage stamps of this SN, get started with the examples in dia_source_object_stamp.ipynb