DIA Analysis: Supernovae from the Run 1.2p Test

Michael Wood-Vasey Last Verified to Run: 2019-07-06

After completing this Notebook, the user will be able to

  1. Get a list of simulated SN from the Run 1.2p truth catalog
  2. Select SNe within the test patch in the Run 1.2p DIA Test
  3. Plot lightcurves from the DIA analysis
  4. Plot lightcurves from the truth catalog
  5. 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. https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/truth_gcr_variables.ipynb

In [1]:
# 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
In [2]:
# Requires issues/337 to fix obshistid query.
sys.path.insert(0, '/global/homes/w/wmwv/local/lsst/gcr-catalogs')
In [3]:
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
In [4]:
import GCRCatalogs
In [5]:
%matplotlib inline

import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
In [6]:
repo = '/global/cscratch1/sd/rearmstr/new_templates/diffim_template'
In [7]:
butler = Butler(repo)
In [8]:
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')
/opt/lsst/software/stack/python/miniconda3-4.5.12/envs/lsst-scipipe-1172c30/lib/python3.7/site-packages/GCR/base.py:373: UserWarning: Native quantity `psFluxMean_{band}` does not exist (required by `magMeanErr_i`, `magMeanErr_z`, `magMeanErr_r`, `magMeanErr_y`, `magMeanErr_u`, `magMeanErr_g`)
  warnings.warn(msg)

(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

In [9]:
columns = ['ra', 'dec', 'redshift', 'uniqueId', 'galaxy_id', 'sn']
truth_all_sne = pd.DataFrame(truth_cat.get_quantities(columns, filters=[f'sn == 1']))
/opt/lsst/software/stack/python/miniconda3-4.5.12/envs/lsst-scipipe-1172c30/lib/python3.7/site-packages/h5py/_hl/dataset.py:313: H5pyDeprecationWarning: dataset.value has been deprecated. Use dataset[()] instead.
  "Use dataset[()] instead.", H5pyDeprecationWarning)

You'll get a dataset.value deprecation warning. Don't worry about this. The Data Access Team will fix this someday.

In [10]:
print(f'We found {len(truth_all_sne)} simulated SNe.')
We found 76689 simulated SNe.

Select SNe from the truth catalog in the DIA test region

In [11]:
tract = 4849
patch = (6, 6)
In [12]:
skymap = butler.get('deepCoadd_skyMap')
tract_info = skymap[tract]
In [13]:
foo = tract_info.getPatchInfo(patch)
In [14]:
bar = foo.getOuterSkyPolygon(tract_info.getWcs())
In [15]:
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])
In [16]:
print(corners)
[[ 54.72776196 -29.78294629]
 [ 52.93572545 -29.78294586]
 [ 52.949112   -28.22767319]
 [ 54.71437636 -28.2276736 ]]
In [17]:
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)
In [18]:
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
In [19]:
columns = ['ra', 'dec', 'redshift', 'uniqueId', 'galaxy_id', 'sn']
truth_sne = pd.DataFrame(truth_cat.get_quantities(columns, filters=all_cuts))
/opt/lsst/software/stack/python/miniconda3-4.5.12/envs/lsst-scipipe-1172c30/lib/python3.7/site-packages/h5py/_hl/dataset.py:313: H5pyDeprecationWarning: dataset.value has been deprecated. Use dataset[()] instead.
  "Use dataset[()] instead.", H5pyDeprecationWarning)
In [20]:
print(f'We found {len(truth_sne)} simulated SNe.')
We found 8508 simulated SNe.
In [21]:
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.

In [22]:
i = 0
sn = truth_sne.iloc[0]
In [23]:
print(sn)
ra           5.300728e+01
uniqueId     3.082506e+14
galaxy_id    3.010259e+11
redshift     2.376888e-01
dec         -2.838755e+01
sn           1.000000e+00
Name: 0, dtype: float64

Oops, we need to fix up the dtypes somewhere. Those shouldn't all be floats.

Get the lightcurve

In [24]:
columns = ['obshistid', 'mjd', 'mag', 'filter']

sn_lc = pd.DataFrame(truth_lc.get_quantities(columns,
                                             native_filters=[f'uniqueId == {sn["uniqueId"]}']))
In [25]:
sn_lc.rename(columns={'filter': 'filter_code'}, inplace=True)
In [26]:
# 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')
In [27]:
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()
In [28]:
plt.figure(figsize=(12, 8))
plot_lightcurve(sn_lc, plot='mag', linestyle='-', marker=None)

Match to DIAObject Catalog

In [29]:
ra, dec = sn['ra'], sn['dec']
print(ra, dec)
53.00727812651823 -28.387552632959427
In [30]:
# 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
In [31]:
diaObjectId = diaObject_cat['diaObjectId'][idx]

How did we do with the lightcurve?

In [32]:
sn_lc
Out[32]:
mjd mag obshistid filter_code filter
0 60554.288025 26.432639 657401 0 u
1 60555.258471 25.581090 658296 0 u
2 60556.250832 25.019064 659233 0 u
3 60557.246086 24.602021 660201 0 u
4 60558.249432 24.266763 660882 0 u
5 60559.245016 23.986753 661527 0 u
6 60567.219180 21.397669 664557 2 r
7 60567.228879 21.461969 664577 1 g
8 60567.234411 21.431454 664587 3 i
9 60567.244111 21.719140 664607 4 z
10 60567.256310 21.978863 664633 5 y
11 60578.271042 21.354127 671659 2 r
12 60578.280324 21.694316 671678 1 g
13 60578.285857 21.346081 671688 3 i
14 60578.295556 21.825918 671708 4 z
15 60578.307755 22.263662 671734 5 y
16 60581.174302 21.450943 674105 2 r
17 60581.189152 21.472849 674134 3 i
18 60581.198851 21.983958 674154 4 z
19 60581.211050 22.407560 674180 5 y
20 60582.204745 24.257999 675069 0 u
21 60583.180252 24.374555 675910 0 u
22 60584.176755 24.494031 676739 0 u
23 60585.171706 24.608092 677595 0 u
24 60586.172558 24.727928 678467 0 u
25 60587.166237 24.855251 679317 0 u
26 60588.163054 24.992185 680167 0 u
27 60593.164657 22.093300 684469 2 r
28 60593.173940 23.101931 684488 1 g
29 60593.179472 22.013917 684498 3 i
... ... ... ... ... ...
43 60608.151733 22.571821 696296 4 z
44 60608.163933 22.559708 696322 5 y
45 60611.119189 23.094521 698622 2 r
46 60611.128471 24.440580 698641 1 g
47 60611.134004 22.555700 698651 3 i
48 60611.143703 22.708548 698671 4 z
49 60611.155902 22.649910 698697 5 y
50 60613.097971 26.392147 700187 0 u
51 60614.092032 26.435882 700968 0 u
52 60615.087792 26.482756 701754 0 u
53 60616.088773 26.532553 702553 0 u
54 60617.085683 26.583687 703320 0 u
55 60620.131032 26.722439 704240 0 u
56 60621.122783 23.508854 704746 2 r
57 60621.132065 24.736442 704765 1 g
58 60621.137598 22.996806 704775 3 i
59 60621.147297 23.206187 704795 4 z
60 60621.159496 23.115017 704821 5 y
61 60624.078591 23.592920 706923 2 r
62 60624.093405 23.095192 706952 3 i
63 60624.103104 23.344373 706972 4 z
64 60624.115303 23.258472 706998 5 y
65 60627.276161 23.667327 709647 2 r
66 60627.291393 23.191451 709677 3 i
67 60627.301092 23.484316 709697 4 z
68 60627.313291 23.379125 709723 5 y
69 60632.036052 23.747787 712881 2 r
70 60632.051295 23.312274 712911 3 i
71 60632.060994 23.667778 712931 4 z
72 60632.073193 23.491572 712957 5 y

73 rows × 5 columns

In [33]:
# 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')                    
In [34]:
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.

In [35]:
sn_lc.query('(60567 < mjd) & (mjd < 60568)')
Out[35]:
mjd mag obshistid filter_code filter
6 60567.219180 21.397669 664557 2 r
7 60567.228879 21.461969 664577 1 g
8 60567.234411 21.431454 664587 3 i
9 60567.244111 21.719140 664607 4 z
10 60567.256310 21.978863 664633 5 y
In [36]:
dia_lc.query('(60567 < mjd) & (mjd < 60568)')
Out[36]:
mjd psFlux mag_err psFluxErr visit mag filter
5 60567.219782 51954.298318 0.009033 430.379216 664557 19.610946 r
6 60567.235013 47920.015598 0.010489 461.672183 664587 19.698708 i
7 60567.244713 49375.749296 0.012521 568.633628 664607 19.666216 z

The MJD for the same visits (recall visit == obshistid) are slightly different between the truth catalog and the DIA lightcurve:

In [37]:
(60567.219180 - 60567.219782) * 24 * 3600
Out[37]:
-52.0128000061959

Why the 52 second difference?

Let's match on obshistid == visit to get a matched set of dataframes that we can compare.

In [38]:
joint_lc = pd.merge(sn_lc, dia_lc, left_on='obshistid', right_on='visit', suffixes=('_sim', '_dia'))
In [39]:
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']
In [40]:
joint_lc
Out[40]:
mjd_sim mag_sim obshistid filter_code filter_sim mjd_dia psFlux mag_err psFluxErr visit mag_dia filter_dia mjd filter delta_mag
0 60567.219180 21.397669 664557 2 r 60567.219782 51954.298318 0.009033 430.379216 664557 19.610946 r 60567.219782 r -1.786722
1 60567.234411 21.431454 664587 3 i 60567.235013 47920.015598 0.010489 461.672183 664587 19.698708 i 60567.235013 i -1.732747
2 60567.244111 21.719140 664607 4 z 60567.244713 49375.749296 0.012521 568.633628 664607 19.666216 z 60567.244713 z -2.052925
3 60578.271042 21.354127 671659 2 r 60578.271644 47085.850480 0.008895 385.539974 671659 19.717774 r 60578.271644 r -1.636353
4 60578.280324 21.694316 671678 1 g 60578.280926 10983.583655 0.020358 205.931007 671678 21.298140 g 60578.280926 g -0.396176
5 60578.285857 21.346081 671688 3 i 60578.286459 53817.660622 0.010054 498.132297 671688 19.572688 i 60578.286459 i -1.773393
6 60578.295556 21.825918 671708 4 z 60578.296158 57613.595367 0.012146 644.364773 671708 19.498688 z 60578.296158 z -2.327230
7 60578.307755 22.263662 671734 5 y 60578.308357 47248.056561 0.020698 900.609948 671734 19.714040 y 60578.308357 y -2.549622
8 60581.174302 21.450943 674105 2 r 60581.174904 38596.129793 0.011141 395.892130 674105 19.933641 r 60581.174904 r -1.517302
9 60593.173940 23.101931 684488 1 g 60593.174542 2152.571205 0.074385 147.471570 684488 23.067606 g 60593.174542 g -0.034324
10 60593.201370 22.608499 684544 5 y 60593.201972 28137.340049 0.026534 687.511028 684544 20.276792 y 60593.201972 y -2.331706
11 60596.178013 22.270623 686469 2 r 60596.178615 10590.229808 0.041427 404.050880 686469 21.337737 r 60596.178615 r -0.932887
12 60596.187712 23.396907 686489 1 g 60596.188314 1853.087931 0.159190 271.697872 686489 23.230260 g 60596.188314 g -0.166647
13 60596.193244 22.072659 686499 3 i 60596.193846 21364.683865 0.025884 509.285207 686499 20.575759 i 60596.193846 i -1.496901
14 60596.202943 22.472107 686519 4 z 60596.203545 29860.670592 0.021702 596.795208 686519 20.212251 z 60596.203545 z -2.259856
15 60596.215142 22.573407 686545 5 y 60596.215744 26717.792021 0.033466 823.473242 686545 20.332999 y 60596.215744 y -2.240409
16 60605.134462 22.778590 693871 2 r 60605.135064 5489.243765 0.048667 246.036487 693871 22.051219 r 60605.135064 r -0.727371
17 60605.171175 22.502005 693946 5 y 60605.171777 23313.259952 0.033979 729.486718 693946 20.480992 y 60605.171777 y -2.021012
18 60608.127220 22.936279 696247 2 r 60608.127822 3967.161602 0.059354 216.864375 696247 22.403800 r 60608.127822 r -0.532479
19 60608.136502 24.292038 696266 1 g 60608.137104 755.784122 0.173746 120.944561 696266 24.204006 g 60608.137104 g -0.088033
20 60608.142034 22.400639 696276 3 i 60608.142636 9545.790759 0.035801 314.726180 696276 21.450470 i 60608.142636 i -0.950169
21 60608.151733 22.571821 696296 4 z 60608.152335 21319.168260 0.022675 445.118346 696296 20.578074 z 60608.152335 z -1.993746
22 60608.163933 22.559708 696322 5 y 60608.164535 19782.167896 0.034725 632.615017 696322 20.659315 y 60608.164535 y -1.900393
23 60611.119189 23.094521 698622 2 r 60611.119791 3932.489517 0.071238 258.012421 698622 22.413331 r 60611.119791 r -0.681190
24 60611.128471 24.440580 698641 1 g 60611.129073 674.906577 0.168484 104.731225 698641 24.326891 g 60611.129073 g -0.113689
25 60611.134004 22.555700 698651 3 i 60611.134606 8362.159167 0.044743 344.553900 698651 21.594204 i 60611.134606 i -0.961496
26 60611.155902 22.649910 698697 5 y 60611.156504 16851.248510 0.038859 603.064822 698697 20.833420 y 60611.156504 y -1.816491
27 60621.122783 23.508854 704746 2 r 60621.123385 3054.910661 0.072668 204.460426 704746 22.687504 r 60621.123385 r -0.821351
28 60621.159496 23.115017 704821 5 y 60621.160098 11238.585132 0.060528 626.510025 704821 21.273221 y 60621.160098 y -1.841796
29 60624.093405 23.095192 706952 3 i 60624.094007 5518.862250 0.068156 346.433388 706952 22.045376 i 60624.094007 i -1.049815
30 60624.103104 23.344373 706972 4 z 60624.103706 12068.193213 0.042003 466.842536 706972 21.195894 z 60624.103706 z -2.148479
31 60624.115303 23.258472 706998 5 y 60624.115905 10180.921902 0.071504 670.470176 706998 21.380532 y 60624.115905 y -1.877940
32 60627.291393 23.191451 709677 3 i 60627.291995 4111.582997 0.140924 533.663023 709677 22.364977 i 60627.291995 i -0.826474
33 60627.313291 23.379125 709723 5 y 60627.313893 8028.566688 0.120211 888.906650 709723 21.638405 y 60627.313893 y -1.740720
34 60632.036052 23.747787 712881 2 r 60632.036654 2517.807350 0.136071 315.546726 712881 22.897444 r 60632.036654 r -0.850344
35 60632.060994 23.667778 712931 4 z 60632.061596 9233.395339 0.071836 610.905922 712931 21.486596 z 60632.061596 z -2.181181
In [41]:
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

In [ ]: