Michael Wood-Vasey Last Verified to Run: 2019-07-17
After completing this Notebook, the user will be able to
See the Truth GCR Variables for a basic overview of the Truth information and how to access it for variables. https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/truth_gcr_variables.ipynb
# 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.)
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.')
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}');
Search for diaObjects that match input simulated SNe.
i = 0
sn = truth_sne.iloc[0]
print(sn)
Oops, we need to fix up the dtype
s 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)
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