06_Comparing_Object_and_Truth_Tables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e9a2fea9",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<br>\n",
    "<b>Comparing Object and Truth Tables</b> <br>\n",
    "Last verified to run on 2022-05-31 with LSST Science Pipelines release w_2022_22<br>\n",
    "Contact author: Jeff Carlin <br>\n",
    "Target audience: All DP0 delegates. <br>\n",
    "Container Size: medium <br>\n",
    "Questions welcome at <a href=\"https://community.lsst.org/c/support/dp0\">community.lsst.org/c/support/dp0</a> <br>\n",
    "Find DP0 documentation and resources at <a href=\"https://dp0-1.lsst.io\">dp0-1.lsst.io</a> <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6222d63c",
   "metadata": {},
   "source": [
    "**Credit:** This tutorial was developed by Jeff Carlin in the context of the DP0. Please consider acknowledging Jeff Carlin in any publications or software releases that make use of this notebook's contents."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b181bce",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "\n",
    "In this short tutorial, users will learn how to extract a table containing data from both the DP0.1 `Object` and `Truth-Match` tables for each object in common between them. The query returns the table `JOIN` of these two catalogs, which enables comparison of the recovered (measured) properties (e.g., fluxes, positions, magnitudes, etc.) to the simulated values that were assigned to each object when creating the DC2 simulations.\n",
    "\n",
    "More information about the DC2 simulations that make up DP0.1 can be found in [the DC2 Data Release Note](https://ui.adsabs.harvard.edu/abs/2021arXiv210104855L/abstract).\n",
    "\n",
    "### Set Up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db0c9660",
   "metadata": {},
   "outputs": [],
   "source": [
    "# What version of the Stack are we using?\n",
    "! echo $IMAGE_DESCRIPTION\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a74653c",
   "metadata": {},
   "source": [
    "### 1. Import Common Python Libraries\n",
    "\n",
    "The [`matplotlib`](https://matplotlib.org/), [`numpy`](http://www.numpy.org/), [`pandas`](https://pandas.pydata.org/docs/), and [`astropy`](http://www.astropy.org/) libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis. We will use these packages below, including the `matplotlib.pyplot` plotting sublibrary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a75a320a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# allow for matplotlib to create inline plots in our notebook\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas                        # imports the pandas data analysis tools\n",
    "import matplotlib.pyplot as plt      # imports matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "762a6ccc-b007-476b-8ac6-f15a8c374693",
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter some unit warnings\n",
    "import warnings\n",
    "from astropy.units import UnitsWarning\n",
    "warnings.simplefilter(\"ignore\", category=UnitsWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "379c34ed",
   "metadata": {},
   "source": [
    "To access tables, we will use the TAP service in a similar manner to what we showed in the [Intro to DP0 notebook](https://github.com/rubin-dp0/tutorial-notebooks/blob/main/01_Intro_to_DP0_Notebooks.ipynb), and explored further in the [TAP tutorial notebook](https://github.com/rubin-dp0/tutorial-notebooks/blob/main/02_Intermediate_TAP_Query.ipynb). See those notebooks for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5803048d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up some options, and import a couple more tools we will need:\n",
    "pandas.set_option('display.max_rows', 200)\n",
    "\n",
    "from lsst.rsp import get_tap_service, retrieve_query\n",
    "\n",
    "service = get_tap_service()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6cc2670",
   "metadata": {},
   "source": [
    "### 2. Loading tables with TAP\n",
    "\n",
    "What tables are available?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc7d10ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT description,\\\n",
    "                          table_name FROM TAP_SCHEMA.tables\")\n",
    "results_tab = results.to_table()\n",
    "results_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca9b3dc9",
   "metadata": {},
   "source": [
    "For our analysis, let's choose the Object table, `dp01_dc2_catalogs.object`, and then we will compare the measurements from this table to the \"truth\" values from `dp01_dc2_catalogs.truth_match`.\n",
    "\n",
    "For later reference, let's print out the table schema (i.e., the list of columns) for each of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63c516fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Object table:\n",
    "\n",
    "results = service.search(\"SELECT column_name, datatype, description,\\\n",
    "                          unit from TAP_SCHEMA.columns\\\n",
    "                          WHERE table_name = 'dp01_dc2_catalogs.object'\")\n",
    "# Note that we use the .to_pandas() method here so that all rows will display.\n",
    "#   Astropy will truncate the table for display, whereas we set the maximum\n",
    "#   number of rows for pandas to display to 200 in a cell above.\n",
    "results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31d870ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Truth-match table\n",
    "\n",
    "results = service.search(\"SELECT column_name, datatype, description,\\\n",
    "                          unit from TAP_SCHEMA.columns\\\n",
    "                          WHERE table_name = 'dp01_dc2_catalogs.truth_match'\")\n",
    "results_tab = results.to_table()\n",
    "results_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "943a1e80",
   "metadata": {},
   "source": [
    "### 3. Extract a table of joined results from two tables, using a single query within ADQL\n",
    "\n",
    "For those who are not familiar with databases, the simplest way of accomplishing the task might seem to be querying the two tables (Object and Truth-Match) separately, then matching them afterwards. This would work, but it is not the best use of the resources that are available. As seen in the [Advanced TAP/ADQL Usage in the Portal Aspect](https://dp0-1.lsst.io/tutorials-examples/index-portal-advanced.html#examples-dp0-1-portal-advanced) tutorial, the table JOIN can be done directly with ADQL. Not only will this save you a few steps, it should also execute much faster.\n",
    "\n",
    "For this exploration, we will select a small region of sky around a random RA, Dec position. The following cell reads data centered on (RA, Dec) = (62.0, -37.0) degrees, within a radius of 0.1 degrees, from both the Object and Truth-Match table, then joins database entries where `match_objectId` from Truth-Match equals `objectId` from the Object table. Note that we are selecting only a subset of the columns seen in the schema above. You can add or remove columns as you wish.\n",
    "\n",
    "Note that for the Object table we select all objects within the cone-shaped region of interest. In the Truth-Match table, we restrict the results to objects satisfying \"match_objectId >= 0 AND is_good_match = 1\". According to the Truth-Match schema above, the `is_good_match` flag is \"True if this object--truth matching pair satisfies all matching criteria\" as laid out in the [DESC DC2 Release Note](https://ui.adsabs.harvard.edu/abs/2021arXiv210104855L/abstract). We'll use that to select \"good\" matches. In the column description for `match_objectId` from above, it says \"objectId of the matching object entry (-1 for unmatched truth entries).\" Thus the criterion \"match_objectId >= 0\" removes the unmatched entries, leaving us with only the truth-table entries that were detected and appear in the Object table.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "441468a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the query\n",
    "query = \"SELECT obj.objectId, obj.ra, obj.dec, obj.mag_g, obj.mag_r, \"\\\n",
    "        \"obj.mag_i, obj.mag_g_cModel, obj.mag_r_cModel, obj.mag_i_cModel, \"\\\n",
    "        \"obj.psFlux_g, obj.psFlux_r, obj.psFlux_i, obj.cModelFlux_g, \"\\\n",
    "        \"obj.cModelFlux_r, obj.cModelFlux_i, obj.tract, obj.patch, \"\\\n",
    "        \"obj.extendedness, obj.good, obj.clean, \"\\\n",
    "        \"truth.mag_r as truth_mag_r, truth.match_objectId, \"\\\n",
    "        \"truth.flux_g, truth.flux_r, truth.flux_i, truth.truth_type,  \"\\\n",
    "        \"truth.match_sep, truth.is_variable \"\\\n",
    "        \"FROM dp01_dc2_catalogs.object as obj \"\\\n",
    "        \"JOIN dp01_dc2_catalogs.truth_match as truth \"\\\n",
    "        \"ON truth.match_objectId = obj.objectId \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', obj.ra, obj.dec), \"\\\n",
    "        \"CIRCLE('ICRS', 62.0, -37.0, 0.10)) = 1 \"\\\n",
    "        \"AND truth.match_objectid >= 0 \"\\\n",
    "        \"AND truth.is_good_match = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "679e8229",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "truth_plus_meas = service.search(query)\n",
    "print('Query returned %s matched objects.' % len(truth_plus_meas))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9588f0b2",
   "metadata": {},
   "source": [
    "Notice the `%%time` [cell magic](https://ipython.readthedocs.io/en/stable/interactive/magics.html) we used. This was included to highlight that selecting more than 14000 objects from two tables and joining them takes less than 2 seconds (typically). Doing two separate queries, then joining them using `pandas`, typically takes more than a full minute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9775666",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the names of the fields in the returned table:\n",
    "\n",
    "truth_plus_meas.fieldnames"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ac08554",
   "metadata": {},
   "source": [
    "### 4. Compare table values by plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1806b66c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up some plotting defaults:\n",
    "\n",
    "params = {'axes.labelsize': 28,\n",
    "          'font.size': 24,\n",
    "          'legend.fontsize': 18,\n",
    "          'xtick.major.width': 3,\n",
    "          'xtick.minor.width': 2,\n",
    "          'xtick.major.size': 12,\n",
    "          'xtick.minor.size': 6,\n",
    "          'xtick.direction': 'in',\n",
    "          'xtick.top': True,\n",
    "          'lines.linewidth': 3,\n",
    "          'axes.linewidth': 3,\n",
    "          'axes.labelweight': 3,\n",
    "          'axes.titleweight': 3,\n",
    "          'ytick.major.width': 3,\n",
    "          'ytick.minor.width': 2,\n",
    "          'ytick.major.size': 12,\n",
    "          'ytick.minor.size': 6,\n",
    "          'ytick.direction': 'in',\n",
    "          'ytick.right': True,\n",
    "          'figure.figsize': [10, 8],\n",
    "          'figure.facecolor': 'White'\n",
    "          }\n",
    "\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b15743bb",
   "metadata": {},
   "source": [
    "#### Compare the measurements from the Object table to the \"true\" values for some objects.\n",
    "\n",
    "To do this, we will separate the \"stars\" and \"galaxies\" using the `truth_type` column from the Truth-Match table. Simulated stars have `truth_type = 2`, and galaxies, `truth_type = 1`.\n",
    "\n",
    "After separating stars and galaxies, we'll compare the recovered flux to the \"true\" value that was simulated for each object (as a ratio of the fluxes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "401ba83e",
   "metadata": {},
   "outputs": [],
   "source": [
    "star = np.where(truth_plus_meas['truth_type'] == 2)\n",
    "gx = np.where(truth_plus_meas['truth_type'] == 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd57ca6a",
   "metadata": {},
   "source": [
    "Just to confirm that things look like we expect, let's plot a color-magnitude (g vs. g-i) and color-color (r-i vs. g-r) diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "884763ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(15, 8))\n",
    "\n",
    "plt.sca(ax[0])  # set the first axis as current\n",
    "\n",
    "plt.plot(truth_plus_meas['mag_g_cModel'][gx] - truth_plus_meas['mag_i_cModel'][gx],\n",
    "         truth_plus_meas['mag_g_cModel'][gx], 'k.', alpha=0.2, label='galaxies')\n",
    "plt.plot(truth_plus_meas['mag_g_cModel'][star] - truth_plus_meas['mag_i_cModel'][star],\n",
    "         truth_plus_meas['mag_g_cModel'][star], 'ro', label='stars')\n",
    "plt.legend(loc='upper left')\n",
    "plt.xlabel(r'$(g-i)$')\n",
    "plt.ylabel(r'$g$')\n",
    "plt.xlim(-1.8, 4.3)\n",
    "plt.ylim(29.3, 16.7)\n",
    "plt.minorticks_on()\n",
    "\n",
    "plt.sca(ax[1])  # set the first axis as current\n",
    "plt.plot(truth_plus_meas['mag_g_cModel'][gx] - truth_plus_meas['mag_r_cModel'][gx],\n",
    "         truth_plus_meas['mag_r_cModel'][gx] - truth_plus_meas['mag_i_cModel'][gx],\n",
    "         'k.', alpha=0.1, label='galaxies')\n",
    "plt.plot(truth_plus_meas['mag_g_cModel'][star] - truth_plus_meas['mag_r_cModel'][star],\n",
    "         truth_plus_meas['mag_r_cModel'][star] - truth_plus_meas['mag_i_cModel'][star],\n",
    "         'ro', label='stars')\n",
    "plt.legend(loc='upper left')\n",
    "plt.xlabel(r'$(g-r)$')\n",
    "plt.ylabel(r'$(r-i)$')\n",
    "plt.xlim(-1.3, 2.3)\n",
    "plt.ylim(-1.3, 2.8)\n",
    "plt.minorticks_on()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "682c5d15",
   "metadata": {},
   "source": [
    "Looks pretty normal - the stellar locus in color-color space is right where one expects it to be, and the galaxies dominate at the faint end of the CMD. \n",
    "\n",
    "Now let's compare the fluxes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87318bec",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize': (11, 10)})\n",
    "\n",
    "plt.plot(truth_plus_meas['truth_mag_r'][gx],\n",
    "         truth_plus_meas['cModelFlux_r'][gx] / truth_plus_meas['flux_r'][gx],\n",
    "         'k.', alpha=0.2, label='galaxies')\n",
    "plt.plot(truth_plus_meas['truth_mag_r'][star],\n",
    "         truth_plus_meas['cModelFlux_r'][star] / truth_plus_meas['flux_r'][star],\n",
    "         'ro', label='stars')\n",
    "plt.legend(loc='upper left')\n",
    "plt.xlabel(r'$r$ magnitude (truth)')\n",
    "plt.ylabel(r'$f_{\\rm meas}/f_{\\rm truth}$')\n",
    "plt.ylim(0.15, 2.15)\n",
    "plt.xlim(17.6, 27.8)\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "667acf5a",
   "metadata": {},
   "source": [
    "Well, that looks good -- the ratio of measured to true fluxes is centered on 1.0. It seems like the fluxes are recovered pretty well, on average.\n",
    "\n",
    "Congratulations! You have now learned how to compare measurements in the DP0.1 catalogs to the \"true\" simulated properties of objects. Have fun exploring more properties!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}