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.')
We found 76689 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])
/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/ipykernel/__main__.py:1: FutureWarning: Call to deprecated function (or staticmethod) Box2D. (Replaced by lsst.geom.Box2D (will be removed before the release of v20.0)) if __name__ == '__main__':
print(corners)
[[ 54.72776196 -29.78294629] [ 52.93572545 -29.78294586] [ 52.949112 -28.22767319] [ 54.71437636 -28.2276736 ]]
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.')
We found 8508 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)
galaxy_id 3.010259e+11 dec -2.838755e+01 sn 1.000000e+00 redshift 2.376888e-01 ra 5.300728e+01 uniqueId 3.082506e+14 Name: 0, dtype: float64
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)
53.00727812651823 -28.387552632959427
# 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')
Index: 22 is 0.028836 arcsec away
diaObjectId = diaObject_cat['diaObjectId'][idx]
How did we do with the lightcurve?
sn_lc
filter_code | mjd | mag | obshistid | filter | |
---|---|---|---|---|---|
0 | 0 | 60554.288025 | 26.432639 | 657401 | u |
1 | 0 | 60555.258471 | 25.581090 | 658296 | u |
2 | 0 | 60556.250832 | 25.019064 | 659233 | u |
3 | 0 | 60557.246086 | 24.602021 | 660201 | u |
4 | 0 | 60558.249432 | 24.266763 | 660882 | u |
5 | 0 | 60559.245016 | 23.986753 | 661527 | u |
6 | 2 | 60567.219180 | 21.397669 | 664557 | r |
7 | 1 | 60567.228879 | 21.461969 | 664577 | g |
8 | 3 | 60567.234411 | 21.431454 | 664587 | i |
9 | 4 | 60567.244111 | 21.719140 | 664607 | z |
10 | 5 | 60567.256310 | 21.978863 | 664633 | y |
11 | 2 | 60578.271042 | 21.354127 | 671659 | r |
12 | 1 | 60578.280324 | 21.694316 | 671678 | g |
13 | 3 | 60578.285857 | 21.346081 | 671688 | i |
14 | 4 | 60578.295556 | 21.825918 | 671708 | z |
15 | 5 | 60578.307755 | 22.263662 | 671734 | y |
16 | 2 | 60581.174302 | 21.450943 | 674105 | r |
17 | 3 | 60581.189152 | 21.472849 | 674134 | i |
18 | 4 | 60581.198851 | 21.983958 | 674154 | z |
19 | 5 | 60581.211050 | 22.407560 | 674180 | y |
20 | 0 | 60582.204745 | 24.257999 | 675069 | u |
21 | 0 | 60583.180252 | 24.374555 | 675910 | u |
22 | 0 | 60584.176755 | 24.494031 | 676739 | u |
23 | 0 | 60585.171706 | 24.608092 | 677595 | u |
24 | 0 | 60586.172558 | 24.727928 | 678467 | u |
25 | 0 | 60587.166237 | 24.855251 | 679317 | u |
26 | 0 | 60588.163054 | 24.992185 | 680167 | u |
27 | 2 | 60593.164657 | 22.093300 | 684469 | r |
28 | 1 | 60593.173940 | 23.101931 | 684488 | g |
29 | 3 | 60593.179472 | 22.013917 | 684498 | i |
... | ... | ... | ... | ... | ... |
43 | 4 | 60608.151733 | 22.571821 | 696296 | z |
44 | 5 | 60608.163933 | 22.559708 | 696322 | y |
45 | 2 | 60611.119189 | 23.094521 | 698622 | r |
46 | 1 | 60611.128471 | 24.440580 | 698641 | g |
47 | 3 | 60611.134004 | 22.555700 | 698651 | i |
48 | 4 | 60611.143703 | 22.708548 | 698671 | z |
49 | 5 | 60611.155902 | 22.649910 | 698697 | y |
50 | 0 | 60613.097971 | 26.392147 | 700187 | u |
51 | 0 | 60614.092032 | 26.435882 | 700968 | u |
52 | 0 | 60615.087792 | 26.482756 | 701754 | u |
53 | 0 | 60616.088773 | 26.532553 | 702553 | u |
54 | 0 | 60617.085683 | 26.583687 | 703320 | u |
55 | 0 | 60620.131032 | 26.722439 | 704240 | u |
56 | 2 | 60621.122783 | 23.508854 | 704746 | r |
57 | 1 | 60621.132065 | 24.736442 | 704765 | g |
58 | 3 | 60621.137598 | 22.996806 | 704775 | i |
59 | 4 | 60621.147297 | 23.206187 | 704795 | z |
60 | 5 | 60621.159496 | 23.115017 | 704821 | y |
61 | 2 | 60624.078591 | 23.592920 | 706923 | r |
62 | 3 | 60624.093405 | 23.095192 | 706952 | i |
63 | 4 | 60624.103104 | 23.344373 | 706972 | z |
64 | 5 | 60624.115303 | 23.258472 | 706998 | y |
65 | 2 | 60627.276161 | 23.667327 | 709647 | r |
66 | 3 | 60627.291393 | 23.191451 | 709677 | i |
67 | 4 | 60627.301092 | 23.484316 | 709697 | z |
68 | 5 | 60627.313291 | 23.379125 | 709723 | y |
69 | 2 | 60632.036052 | 23.747787 | 712881 | r |
70 | 3 | 60632.051295 | 23.312274 | 712911 | i |
71 | 4 | 60632.060994 | 23.667778 | 712931 | z |
72 | 5 | 60632.073193 | 23.491572 | 712957 | y |
73 rows × 5 columns
# 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)')
filter_code | mjd | mag | obshistid | filter | |
---|---|---|---|---|---|
6 | 2 | 60567.219180 | 21.397669 | 664557 | r |
7 | 1 | 60567.228879 | 21.461969 | 664577 | g |
8 | 3 | 60567.234411 | 21.431454 | 664587 | i |
9 | 4 | 60567.244111 | 21.719140 | 664607 | z |
10 | 5 | 60567.256310 | 21.978863 | 664633 | y |
dia_lc.query('(60567 < mjd) & (mjd < 60568)')
mjd | visit | mag_err | mag | psFlux | psFluxErr | filter | |
---|---|---|---|---|---|---|---|
1 | 60567.219782 | 664557 | 0.009033 | 19.610946 | 51954.298318 | 430.379216 | r |
2 | 60567.235013 | 664587 | 0.010489 | 19.698708 | 47920.015598 | 461.672183 | i |
3 | 60567.244713 | 664607 | 0.012521 | 19.666216 | 49375.749296 | 568.633628 | z |
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
-52.0128000061959
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
filter_code | mjd_sim | mag_sim | obshistid | filter_sim | mjd_dia | visit | mag_err | mag_dia | psFlux | psFluxErr | filter_dia | mjd | filter | delta_mag | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2 | 60567.219180 | 21.397669 | 664557 | r | 60567.219782 | 664557 | 0.009033 | 19.610946 | 51954.298318 | 430.379216 | r | 60567.219782 | r | -1.786722 |
1 | 3 | 60567.234411 | 21.431454 | 664587 | i | 60567.235013 | 664587 | 0.010489 | 19.698708 | 47920.015598 | 461.672183 | i | 60567.235013 | i | -1.732747 |
2 | 4 | 60567.244111 | 21.719140 | 664607 | z | 60567.244713 | 664607 | 0.012521 | 19.666216 | 49375.749296 | 568.633628 | z | 60567.244713 | z | -2.052925 |
3 | 2 | 60578.271042 | 21.354127 | 671659 | r | 60578.271644 | 671659 | 0.008895 | 19.717774 | 47085.850480 | 385.539974 | r | 60578.271644 | r | -1.636353 |
4 | 1 | 60578.280324 | 21.694316 | 671678 | g | 60578.280926 | 671678 | 0.020358 | 21.298140 | 10983.583655 | 205.931007 | g | 60578.280926 | g | -0.396176 |
5 | 3 | 60578.285857 | 21.346081 | 671688 | i | 60578.286459 | 671688 | 0.010054 | 19.572688 | 53817.660622 | 498.132297 | i | 60578.286459 | i | -1.773393 |
6 | 4 | 60578.295556 | 21.825918 | 671708 | z | 60578.296158 | 671708 | 0.012146 | 19.498688 | 57613.595367 | 644.364773 | z | 60578.296158 | z | -2.327230 |
7 | 5 | 60578.307755 | 22.263662 | 671734 | y | 60578.308357 | 671734 | 0.020698 | 19.714040 | 47248.056561 | 900.609948 | y | 60578.308357 | y | -2.549622 |
8 | 2 | 60581.174302 | 21.450943 | 674105 | r | 60581.174904 | 674105 | 0.011141 | 19.933641 | 38596.129793 | 395.892130 | r | 60581.174904 | r | -1.517302 |
9 | 1 | 60593.173940 | 23.101931 | 684488 | g | 60593.174542 | 684488 | 0.074385 | 23.067606 | 2152.571205 | 147.471570 | g | 60593.174542 | g | -0.034324 |
10 | 5 | 60593.201370 | 22.608499 | 684544 | y | 60593.201972 | 684544 | 0.026534 | 20.276792 | 28137.340049 | 687.511028 | y | 60593.201972 | y | -2.331706 |
11 | 2 | 60596.178013 | 22.270623 | 686469 | r | 60596.178615 | 686469 | 0.041427 | 21.337737 | 10590.229808 | 404.050880 | r | 60596.178615 | r | -0.932887 |
12 | 1 | 60596.187712 | 23.396907 | 686489 | g | 60596.188314 | 686489 | 0.159190 | 23.230260 | 1853.087931 | 271.697872 | g | 60596.188314 | g | -0.166647 |
13 | 3 | 60596.193244 | 22.072659 | 686499 | i | 60596.193846 | 686499 | 0.025884 | 20.575759 | 21364.683865 | 509.285207 | i | 60596.193846 | i | -1.496901 |
14 | 4 | 60596.202943 | 22.472107 | 686519 | z | 60596.203545 | 686519 | 0.021702 | 20.212251 | 29860.670592 | 596.795208 | z | 60596.203545 | z | -2.259856 |
15 | 5 | 60596.215142 | 22.573407 | 686545 | y | 60596.215744 | 686545 | 0.033466 | 20.332999 | 26717.792021 | 823.473242 | y | 60596.215744 | y | -2.240409 |
16 | 2 | 60605.134462 | 22.778590 | 693871 | r | 60605.135064 | 693871 | 0.048667 | 22.051219 | 5489.243765 | 246.036487 | r | 60605.135064 | r | -0.727371 |
17 | 5 | 60605.171175 | 22.502005 | 693946 | y | 60605.171777 | 693946 | 0.033979 | 20.480992 | 23313.259952 | 729.486718 | y | 60605.171777 | y | -2.021012 |
18 | 2 | 60608.127220 | 22.936279 | 696247 | r | 60608.127822 | 696247 | 0.059354 | 22.403800 | 3967.161602 | 216.864375 | r | 60608.127822 | r | -0.532479 |
19 | 1 | 60608.136502 | 24.292038 | 696266 | g | 60608.137104 | 696266 | 0.173746 | 24.204006 | 755.784122 | 120.944561 | g | 60608.137104 | g | -0.088033 |
20 | 3 | 60608.142034 | 22.400639 | 696276 | i | 60608.142636 | 696276 | 0.035801 | 21.450470 | 9545.790759 | 314.726180 | i | 60608.142636 | i | -0.950169 |
21 | 4 | 60608.151733 | 22.571821 | 696296 | z | 60608.152335 | 696296 | 0.022675 | 20.578074 | 21319.168260 | 445.118346 | z | 60608.152335 | z | -1.993746 |
22 | 5 | 60608.163933 | 22.559708 | 696322 | y | 60608.164535 | 696322 | 0.034725 | 20.659315 | 19782.167896 | 632.615017 | y | 60608.164535 | y | -1.900393 |
23 | 2 | 60611.119189 | 23.094521 | 698622 | r | 60611.119791 | 698622 | 0.071238 | 22.413331 | 3932.489517 | 258.012421 | r | 60611.119791 | r | -0.681190 |
24 | 1 | 60611.128471 | 24.440580 | 698641 | g | 60611.129073 | 698641 | 0.168484 | 24.326891 | 674.906577 | 104.731225 | g | 60611.129073 | g | -0.113689 |
25 | 3 | 60611.134004 | 22.555700 | 698651 | i | 60611.134606 | 698651 | 0.044743 | 21.594204 | 8362.159167 | 344.553900 | i | 60611.134606 | i | -0.961496 |
26 | 5 | 60611.155902 | 22.649910 | 698697 | y | 60611.156504 | 698697 | 0.038859 | 20.833420 | 16851.248510 | 603.064822 | y | 60611.156504 | y | -1.816491 |
27 | 2 | 60621.122783 | 23.508854 | 704746 | r | 60621.123385 | 704746 | 0.072668 | 22.687504 | 3054.910661 | 204.460426 | r | 60621.123385 | r | -0.821351 |
28 | 5 | 60621.159496 | 23.115017 | 704821 | y | 60621.160098 | 704821 | 0.060528 | 21.273221 | 11238.585132 | 626.510025 | y | 60621.160098 | y | -1.841796 |
29 | 3 | 60624.093405 | 23.095192 | 706952 | i | 60624.094007 | 706952 | 0.068156 | 22.045376 | 5518.862250 | 346.433388 | i | 60624.094007 | i | -1.049815 |
30 | 4 | 60624.103104 | 23.344373 | 706972 | z | 60624.103706 | 706972 | 0.042003 | 21.195894 | 12068.193213 | 466.842536 | z | 60624.103706 | z | -2.148479 |
31 | 5 | 60624.115303 | 23.258472 | 706998 | y | 60624.115905 | 706998 | 0.071504 | 21.380532 | 10180.921902 | 670.470176 | y | 60624.115905 | y | -1.877940 |
32 | 3 | 60627.291393 | 23.191451 | 709677 | i | 60627.291995 | 709677 | 0.140924 | 22.364977 | 4111.582997 | 533.663023 | i | 60627.291995 | i | -0.826474 |
33 | 5 | 60627.313291 | 23.379125 | 709723 | y | 60627.313893 | 709723 | 0.120211 | 21.638405 | 8028.566688 | 888.906650 | y | 60627.313893 | y | -1.740720 |
34 | 2 | 60632.036052 | 23.747787 | 712881 | r | 60632.036654 | 712881 | 0.136071 | 22.897444 | 2517.807350 | 315.546726 | r | 60632.036654 | r | -0.850344 |
35 | 4 | 60632.060994 | 23.667778 | 712931 | z | 60632.061596 | 712931 | 0.071836 | 21.486596 | 9233.395339 | 610.905922 | z | 60632.061596 | z | -2.181181 |
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