DP02_02b_Catalog_Queries_with_TAP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Catalog Queries with the TAP Service\n",
    "\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\">\n",
    "<br>\n",
    "Contact authors: Leanne Guy and Melissa Graham <br>\n",
    "Last verified to run: 2024-09-13 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_37 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: intermediate <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Description:** Execute complex ADQL queries with the TAP service. Visualize catalog data sets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Skills:** Use advanced ADQL and TAP functionality. Visualize retrieved datasets with bokeh and holoviews."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**LSST Data Products:** Object, ForcedSource, CcdVisit tables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Packages:** lsst.rsp, bokeh, pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Credit:**\n",
    "Originally developed by Leanne Guy in the context of the Rubin DP0.1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>. Questions are welcome as new topics in the <a href=\"https://community.lsst.org/c/support/dp0\">Support - Data Preview 0 Category</a> of the Rubin Community Forum. Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notebook provides an intermediate-level demonstration of how to use the Table Access Protocol (TAP) server and ADQL (Astronomy Data Query Language) to query and retrieve data from the DP0.2 catalogs.\n",
    "\n",
    "TAP provides standardized access to catalog data for discovery, search, and retrieval.\n",
    "Full <a href=\"http://www.ivoa.net/documents/TAP\">documentation for TAP</a> is provided by the International Virtual Observatory Alliance (IVOA).\n",
    "ADQL is similar to SQL (Structured Query Langage).\n",
    "The <a href=\"http://www.ivoa.net/documents/latest/ADQL.html\">documentation for ADQL</a> includes more information about syntax and keywords.\n",
    "Note that not all ADQL functionality is supported yet in the DP0-era RSP.\n",
    "\n",
    "**See the recommendations for TAP queries in DP0.2 tutorial 02a \"Introduction to the TAP Service\".**\n",
    "\n",
    "The [documentation for Data Preview 0.2](https://dp0-2.lsst.io/) includes definitions\n",
    "of the data products, descriptions of catalog contents, and ADQL recipes.\n",
    "\n",
    "### 1.1. Package imports\n",
    "\n",
    "Import general python packages, the Rubin TAP service utilities, and bokeh and holoviews for interactive visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "import bokeh\n",
    "from bokeh.io import output_file, output_notebook, show\n",
    "from bokeh.layouts import gridplot\n",
    "from bokeh.models import ColumnDataSource, CDSView, GroupFilter, HoverTool\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.transform import factor_cmap\n",
    "\n",
    "import holoviews as hv\n",
    "\n",
    "from lsst.rsp import get_tap_service, retrieve_query"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and parameters\n",
    "\n",
    "Instantiate the TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")\n",
    "assert service is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the maximum number of rows to display from pandas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "pandas.set_option('display.max_rows', 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Configure bokeh to generate output in notebook cells when show() is called."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create dictionaries for standard symbols, colors, and labels to use when\n",
    "plotting data in the six LSST filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_filter_labels = ['u', 'g', 'r', 'i', 'z', 'y']\n",
    "plot_filter_colors = {'u': '#0c71ff', 'g': '#49be61', 'r': '#c61c00',\n",
    "                      'i': '#ffc200', 'z': '#f341a2', 'y': '#5d0000'}\n",
    "plot_filter_symbols = {'u': 'o', 'g': '^', 'r': 'v', 'i': 's', 'z': '*', 'y': 'p'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Column manipulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the coordinates and radius to use for the example queries in Sections 2 and 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center_ra = 62\n",
    "center_dec = -37\n",
    "radius = 0.01\n",
    "\n",
    "str_center_coords = str(center_ra) + \", \" + str(center_dec)\n",
    "str_radius = str(radius)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start with the same query as used in the beginner TAP tutorial notebook 02a."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT coord_ra, coord_dec \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_center_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the query job asynchronously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Return the results as a `pandas` dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display `results`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. Rename columns\n",
    "\n",
    "Columns can be renamed by using `AS`.\n",
    "For example, use the same query as above but rename `coord_ra` and `coord_dec` \n",
    "as `ra` and `dec`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT coord_ra AS ra, \"\\\n",
    "        \"coord_dec AS dec \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_center_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, see the columns have been renamed in the `results` dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2. Apply functions\n",
    "\n",
    "A variety of mathematical, trigonometrical, and geometrical functions can be applied to columns,\n",
    "with the results stored as columns in the returned array.\n",
    "See the [ADQL documentation](https://www.ivoa.net/documents/latest/ADQL.html) for the full list.\n",
    "\n",
    "Below, use the same query as above but convert the coordinates from degrees to radians,\n",
    "and rename the columns `ra_radians` and `dec_radians`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT RADIANS(coord_ra) AS ra_radians, \"\\\n",
    "        \"RADIANS(coord_dec) AS dec_radians \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_center_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. Convert fluxes to magnitudes\n",
    "\n",
    "Photometric measurements for LSST objects are stored in the catalogs as fluxes with units of nanoJanskys (nJy).\n",
    "\n",
    "The conversion from nJy to [AB magnitude](https://en.wikipedia.org/wiki/AB_magnitude) is\n",
    "$m_{\\rm AB} = -2.5 \\log(f_{\\rm nJy}) + 31.4$.\n",
    "\n",
    "> **Warning:** Photometric measurements performed on difference images (images for which a \"template\" or \"reference\" image\n",
    "> of the same field has been subtracted) can be negative. Take care if converting fluxes from catalogs that are the result of difference image analysis (DIA).\n",
    "\n",
    "A function has been created to return fluxes as AB magnitudes: `scisql_nanojanskyToAbMab`.\n",
    "\n",
    "The query below converts the $g$-band `cModelFlux` to the $g$-band AB magnitude two ways,\n",
    "first using the equation (`gmag1`) and also using the function (`gmag2`), and the function\n",
    "for the error.\n",
    "See Section 5.1 for a description of the `cModel` flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT coord_ra, coord_dec, \"\\\n",
    "        \"-2.5 * LOG10(g_cModelFlux) + 31.4 AS gmag1, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(g_cModelFlux) AS gmag2, \"\\\n",
    "        \"scisql_nanojanskyToAbMagSigma(g_cModelFlux, g_cModelFluxErr) AS gmag2err \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_center_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4. Retrieve derived columns\n",
    "\n",
    "Here, a \"derived column\" refers to a column which has been create from other columns.\n",
    "\n",
    "The query below retrives the $gri$-band AB magnitudes, and\n",
    "derives the $g-r$ and $r-i$ colors and returns them as separate columns `grclr` and `riclr`.\n",
    "\n",
    "It is the convention in astronomy to define colors as the subtraction of the magnitudes, $c = m_1 - m_2$, not the fluxes.\n",
    "It is also the convention that the longer-wavelength filter is always subtracted from shorter-wavelength filter.\n",
    "In this convention, $m_1 - m_2 = c > 0$ indicates a \"redder\" object, and $c < 0$ is a \"bluer\" object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT scisql_nanojanskyToAbMag(g_cModelFlux) AS gmag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(r_cModelFlux) AS rmag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(i_cModelFlux) AS imag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(g_cModelFlux) - \"\\\n",
    "        \"scisql_nanojanskyToAbMag(r_cModelFlux) AS grclr, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(r_cModelFlux) - \"\\\n",
    "        \"scisql_nanojanskyToAbMag(i_cModelFlux) AS riclr \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_center_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a scatter plot to show the colors. Plot horizontal and vertical line at color 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 3))\n",
    "plt.plot(results['grclr'], results['riclr'],\n",
    "         'o', mew=0, ms=3, alpha=0.5, color='black')\n",
    "plt.axvline(0.0, lw=1, color='grey')\n",
    "plt.axhline(0.0, lw=1, color='grey')\n",
    "plt.xlim([-1, 2])\n",
    "plt.ylim([-2, 2])\n",
    "plt.xlabel('g-r color')\n",
    "plt.ylabel('r-i color')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 1: A color-color scatter plot with $g-r$ on the x-axis and $r-i$ on the y-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Column constraints\n",
    "\n",
    "In the query above, two constraints are used: the spatial constraint that a catalog object's\n",
    "coordinates must be within a certain radius of a defined point of interest, and that the\n",
    "object's `detect_isPrimary` flag must be `True`.\n",
    "\n",
    "This section demonstrates other types of column constraints.\n",
    "\n",
    "### 3.1. Comparison operators\n",
    "\n",
    "In the queries above, only the `=` operator is used,\n",
    "but others are available, e.g., `!=`, `<`, `>`, `<=`, and `>=`.\n",
    "\n",
    "Rewrite `detect_isPrimary = ` as `detect_isPrimary != 0`,\n",
    "and add a constraint that the $i$-band magnitude is less than or equal to 25 mag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT scisql_nanojanskyToAbMag(g_cModelFlux) AS gmag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(r_cModelFlux) AS rmag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(i_cModelFlux) AS imag \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_center_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary != 0 \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(i_cModelFlux) <= 25\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, see that adding the magnitude constraint has decreased the number of rows returned by the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Operators applied to functions\n",
    "\n",
    "The `DISTANCE` function will return the distance, in degrees, between two points.\n",
    "\n",
    "Use the coordinates of one of the objects returned by the queries above to \n",
    "define a point: `POINT('ICRS', 61.989258, -37.005119))`.\n",
    "Retrieve the distance of every other returned object from this object as `dist`.\n",
    "\n",
    "To furthermore demonstrate that operators can be used in both the `SELECT` and `WHERE`\n",
    "statements of an ADQL query, apply a constraint that objects must be less\n",
    "than 0.005 degrees away from the defined point.\n",
    "\n",
    "Sort the results based on `dist` by including an `ORDER BY` statement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT coord_ra, coord_dec, \"\\\n",
    "        \"DISTANCE(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"POINT('ICRS', 61.989258, -37.005119)) AS dist \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_center_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1 \"\\\n",
    "        \"AND DISTANCE(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"POINT('ICRS', 61.989258, -37.005119)) < 0.005 \"\\\n",
    "        \"ORDER BY dist ASC\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above, see that the first row has a very very small distance (a rounding error) from\n",
    "the defined point because it is the same object as was used to define the point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. Use of `WHERE ... BETWEEN ... AND ...`\n",
    "\n",
    "A constraint that one column's value is between two other columns' values\n",
    "can be applied with a statement like `<col1> BETWEEN <col2> AND <col3>`.\n",
    "The order matters.\n",
    "The statement `WHERE a BETWEEN b AND c` means $b < a < c$.\n",
    "Whereas `WHERE a BETWEEN c AND b` means $c < a < b$.\n",
    "\n",
    "Create a query with the constraint that the $r$-band `cModelFlux` value\n",
    "is in between the $g$- and $i$-band fluxes: $g < r < i$.\n",
    "This query will return red objects which have more flux in the redward filters ($r$ and $i$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT scisql_nanojanskyToAbMag(g_cModelFlux) AS gmag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(r_cModelFlux) AS rmag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(i_cModelFlux) AS imag \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_center_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1 \"\\\n",
    "        \"AND r_cModelFlux BETWEEN g_cModelFlux AND i_cModelFlux\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `pandas` dataframes, columns can be derived from other columns after the TAP query.\n",
    "\n",
    "Create columns for the $g-r$ and $r-i$ colors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results['grclr'] = results['gmag'] - results['rmag']\n",
    "results['riclr'] = results['rmag'] - results['imag']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a color-color scatter plot, as in Section 2.4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 3))\n",
    "plt.plot(results['grclr'], results['riclr'],\n",
    "         'o', mew=0, ms=3, alpha=0.5, color='black')\n",
    "plt.axvline(0.0, lw=1, color='grey')\n",
    "plt.axhline(0.0, lw=1, color='grey')\n",
    "plt.xlim([-1, 2])\n",
    "plt.ylim([-2, 2])\n",
    "plt.xlabel('g-r color')\n",
    "plt.ylabel('r-i color')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 2: The same as Figure 1, but only for objects with fluxes $f_g < f_r < f_i$.\n",
    "In other words, only \"red\" objects with $m_g - m_r > 0$ and $m_r - m_i > 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4. Use of `WHERE ... IN (...)`\n",
    "\n",
    "The constraint option of `WHERE <col1> IN (value1, value2, value3)` is\n",
    "a simpler way of doing, e.g., \n",
    "`WHERE <col1> = value1 OR <col1> = value2 OR <col1 = value3`.\n",
    "\n",
    "A scenario where the `WHERE ... IN ()` constraint might be used is when a\n",
    "colleague has provided a list of `objectId` that they have identified as\n",
    "targets of interest, and only the rows for those objects are desired.\n",
    "\n",
    "For example, the `objectId` for 12 objects which were identified as bright stars with similar\n",
    "$g-r$ and $i-z$ colors are used to define `id_array`.\n",
    "In the scenario where a colleague provided the list of `objectId`, the array\n",
    "could have been read in from a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_array = np.asarray((1249537790362809267, 1252528461990360512, 1248772530269893180,\n",
    "                       1251728017525343554, 1251710425339299404, 1250030371572068167,\n",
    "                       1253443255664678173, 1251807182362538413, 1252607626827575504,\n",
    "                       1249784080967440401, 1253065023664713612, 1325835101237446771),\n",
    "                      dtype='long')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_list = \"(\" + \", \".join(str(value) for value in id_array) + \")\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(id_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The query below retrieves the $griz$ `psfFlux` values, converted to AB magnitudes.\n",
    "The PSF fluxes are retrieved instead of `cModel` fluxes because these objects were\n",
    "already identified as stars, point sources, and so fluxes which assume the object\n",
    "has the shape of the point spread function (PSF) are appropriate.\n",
    "\n",
    "> **Warning:** Very long lists (tens of thousands) can be problematic for the TAP service, and are best avoided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT objectId, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(g_psfFlux) AS gmag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(r_psfFlux) AS rmag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(r_psfFlux) AS imag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(i_psfFlux) AS zmag \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE objectId IN \" + id_list\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Table joins\n",
    "\n",
    "Tables can be joined when they have an index or identifier column in common.\n",
    "\n",
    "There are four main types of table joins in ADQL/SQL:\n",
    " * inner: will return rows that appear in both tables and meet constraints (intersection set)\n",
    " * outer: will return rows from either table that meet constraints (union set)\n",
    " * left: will return rows from the first table that meet constraints, plus matching rows from the second table\n",
    " * right: will return rows from the second table that meet constraints, plus matching rows from the first table\n",
    "\n",
    "When `JOIN` is used in an ADQL query, it is an inner join by default.\n",
    "\n",
    "**Generic example:**\n",
    "To include a table join in an ADQL query requires the specification of which column to be \"joined on\".\n",
    "For example, the query below retrieves all rows from \"table1\" that are within\n",
    "0.002 degrees of coordinates 62.123, -37.456, but only if there is a row\n",
    "in \"table2\" which has the same value of `id` (i.e., an inner join, which is default). \n",
    "Columns `mag` and `band` from \"table2\" are also returned for these matching rows.\n",
    "\n",
    "```\n",
    "SELECT t1.id, t1.ra, t1.dec, t2.mag, t2.band\n",
    "FROM catalog.table1 AS t1\n",
    "JOIN catalog.table2 AS t2 ON t1.id = t2.id\n",
    "WHERE CONTAINS(POINT('IRCS', t1.ra, t1.dec), CIRCLE('ICRS', 62.123, -37.456 0.002)) = 1\n",
    "```\n",
    "\n",
    "**Object and Source examples:**\n",
    "a join between the DP0.2 `Object` and `ForcedSource` tables\n",
    "is demonstrated in Section 4.1.\n",
    "\n",
    "**Difference Image Analysis (DIA) examples:**\n",
    "table joins to create difference-image lightcurves for transient\n",
    "and variable objects are demonstrated in DP0.2 tutorial notebooks 07a and 07b.\n",
    "\n",
    "**DP0.2 Truth Match examples:**\n",
    "table joins with the DP0.2 `MatchesTruth` and `TruthSummary` tables are\n",
    "demonstrated in DP0.2 tutorial notebook 08."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1. Variable star lightcurve\n",
    "\n",
    "A common reason for a table join with DP0.2 is to create lightcurves.\n",
    "\n",
    "Examples of common joins:\n",
    " * `Object` and `ForcedSource` tables on column `objectId`\n",
    " * `Source` or `ForcedSource` and `CcdVisit` tables on column `ccdVisitId`\n",
    "\n",
    "Recall the contents of these tables:\n",
    " * `Source` table: measurements for detections in the processed visit images (PVIs)\n",
    " * `Object` table: measurements for detections in the deeply stacked images (coadds)\n",
    " * `ForcedSource` table: forced photometry on the PVIs at the location of all `Object`s\n",
    " * `CcdVisit` table: PVI metadata such a filter, date, exposure time\n",
    "\n",
    "Note that the `Source` table cannot be joined with the `ForcedSource` or `Object` tables,\n",
    "because it is the result of independent processing.\n",
    "\n",
    "First, define the coordinates of a known RR Lyrae variable star.\n",
    "Use a search radius of $\\sim2$\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_ra = 62.1479031\n",
    "target_dec = -35.799138\n",
    "radius = 0.0006\n",
    "\n",
    "str_target_coords = str(target_ra) + \", \" + str(target_dec)\n",
    "str_radius = str(radius)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1.1. Single-table query to find the objectId\n",
    "\n",
    "Find the `objectId` for the variable star with a simple, single-table query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT objectId \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_target_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retrieve the `results` and print the length of the dataframe.\n",
    "There should be only one row (one match within the small search radius)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the `objectId` of the first (and only) row of the `results` dataframe as `str_target_objectId`\n",
    "for use in Section 4.1.2 below.\n",
    "Also print it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "str_target_objectId = str(results['objectId'][0])\n",
    "print(str_target_objectId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1.2. Two-table join to retrieve the light curve\n",
    "\n",
    "Find all rows of the `ForcedSource` catalog which are measurements of this object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT scisql_nanojanskyToAbMag(fs.psfFlux) AS mag, \"\\\n",
    "        \"scisql_nanojanskyToAbMagSigma(fs.psfFlux, fs.psfFluxErr) AS magerr, \"\\\n",
    "        \"cv.band, cv.expMidptMJD \"\\\n",
    "        \"FROM dp02_dc2_catalogs.ForcedSource AS fs \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.CcdVisit AS cv ON fs.ccdVisitId = cv.ccdVisitId \"\\\n",
    "        \"WHERE fs.objectId = \" + str_target_objectId\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retrieve the `results` and print the length of the dataframe.\n",
    "The length will be the total number of individual observations of the object, in all filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the number of observations in each of the six filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values, counts = np.unique(results['band'], return_counts=True)\n",
    "for value, count in zip(values, counts):\n",
    "    print(value, count)\n",
    "del values, counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the $r$- and $i$-band light curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 3))\n",
    "for filt in ['r', 'i']:\n",
    "    tx = np.where(results['band'] == filt)[0]\n",
    "    plt.plot(results['expMidptMJD'][tx], results['mag'][tx],\n",
    "             plot_filter_symbols[filt], alpha=0.5, mew=0,\n",
    "             color=plot_filter_colors[filt], label=filt)\n",
    "plt.gca().invert_yaxis()\n",
    "plt.xlabel('Modified Julian Date')\n",
    "plt.ylabel('Apparent Magnitude')\n",
    "plt.legend(loc='upper left', handletextpad=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 3: The light curve for `objectId` = `1651589610221899038`, a known RR Lyrae variable star,\n",
    "in the $r$ (orange triangles) and $i$ (brown squares) filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. Three-table join to retrieve the light curve\n",
    "\n",
    "The single-table query for the `objectId` and the two-table join to retrieve\n",
    "the `ForcedSource` lightcurve can be combined into the following three-table join."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT scisql_nanojanskyToAbMag(fs.psfFlux) AS mag, \"\\\n",
    "        \"scisql_nanojanskyToAbMagSigma(fs.psfFlux, fs.psfFluxErr) AS magerr, \"\\\n",
    "        \"cv.band, cv.expMidptMJD \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object AS o \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.ForcedSource AS fs ON o.objectId = fs.objectId \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.CcdVisit AS cv ON fs.ccdVisitId = cv.ccdVisitId \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', o.coord_ra, o.coord_dec), \"\\\n",
    "        \"CIRCLE('ICRS', \" + str_target_coords + \", \" + str_radius + \")) = 1 \"\\\n",
    "        \"AND o.detect_isPrimary = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As these are the exact same data as retrieved in Section 4.1.2 above,\n",
    "plotting the light curve is left as an exercise for the learner (Section 6)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()\n",
    "del query, results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Interactive plots for large data sets\n",
    "\n",
    "In this section, the **bokeh** package is used to create interactive plots to explore a dataset.\n",
    "Multiple panels are used to show different representations of the same dataset.\n",
    "A selection applied to either panel will highlight the selected points in the other panel.\n",
    "\n",
    "<a href=\"https://bokeh.org/\">Bokeh Documentation</a> <br>\n",
    "<a href=\"https://holoviews.org/\">Holoviews Documentation</a>\n",
    "\n",
    "This section uses the same center coordiantes as Section 2 and 3, but instead of a small radius\n",
    "of 0.01 degrees, a full degree radius is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_center_coords = \"62, -37\"\n",
    "use_radius = \"1.0\"\n",
    "\n",
    "center_coords = SkyCoord(62, -37, frame='icrs', unit='deg')\n",
    "search_radius = 1.0*u.deg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1. Create the query\n",
    "\n",
    "Select rows of the `Object` table within the central search area defined above.\n",
    "As always with the `Object` catalog, enforce that the `detect_isPrimary` flag is \"True\" ($=1$).\n",
    "Only return objects which have $g$, $r$, and $i$-band magnitudes between 17 and 23 mag (using the $<$ and $>$ operators),\n",
    "and which have a measured r-band extendedness parameter (i.e., it is not NaN or NULL).\n",
    "Extendedness is typically a parameter between 0 (point-like) and 1 (extended); for DP0.2 it is either 0 or 1.\n",
    "\n",
    "Return the `cModelFlux` measurements as magnitudes.\n",
    "Composite Model (CModel) fluxes are similar in nature to [those measured for SDSS](https://www.sdss3.org/dr8/algorithms/magnitudes.php#cmodel).\n",
    "In short, it is the linear combination of the best fit exponential and de Vaucouleurs profiles.\n",
    "The `cModel` type of flux is appropriate to use for extended objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT objectId, detect_isPrimary, \" + \\\n",
    "        \"coord_ra AS ra, coord_dec AS dec, \" + \\\n",
    "        \"scisql_nanojanskyToAbMag(g_cModelFlux) AS mag_g_cModel, \" + \\\n",
    "        \"scisql_nanojanskyToAbMag(r_cModelFlux) AS mag_r_cModel, \" + \\\n",
    "        \"scisql_nanojanskyToAbMag(i_cModelFlux) AS mag_i_cModel, \" + \\\n",
    "        \"r_extendedness \" + \\\n",
    "        \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \" + \\\n",
    "        \"CIRCLE('ICRS', \" + use_center_coords + \", \" + use_radius + \")) = 1 \" + \\\n",
    "        \"AND detect_isPrimary = 1 \" + \\\n",
    "        \"AND scisql_nanojanskyToAbMag(g_cModelFlux) > 17.0 \" + \\\n",
    "        \"AND scisql_nanojanskyToAbMag(g_cModelFlux) < 23.0 \" + \\\n",
    "        \"AND scisql_nanojanskyToAbMag(r_cModelFlux) > 17.0 \" + \\\n",
    "        \"AND scisql_nanojanskyToAbMag(r_cModelFlux) < 23.0 \" + \\\n",
    "        \"AND scisql_nanojanskyToAbMag(i_cModelFlux) > 17.0 \" + \\\n",
    "        \"AND scisql_nanojanskyToAbMag(i_cModelFlux) < 23.0 \" + \\\n",
    "        \"AND r_extendedness IS NOT NULL \"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2. Prepare the results for plotting\n",
    "\n",
    "The basis for any data visualization is the underlying data.\n",
    "\n",
    "The `ColumnDataSource` (CDS) is the core of bokeh plots.\n",
    "A CDS can be prepared from the data returned by a query, enabling it to be passed directly to bokeh.\n",
    "Bokeh automatically creates a CDS from data passed as python lists or numpy arrays.\n",
    "CDSes are useful as they allow data to be shared between multiple plots and renderers, enabling brushing and linking.\n",
    "A CDS is essentially a collection of sequences of data that have their own unique column name. \n",
    "\n",
    "Getting the data preparation phase right is the key to creating powerful visualizations. \n",
    "\n",
    "Recall the center coordinates that were defined with `astropy` `SkyCoord`.\n",
    "Get the central RA and Dec as values (datatype float) and print them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "center_ra = center_coords.ra.deg\n",
    "center_dec = center_coords.dec.deg\n",
    "print(center_ra, center_dec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a python dictionary that associates the `r_extendedness` column with the Object being a star or a galaxy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "object_map = {0.0: 'star', 1.0: 'galaxy'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a python dictionary to store the data from the query and pass to the `ColumnDataSource`.\n",
    "All columns in a CDS must have the same length."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data = dict(ra=results['ra'], dec=results['dec'],\n",
    "            target_ra=results['ra']-center_ra,\n",
    "            target_dec=results['dec']-center_dec,\n",
    "            gmi=results['mag_g_cModel']-results['mag_i_cModel'],\n",
    "            gmag=results['mag_g_cModel'],\n",
    "            rmag=results['mag_r_cModel'],\n",
    "            imag=results['mag_i_cModel'])\n",
    "source = ColumnDataSource(data=data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additional data can be added to the CDS after creation, as\n",
    "demonstrated in the cell below with the addition of `objectId` and `extendedness` in the r-band.\n",
    "\n",
    "Note that the `objectId` needs to be converted to a string, \n",
    "because bokeh can't interpret an integer value that high."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objIdString = np.array(results['objectId']).astype('str')\n",
    "source.data['objectId'] = objIdString\n",
    "source.data['r_extendedness'] = results['r_extendedness']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `object_map` dictionary to add an `object_type` column to the CDS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "source.data['object_type'] = results['r_extendedness'].map(object_map)\n",
    "source.data['object_type']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3. Plot a color-magnitude diagram\n",
    "\n",
    "Use bokeh to plot a color-magnitude (g vs. g-i) diagram making use of the `cModel` magnitudes.\n",
    "Hover over the points in the plot to see their values.\n",
    "\n",
    "Define the plot aesthetics and tools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_options = {'height': 400, 'width': 400,\n",
    "                'tools': ['box_select', 'reset', 'box_zoom', 'help']}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the hover tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tooltips = [\n",
    "    (\"Col (g-i)\", \"@gmi\"),\n",
    "    (\"Mag (g)\", \"@gmag\"),\n",
    "    (\"Mag (r)\", \"@rmag\"),\n",
    "    (\"Mag (i)\", \"@imag\"),\n",
    "    (\"Type\", \"@object_type\")\n",
    "]\n",
    "hover_tool_cmd = HoverTool(tooltips=tooltips)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a color-magnitude diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "p = figure(title=\"Color - Magnitude Diagram\",\n",
    "           x_axis_label='g-i', y_axis_label='g',\n",
    "           x_range=(-1.8, 4.3), y_range=(23.5, 16),\n",
    "           **plot_options)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Color-code the different object types by defining a color map palette."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "object_type_palette = ['darkred', 'lightgrey']\n",
    "p.add_tools(hover_tool_cmd)\n",
    "p.scatter(x='gmi', y='gmag', source=source,\n",
    "          size=3, alpha=0.6,\n",
    "          legend_field=\"object_type\",\n",
    "          color=factor_cmap('object_type',\n",
    "                            palette=object_type_palette,\n",
    "                            factors=['star', 'galaxy']),\n",
    "          hover_color=\"darkblue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the interactive plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 4: The $g$-band apparent magnitude vs. the $g-i$ color for galaxies (grey) and stars (maroon). The hover text includes the color, magnitudes, and type for each point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.4. Plot a color-color (r-i vs. g-r) diagram\n",
    "\n",
    "Add a color-color (r-i vs. g-r) diagram and make use of the advanced linking features of bokeh \n",
    "to enable brushing and linking between the the color-magnitude diagram and this color-color plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "source.data['rmi'] = results['mag_r_cModel'] - results['mag_i_cModel']\n",
    "source.data['gmr'] = results['mag_g_cModel'] - results['mag_r_cModel']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CMD above is very crowded. Below, filter on object type and plot stars only.\n",
    "\n",
    "Use a GroupFilter to select rows from the CDS that satisfy `object_type` = \"star\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "stars = CDSView()\n",
    "stars.filter = GroupFilter(column_name='object_type', group=\"star\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the various options for the plot, and create the hover tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_options = {'height': 350, 'width': 350,\n",
    "                'tools': ['box_zoom', 'box_select',\n",
    "                          'lasso_select', 'reset', 'help']}\n",
    "\n",
    "hover_tool = HoverTool(tooltips=[(\"(RA,DEC)\", \"(@ra, @dec)\"),\n",
    "                                 (\"(g-r,g)\", \"(@gmr, @gmag)\"),\n",
    "                                 (\"objectId\", \"@objectId\"),\n",
    "                                 (\"type\", \"@object_type\")])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the three plots: spatial, color-magnitude, and color-color, then plot all three together on a grid.\n",
    "\n",
    "Create the spatial plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "title_spatial = 'Spatial centred on (RA,DEC) = '+use_center_coords\n",
    "\n",
    "fig_spatial = figure(title=title_spatial,\n",
    "                     x_axis_label=\"Delta RA\", y_axis_label=\"Delta DEC\",\n",
    "                     **plot_options)\n",
    "fig_spatial.scatter(x='target_ra', y='target_dec',\n",
    "                    source=source, view=stars,\n",
    "                    size=4.0, alpha=0.6,\n",
    "                    color='teal', hover_color='firebrick')\n",
    "fig_spatial.add_tools(hover_tool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the color-magnitude plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_cmag = figure(title=\"Color-Magnitude Diagram\",\n",
    "                  x_axis_label=\"g-r\", y_axis_label=\"g\",\n",
    "                  x_range=(-1.0, 2.0), y_range=(23.5, 16),\n",
    "                  **plot_options)\n",
    "fig_cmag.scatter(x='gmr', y='gmag', source=source, view=stars,\n",
    "                 size=4.0, alpha=0.6,\n",
    "                 color='teal', hover_color='firebrick')\n",
    "fig_cmag.add_tools(hover_tool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the color-color plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_cc = figure(title=\"Color-Color Diagram\",\n",
    "                x_axis_label=\"g-r\", y_axis_label=\"r-i\",\n",
    "                x_range=(-1.0, 2.0), y_range=(-1.0, 2.5),\n",
    "                **plot_options)\n",
    "fig_cc.scatter(x='gmr', y='rmi', source=source, view=stars,\n",
    "               size=4.0, alpha=0.6,\n",
    "               color='teal', hover_color='firebrick')\n",
    "fig_cc.add_tools(hover_tool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show all three plots together in a row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "p = gridplot([[fig_spatial, fig_cmag, fig_cc]])\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 5: Three side-by-side panels showing points for stars only (teal). Left: the difference in right ascension (Delta RA) vs. the difference in declination (Delta DEC) between object coordinates and 62, -37 degrees. Center: the $g$-band apparent magnitude vs. the $g-r$ color. Right: the $r-i$ color vs. the $g-r$ color. The hover text is the same for all panels: the coordinates, the $g-r$ color and $g$ mag, `objectId`, and type. \n",
    "\n",
    "Hover the mouse over the data points in any of the plots in order to \n",
    "use the hover tool to see information about individual data points (e.g., the `objectId`).\n",
    "Notice the data points highlighted in red on one panel with the hover tool are also highlighted on the other panels.\n",
    "\n",
    "Click on the box select or lasso select icons in the upper right corner of the figure. \n",
    "Use the box select or lasso select functionality to make a selection in any panel by clicking and dragging to define a region. \n",
    "The selected data points will be displayed in darker teal in all panels.\n",
    "\n",
    "Use the refresh icon to reset the plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del results, data, p, stars, source"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Exercises for the learner\n",
    "\n",
    "(1) In Section 3.3, use `WHERE ... BETWEEN ... AND ...` to create a query to retrieve\n",
    "blue objects (instead of red) which have $g-r<0$ and $r-i<0$.\n",
    "Hint: `r_cModelFlux BETWEEN i_cModelFlux AND g_cModelFlux`\n",
    "\n",
    "(2) In Section 5.1, rewrite the query to use `WHERE ... BETWEEN ... AND ...` in order\n",
    "to restrict the magnitudes of the rows returned. Use the demonstration in Section 2.4\n",
    "to retrieve the $g-r$ and $r-i$ colors as part of the query, instead of creating\n",
    "them in Section 5.4.\n",
    "\n",
    "(3) The query below has three errors. Try running it, and see that the job fails with\n",
    "phase ERROR. Use `job.raise_if_error()` to figure out what the three errors are.\n",
    "This is the same three-table join query as was used in Section 4.2, just with \n",
    "three errors added.\n",
    "\n",
    "```\n",
    "SELECT scisql_nanojanskyToAbMag(fs.psfFlux) AS mag, scisql_nanojanskyToAbMag(fs.psfFlux, fs.psfFluxErr) AS magerr, cv.band, cv.expMidptMJD FROM dp02_dc2_catalogs.Object AS o JOIN dp02_dc2.catalogs.ForcedSource AS fs ON o.objectId = fs.objectId JOIN dp02_dc2_catalogs.CcdVisit AS cv ON fs.ccdVisitId = cv.visitId WHERE CONTAINS(POINT('ICRS', o.coord_ra, o.coord_dec), CIRCLE('ICRS', 62.1479031, -35.799138, 0.0006)) = 1 AND o.detect_isPrimary = 1\n",
    "```\n",
    "\n",
    "(4) In Section 4.2, the light curve is not plotted because the same data was retrieved\n",
    "as in Section 4.1.2. Go back and plot the data retrieved in Section 4.2, but this time\n",
    "show all six LSST filters in the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}