Rubin DPO

TAP

Much of this code is taken from the DP0 notebooks.

import os
import requests
from astropy import units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
import pyvo

params = {'axes.labelsize': 28,
          'font.size': 24,
          'legend.fontsize': 18,
          'xtick.major.width': 3,
          'xtick.minor.width': 2,
          'xtick.major.size': 12,
          'xtick.minor.size': 6,
          'xtick.direction': 'in',
          'xtick.top': True,
          'lines.linewidth': 3,
          'axes.linewidth': 3,
          'axes.labelweight': 3,
          'axes.titleweight': 3,
          'ytick.major.width': 3,
          'ytick.minor.width': 2,
          'ytick.major.size': 12,
          'ytick.minor.size': 6,
          'ytick.direction': 'in',
          'ytick.right': True,
          'figure.figsize': [10, 8],
          'figure.facecolor': 'White'
          }
plt.rcParams.update(params)

def _get_auth(tap_url, access_token):
    """https://github.com/lsst-sqre/rubin-jupyter-lab/blob/master/rubin_jupyter_utils/lab/notebook/catalog.py"""
    # tap_url = _get_tap_url()
    s = requests.Session()
    s.headers["Authorization"] = "Bearer " + access_token
    auth = pyvo.auth.authsession.AuthSession()
    auth.credentials.set("lsst-token", s)
    auth.add_security_method_for_url(tap_url, "lsst-token")
    auth.add_security_method_for_url(tap_url + "/sync", "lsst-token")
    auth.add_security_method_for_url(tap_url + "/async", "lsst-token")
    auth.add_security_method_for_url(tap_url + "/tables", "lsst-token")
    return auth

tap_url = "https://data.lsst.cloud/api/tap"
# username = "x-oauth-basic"
auth = _get_auth(tap_url, os.getenv("DP0token"))
service = pyvo.dal.TAPService(tap_url, auth)

# get table descriptions
results = service.search("SELECT description, table_name FROM TAP_SCHEMA.tables")
resultsdf = results.to_table().to_pandas()

# --- CMD and CCD from cone search
query = "SELECT obj.objectId, obj.ra, obj.dec, obj.mag_g, obj.mag_r, "\
        "obj.mag_i, obj.mag_g_cModel, obj.mag_r_cModel, obj.mag_i_cModel, "\
        "obj.psFlux_g, obj.psFlux_r, obj.psFlux_i, obj.cModelFlux_g, "\
        "obj.cModelFlux_r, obj.cModelFlux_i, obj.tract, obj.patch, "\
        "obj.extendedness, obj.good, obj.clean, "\
        "truth.mag_r as truth_mag_r, truth.match_objectId, "\
        "truth.flux_g, truth.flux_r, truth.flux_i, truth.truth_type,  "\
        "truth.match_sep, truth.is_variable "\
        "FROM dp01_dc2_catalogs.object as obj "\
        "JOIN dp01_dc2_catalogs.truth_match as truth "\
        "ON truth.match_objectId = obj.objectId "\
        "WHERE CONTAINS(POINT('ICRS', obj.ra, obj.dec), "\
        "CIRCLE('ICRS', 62.0, -37.0, 0.10)) = 1 "\
        "AND truth.match_objectid >= 0 "\
        "AND truth.is_good_match = 1"
results = service.search(query)
df = results.to_table().to_pandas()
# separate stars and galaxies
star = np.where(results['truth_type'] == 2)
gx = np.where(results['truth_type'] == 1)

ccd_cmd(results, star, gx)
compare_fluxes(results, star, gx)

def ccd_cmd(results, star, gx):
    fig, ax = plt.subplots(1, 2, figsize=(15, 8))

    plt.sca(ax[0])  # set the first axis as current

    plt.plot(results['mag_g_cModel'][gx] - results['mag_i_cModel'][gx],
             results['mag_g_cModel'][gx], 'k.', alpha=0.2, label='galaxies')
    plt.plot(results['mag_g_cModel'][star] - results['mag_i_cModel'][star],
             results['mag_g_cModel'][star], 'ro', label='stars')
    plt.legend(loc='upper left')
    plt.xlabel(r'$(g-i)$')
    plt.ylabel(r'$g$')
    plt.xlim(-1.8, 4.3)
    plt.ylim(29.3, 16.7)
    plt.minorticks_on()

    plt.sca(ax[1])  # set the first axis as current
    plt.plot(results['mag_g_cModel'][gx] - results['mag_r_cModel'][gx],
             results['mag_r_cModel'][gx] - results['mag_i_cModel'][gx],
             'k.', alpha=0.1, label='galaxies')
    plt.plot(results['mag_g_cModel'][star] - results['mag_r_cModel'][star],
             results['mag_r_cModel'][star] - results['mag_i_cModel'][star],
             'ro', label='stars')
    plt.legend(loc='upper left')
    plt.xlabel(r'$(g-r)$')
    plt.ylabel(r'$(r-i)$')
    plt.xlim(-1.3, 2.3)
    plt.ylim(-1.3, 2.8)
    plt.minorticks_on()

    plt.tight_layout()
    plt.show(block=False)

def compare_fluxes(results, star, gx):
    plt.rcParams.update({'figure.figsize': (11, 10)})

    plt.plot(results['truth_mag_r'][gx],
             results['cModelFlux_r'][gx] / results['flux_r'][gx],
             'k.', alpha=0.2, label='galaxies')
    plt.plot(results['truth_mag_r'][star],
             results['cModelFlux_r'][star] / results['flux_r'][star],
             'ro', label='stars')
    plt.legend(loc='upper left')
    plt.xlabel(r'$r$ magnitude (truth)')
    plt.ylabel(r'$f_{\rm meas}/f_{\rm truth}$')
    plt.ylim(0.15, 2.15)
    plt.xlim(17.6, 27.8)
    plt.minorticks_on()
    plt.show(block=False)


# --- Cone search
coord = SkyCoord(ra=62.0*u.degree, dec=-37.0*u.degree, frame='icrs')
radius = 0.1 * u.deg
query = "SELECT ra, dec, mag_g, mag_i " \
        "mag_i, mag_g_cModel, mag_r_cModel, mag_i_cModel, " \
        "psFlux_g, psFlux_r, psFlux_i, " \
        "cModelFlux_g, cModelFlux_r, cModelFlux_i, " \
        "tract, patch, extendedness, good, clean " \
        "FROM dp01_dc2_catalogs.object " \
        "WHERE CONTAINS(POINT('ICRS', ra, dec),CIRCLE('ICRS', " \
        + str(coord.ra.value) + ", " + str(coord.dec.value) + ", " \
        + str(radius.value) + " )) = 1"

df = service.search(query).to_table().to_pandas()

# join with truth
query = "SELECT obj.objectId, obj.ra, obj.dec, obj.mag_g, obj.mag_r, " \
        " obj.mag_i, obj.mag_g_cModel, obj.mag_r_cModel, obj.mag_i_cModel," \
        "obj.psFlux_g, obj.psFlux_r, obj.psFlux_i, obj.cModelFlux_g, " \
        "obj.cModelFlux_r, obj.cModelFlux_i, obj.tract, obj.patch, " \
        "obj.extendedness, obj.good, obj.clean, " \
        "truth.mag_r as truth_mag_r, truth.match_objectId, "\
        "truth.flux_g, truth.flux_r, truth.flux_i, truth.truth_type, " \
        "truth.match_sep, truth.is_variable " \
        "FROM dp01_dc2_catalogs.object as obj " \
        "JOIN dp01_dc2_catalogs.truth_match as truth " \
        "ON truth.match_objectId = obj.objectId " \
        "WHERE CONTAINS(POINT('ICRS', obj.ra, obj.dec),CIRCLE('ICRS', " \
        + str(coord.ra.value) + ", " + str(coord.dec.value) + ", " \
        + str(radius.value) + " )) = 1 " \
        "AND truth.match_objectid >= 0 "\
        "AND truth.is_good_match = 1"
df = service.search(query).to_table().to_pandas()
# How many of each type in the dataset
# The 'truth_type' in the truth_match table is 1= galaxies, 2=stars, 3=SNe
n_stars = df[df["truth_type"] == 2].shape[0]
print(f'There are {n_stars} stars out of a total of {len(df)}')
print(f'There are {df[df["truth_type"] == 1].shape[0]} galaxies')
print(f'There are {df[df["truth_type"] == 3].shape[0]} SNe')