02_Intermediate_TAP_Query.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250, style=\"padding: 10px\"> \n",
    "<b>Intermediate TAP Queries for DP0 catalogs </b> <br>\n",
    "Last verified to run on 2021-10-12 with LSST Science Pipelines release w_2021_40 <br>\n",
    "Contact authors: Leanne Guy <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",
   "metadata": {},
   "source": [
    "**Credit:** Originally developed by Leanne Guy in the context of the Rubin DP0.1. Please consider acknowledging Leanne Guy if this notebook is used for the preparation of journal articles or software releases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "\n",
    "The Rubin Science Platform provides QUERY access to the DP0.1 catalogs via TAP from jupyter notebooks. TAP is a Virtual Observatory protocol for access to catalog data. In this tutorial, we will learn how to explore the DP0.1 archive via TAP and execute complex queries to retrieve data. Full TAP documentation can be found [here](https://www.ivoa.net/documents/TAP/).\n",
    "\n",
    "This notebook demonstrates how to: <br>\n",
    "1. Explore the DP0.1 schema and catalogs using the Rubin TAP service<br>\n",
    "2. Query the DP0.1 Object and Truth Match catalogs, and retrieve data for analysis<br>\n",
    "3. Visualize and interact with the retrieved data set<br>\n",
    "4. Work with asynchronous TAP queries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Up "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import general python packages\n",
    "import numpy as np\n",
    "import re\n",
    "import pandas\n",
    "from pandas.testing import assert_frame_equal\n",
    "\n",
    "# Import the Rubin TAP service utilities\n",
    "from lsst.rsp import get_tap_service, retrieve_query\n",
    "\n",
    "# Astropy\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "# Bokeh for interactive visualization\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",
    "# Set the maximum number of rows to display from pandas\n",
    "pandas.set_option('display.max_rows', 20)\n",
    "\n",
    "# Configure bokeh to generate output in notebook cells when show() is called.\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ignore warnings\n",
    "import warnings\n",
    "from astropy.units import UnitsWarning\n",
    "warnings.simplefilter(\"ignore\", category=UnitsWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general the order of results from database queries cannot be assumed to be the same every time.\n",
    "This function sorts the data so we can compare the result dataframes even if the records are not in the same order from the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sort_dataframe(df, sort_key='objectid'):\n",
    "    df = df.sort_values('objectId')\n",
    "    df.set_index(np.array(range(len(df))), inplace=True)  # Since we are sorting, we need to reset the incremental index as well\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Explore the DP0.1 schema "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1 Create the Rubin TAP Service client\n",
    "\n",
    "Table Access Procotol (TAP) provides standardized access to catalog data for discovery, search, and retrieval. Full <a href=\"http://www.ivoa.net/documents/TAP\">documentation for TAP</a> is provided by the International Virtual Observatory Alliance (IVOA).\n",
    "\n",
    "The TAP service uses a query language similar to SQL (Structured Query Langage) called ADQL (Astronomical Data Query Language). The <a href=\"http://www.ivoa.net/documents/latest/ADQL.html\">documentation for ADQL</a> includes more information about syntax and keywords.\n",
    "\n",
    "**Hazard Warning:** Not all ADQL functionality is supported yet in the DP0 RSP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get an instance of the TAP service\n",
    "service = get_tap_service()\n",
    "assert service is not None\n",
    "assert service.baseurl == \"https://data.lsst.cloud/api/tap\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2 Schema discovery\n",
    "\n",
    "To find out what schemas, tables and columns exist, we will query the Rubin TAP schema. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Query to find out what schemas are in the Rubin TAP_SCHEMA\n",
    "query = \"SELECT * FROM tap_schema.schemas\"\n",
    "\n",
    "# Execute the query\n",
    "results = service.search(query)\n",
    "\n",
    "# A TAP Results object is returned\n",
    "print(type(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the results to an astropy table and display\n",
    "results = service.search(query).to_table()\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.3 DC2 tables (catalogs) in DP0.1\n",
    "\n",
    "All the DP0 tables (catalogs) are in the \"dp01_dc2_catalogs\" schema (table collection). We can programatically extract the DP0.1 schema_name into a variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the DP0 schema name and store as a variable\n",
    "schema_names = results['schema_name']\n",
    "for name in schema_names:\n",
    "    if re.search('dp01', name):\n",
    "        dp01_schema_name = name\n",
    "        break\n",
    "print(\"DP0.1 schema is \" + dp01_schema_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's explore tables in the DP0.1 schema, ordering them by their database.  This is the order in which they will appear presented to the user in the RSP Portal. We see the five tables in the DP0.1 schema, the same five tables that are presented via the Portal GUI, together with a description of each. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the query to explore the tables in the DP0.1 schema\n",
    "query = \"SELECT * FROM tap_schema.tables \"\\\n",
    "        \"WHERE tap_schema.tables.schema_name = '\" \\\n",
    "        + dp01_schema_name + \"' order by table_index ASC\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the query\n",
    "results = service.search(query)\n",
    "results = results.to_table()\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "Here are some definitions to help delegates understand the contents of the TAP schema. \n",
    "\n",
    "* `schema` - database terminology for the abstract design that represents the storage of data in a database. \n",
    "* `tap_schema` - the specific schema describing the TAP service. All TAP services must support a set of tables in a schema named TAP_SCHEMA that describe the tables and columns included in the service.\n",
    "* `table` - a collection of related data held in a table format in a database, e.g., the object(dp01_dc2_catalogs.object) or position (dp01_dc2_catalogs.position) tables \n",
    "* `table collection` - a collection of tables. e.g., `dp01_dc2_catalogs`\t\n",
    "* `results` - the query result set. The TAP service returns data from a query as a `TAPResults` object. Find more about `TAPResults` [here](https://pyvo.readthedocs.io/en/latest/api/pyvo.dal.TAPResults.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Querying the DP0.1 Object and Truth Match catalogs\n",
    "\n",
    "The Object catalog (dp01_dc2_catalogs.object) contains sources detected in the coadded images (also called stacked or combined images). The Object catalog is likely to be the catalog that is of the most interest to DP0 delgates. \n",
    "\n",
    "The `object` catalog is described in the <a href=\"https://arxiv.org/abs/2101.04855\">DESC's DC2 data release note</a>, and more information about the simulated data can be found in the <a href=\"https://ui.adsabs.harvard.edu/abs/2021ApJS..253...31L/abstract\">DESC's DC2 paper</a>. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1 Specifying the maximum number of records to return\n",
    "For debugging and testing queries, it is often useful to only return a few records for expediency. This can be done in one of two ways, setting the `TOP` field in a query, or setting the `maxrec` parameter in the TAP service query. The two methods are identical. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the maximum records to return\n",
    "max_rec = 5\n",
    "\n",
    "# Build a query to find object with extendedness = 0 and sort the returned\n",
    "# result set by decreasing magnitude in the r band.\n",
    "# Only return the first 5 results\n",
    "query = \"SELECT TOP \" + str(max_rec) + \\\n",
    "        \" objectId, ra, dec, extendedness, mag_r, magerr_r, good \" \\\n",
    "        \"FROM dp01_dc2_catalogs.object \" \\\n",
    "        \"WHERE extendedness = 0 \" \\\n",
    "        \"AND mag_r < 24 \" \\\n",
    "        \"ORDER by mag_r DESC\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the query\n",
    "results = service.search(query)\n",
    "assert len(results) == max_rec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the same query using the maxrec parameter instead of the TOP\n",
    "query = \"SELECT objectId, ra, dec, extendedness, mag_r, magerr_r, good \" \\\n",
    "        \"FROM dp01_dc2_catalogs.object \" \\\n",
    "        \"WHERE extendedness = 0 \" \\\n",
    "        \"AND mag_r < 24 \" \\\n",
    "        \"ORDER by mag_r DESC\"\n",
    "results1 = service.search(query, maxrec=max_rec)\n",
    "assert len(results1) == max_rec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the results to pandas data frames and  assert that the\n",
    "# contents of the two tables are identical\n",
    "assert_frame_equal(sort_dataframe(results.to_table().to_pandas()),\n",
    "                   sort_dataframe(results1.to_table().to_pandas()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2 Cone search around a point with specified radius\n",
    "\n",
    "We will execute a cone search on the Object table centered on (RA, Dec) = (62.0, -37.0) with a radius of 0.1 degrees and applying a cut on magnitude.\n",
    "We expect to get 15,670 results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define our reference position on the sky and cone radius in arcseconds\n",
    "# to use in all following examples\n",
    "coord = SkyCoord(ra=62.0*u.degree, dec=-37.0*u.degree, frame='icrs')\n",
    "radius = 0.1 * u.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT ra, dec, mag_g, mag_i \" \\\n",
    "        \"mag_i, mag_g_cModel, mag_r_cModel, mag_i_cModel, \" \\\n",
    "        \"psFlux_g, psFlux_r, psFlux_i, \" \\\n",
    "        \"cModelFlux_g, cModelFlux_r, cModelFlux_i, \" \\\n",
    "        \"tract, patch, extendedness, good, clean \" \\\n",
    "        \"FROM dp01_dc2_catalogs.object \" \\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', ra, dec),CIRCLE('ICRS', \" \\\n",
    "        + str(coord.ra.value) + \", \" + str(coord.dec.value) + \", \" \\\n",
    "        + str(radius.value) + \" )) = 1\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For more detailed analysis of results, converting\n",
    "# to a pandas dataframe is often very useful\n",
    "results = service.search(query).to_table().to_pandas()\n",
    "\n",
    "# Use an assertion to make sure we got the correct number of results.\n",
    "assert len(results) == 15670"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Cone search joining the results with the truth infomation \n",
    "\n",
    "We will now join the results from the above query with the truth information.  We will also add in some quality filters on the match. Objects in the truth-match table that do not have matches in the object table have “match_objectId = -1,” while those with legitimate matches contain the ‘objectId’ of the corresponding object from the object table in “match_objectId.” By requiring this to be greater than or equal to zero, we extract only objects with matches. We also keep only sources satisfying the “is_good_match” flag, which is described in the schema as being “True if this object–truth matching pair satisfies all matching criteria.” (Note that “1” and “TRUE” are equivalent in ADQL.)\n",
    "\n",
    "With these additional quality filters applied to the matching with the truth information, we only get 14424 results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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),CIRCLE('ICRS', \" \\\n",
    "        + str(coord.ra.value) + \", \" + str(coord.dec.value) + \", \" \\\n",
    "        + str(radius.value) + \" )) = 1 \" \\\n",
    "        \"AND truth.match_objectid >= 0 \"\\\n",
    "        \"AND truth.is_good_match = 1\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = service.search(query).to_table().to_pandas()\n",
    "assert len(results) == 14424"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> If you see a pink-highlighted \"Note: NumExpr detected...\" message after executing the next cell, know that it is not a warning and it is safe to ingore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many of each type in the dataset\n",
    "# The 'truth_type' in the truth_match table is 1= galaxies, 2=stars, 3=SNe\n",
    "n_stars = results[results[\"truth_type\"] == 2].shape[0]\n",
    "print(f'There are {n_stars} stars out of a total of {len(results)}')\n",
    "print(f'There are {results[results[\"truth_type\"] == 1].shape[0]} galaxies')\n",
    "print(f'There are {results[results[\"truth_type\"] == 3].shape[0]} SNe')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Simple histogram to bin categorical data. \n",
    "\n",
    "Now we will create a simple categorical histogram of the number of each truth_type in the dataset. We will use the 'GROUP BY' ADQL command to group the Objects in the truth_match catalog by type (1: galaxies, 2:stars, 3: SNe), and the 'COUNT' command to count the number of Objects in each category. Finally we will use  the 'ORDER' command to order the results by ascending order of truth_type. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "query_histogram = \"SELECT truth_type, count(truth_type) \" \\\n",
    "                  \" FROM dp01_dc2_catalogs.truth_match \" \\\n",
    "                  \" GROUP BY truth_type \" \\\n",
    "                  \" ORDER BY truth_type\"\n",
    "print(query_histogram)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "object_type_histogram = service.search(query_histogram).to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Map the numerical values for each type to a more descriptive name\n",
    "object_map = {1: 'galaxy', 2: 'star', 3: 'SNe'}\n",
    "object_type_histogram['truth_type'] = \\\n",
    "    object_type_histogram['truth_type'].map(object_map)\n",
    "object_type_histogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Visualize and analyse the results\n",
    "\n",
    "Now we will do some interactive analysis with the data we have above. We will use bokeh to create interactive plots so that we can explore the dataset, using multiple panels showing different representations of the same dataset. 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>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1 Data preparation\n",
    "The basis for any data visualization is the underlying data. We will prepare ColumnDataSource (CDS) from the data returned by the query above that can be passed directly to bokeh. The CDS is the core of bokeh plots. Bokeh automatically creates a CDS from data passed as python lists or numpy arrays.  CDS are useful as they allow data to be shared between multiple plots and renderers, enabling brushing and linking.  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 key to creating powerful visualizations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a python dictionary to store the data from the\n",
    "# query and pass to the ColumnDataSource\n",
    "# All columns in a CDS must have the same length\n",
    "data = dict(ra=results['ra'], dec=results['dec'],\n",
    "            target_ra=results['ra']-coord.ra.value,\n",
    "            target_dec=results['dec']-coord.dec.value,\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",
    "            )\n",
    "source = ColumnDataSource(data=data)\n",
    "\n",
    "# Additional data can be added to the Column Data Source after creation\n",
    "source.data['objectId'] = results['objectId']\n",
    "\n",
    "# We will want to filter on the truth type later\n",
    "# We will convert the truth_type integer to a more descriptive string\n",
    "source.data['truth_type'] = results['truth_type'].map(object_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the truth_type has been updated\n",
    "source.data['truth_type']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2 Color-Magnitude Diagram \n",
    "We will use bokeh to plot a color-magnitude (g vs. g-i) diagram making use of the cModel magnitudes. Hover over the points \n",
    "in the plot to see their values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the plot asthetics and tools\n",
    "plot_options = {'plot_height': 400, 'plot_width': 400,\n",
    "                'tools': ['box_select', 'reset', 'box_zoom', 'help']}\n",
    "# Define the hover tool\n",
    "tooltips = [\n",
    "    (\"Col (g-i)\", \"@gmi\"),\n",
    "    (\"Mag (g)\", \"@gmag\"),\n",
    "    (\"Mag (r)\", \"@rmag\"),\n",
    "    (\"Mag (i)\", \"@imag\"),\n",
    "    (\"Type\", \"@truth_type\")\n",
    "]\n",
    "hover_tool_cmd = HoverTool(tooltips=tooltips)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Colour-Magnitude Diagram, color coding the different object types\n",
    "p = figure(title=\"Colour - Magnitude Diagram\",\n",
    "           x_axis_label='g-i', y_axis_label='g',\n",
    "           x_range=(-1.8, 4.3), y_range=(29.5, 16),\n",
    "           **plot_options)\n",
    "\n",
    "# Define a palette for the truth types\n",
    "truth_type_palette = ['darkred', 'lightgrey', 'blue']\n",
    "p.add_tools(hover_tool_cmd)\n",
    "p.circle(x='gmi', y='gmag', source=source,\n",
    "         size=3, alpha=0.6,\n",
    "         legend_field=\"truth_type\",\n",
    "         color=factor_cmap('truth_type',\n",
    "                           palette=truth_type_palette,\n",
    "                           factors=['star', 'galaxy', 'SNe']),\n",
    "         hover_color=\"firebrick\")\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3 Color-color (r-i vs. g-r) diagram. \n",
    "We will add a color-color (r-i vs. g-r) diagram and make use of the advanced linking features of bokeh to enable brushing and linking between the the color-magnitude diagram and this color-color plot. The CMD in 4.2 is very crowded as it contains 14424 data points. We will now filter on the truth-type to plot stars only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will now add some additional data to our data structure to\n",
    "# create a colour-colour diagram\n",
    "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use a GroupFilter to select rows from the\n",
    "# CDS that satisfy 'truth_type' stars\n",
    "stars = CDSView(source=source,\n",
    "                filters=[GroupFilter(column_name='truth_type', group=\"star\")])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define various options for the plot\n",
    "plot_options = {'plot_height': 350, 'plot_width': 350,\n",
    "                'tools': ['box_zoom', 'box_select',\n",
    "                          'lasso_select', 'reset', 'help']}\n",
    "\n",
    "\n",
    "# Create the hover tool for these plots\n",
    "hover_tool = HoverTool(tooltips=[(\"(RA,DEC)\", \"(@ra, @dec)\"),\n",
    "                                 (\"(g-r,g)\", \"(@gmr, @gmag)\"),\n",
    "                                 (\"objectId\", \"@objectId\"),\n",
    "                                 (\"type\", \"@truth_type\")])\n",
    "\n",
    "\n",
    "# Spatial plot\n",
    "title_spatial = f'Spatial centred on (RA,DEC) = \\\n",
    "({coord.ra.value},{coord.dec.value})'\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.circle(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)\n",
    "\n",
    "# Colour magnitude plot\n",
    "fig_cmag = figure(title=\"Colour-Magnitude Diagram\",\n",
    "                  x_axis_label=\"g-r\", y_axis_label=\"g\",\n",
    "                  x_range=(-1.0, 3.5), y_range=(29.5, 16),\n",
    "                  **plot_options)\n",
    "fig_cmag.circle(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)\n",
    "\n",
    "# Colour colour plot\n",
    "fig_cc = figure(title=\"Colour-Colour Diagram\",\n",
    "                x_axis_label=\"g-r\", y_axis_label=\"r-i\",\n",
    "                x_range=(-1.0, 3.5), y_range=(-1.0, 3.5),\n",
    "                **plot_options)\n",
    "fig_cc.circle(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)\n",
    "\n",
    "# Plot all three on a grid\n",
    "p = gridplot([[fig_spatial, fig_cmag, fig_cc]])\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the hover tool to see information about individual datapoints (e.g., the object_id). \n",
    "This information should appear automatically as you hover the mouse over the datapoints in any of the plots.\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 selection box icon (with a \"+\" sign) or the selection lasso icon found in the upper right corner of the figure. \n",
    "Use the selection box and selection lasso to make various selections in either panel by clicking and dragging on either panel. \n",
    "The selected data points will be displayed in the other panel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.4 Comparison with truth information \n",
    "Let's compare the measurements from the Object table to the “true” values for some objects. \n",
    "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,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another way that data can be added to the CDS\n",
    "source.add(results['cModelFlux_r'] / results['flux_r'], name='flux_ratio')\n",
    "source.add(results['truth_mag_r'], name='truth_mag_r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a plot of the distribution of objects on\n",
    "# sky using the ColumnDataSource's two columns\n",
    "plot_options = {'plot_height': 300, 'plot_width': 800,\n",
    "                'tools': ['box_zoom', 'box_select', 'reset', 'help']}\n",
    "\n",
    "# Define the tooltips and create a HoverTool instance\n",
    "tooltips = [\n",
    "    (\"Measured/True flux\", \"@flux_ratio\"),\n",
    "    (\"Mag (r)\", \"@truth_mag_r\"),\n",
    "    (\"Type\", \"@truth_type\")\n",
    "]\n",
    "hover_tool = HoverTool(tooltips=tooltips)\n",
    "\n",
    "# Create a new figure\n",
    "p = figure(title=\"Measured/true flux vs true magnitude\",\n",
    "           x_range=(16, 29.5), y_range=(0.15, 2.4),\n",
    "           x_axis_label='r magnititude (truth)',\n",
    "           y_axis_label='Measured flux / True flux (r band)',\n",
    "           **plot_options)\n",
    "# Add a circle renderer defining several attributes\n",
    "p.circle(x='truth_mag_r', y='flux_ratio',\n",
    "         size=3, alpha=0.5, source=source,\n",
    "         legend_field=\"truth_type\",\n",
    "         color=factor_cmap('truth_type',\n",
    "                           palette=truth_type_palette,\n",
    "                           factors=['star', 'galaxy', 'SNe']))\n",
    "\n",
    "# Add the hover tool created above to the figure\n",
    "p.add_tools(hover_tool)\n",
    "\n",
    "# Display\n",
    "show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.0 Asynchronous TAP queries\n",
    "\n",
    "So far, we have executed all queries as synchronous queries. This means that the query will continue executing in the notebook until it is finished. You can see when the Jupyter cell is running by the asterisk to the left of the  cell. For synchronous queries, the cell will continue to run until the query completes and the results are returned. The asterisk will then become a number. This is a good option for short queries that take order seconds to minutes.\n",
    "\n",
    "For longer queries, or for running multiple queries at the same time, an asynchronous query may be more suitable. Asynchronous queries allow you to execute more python while the query runs on the database. Results can be retrieved later on. This is especially important for queries that are long or may return a lot of results. It is also safeguards long queries against network outages or timeouts. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.1 Submit and run a job\n",
    "We will use the cone search joining the results with the truth infomation query from 2.3 and confirm that the results from the asynchronous query are the same as from the synchronous query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create and submit the job. This step does not run the query yet\n",
    "job = service.submit_job(query)\n",
    "\n",
    "# # Get the job URL\n",
    "print('Job URL is', job.url)\n",
    "\n",
    "# Get the job phase. It will be pending as we have not yet started the job\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the job. You will see that the the cell completes executing,\n",
    "# even though the query is still running\n",
    "job.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use this to tell python to wait for the job to finish if\n",
    "# you don't want to run anything else while waiting\n",
    "# The cell will continue executing until the job is finished\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A usefull funtion to raise an exception if there was a problem with the query\n",
    "job.raise_if_error()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Once the job completes successfully, you can fetch the results\n",
    "async_results = job.fetch_result()\n",
    "\n",
    "# Assert that the results are the same as obtained from\n",
    "# executing synchronous queries\n",
    "assert len(async_results) == 14424 \n",
    "assert_frame_equal(sort_dataframe(results), sort_dataframe(async_results.to_table().to_pandas()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.2 Retrieving the results from a previous asynchronous job\n",
    "Job results may still be available from previously run queries. You can retrieve these results if you know the URL of the job.\n",
    "This includes jobs executed in the Portal. You can retrieve the URL for the query and retrieve the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "retrieved_job = retrieve_query(job.url)\n",
    "retrieved_results = retrieved_job.fetch_result()\n",
    "assert len(retrieved_results) == 14424\n",
    "assert_frame_equal(sort_dataframe(retrieved_results.to_table().to_pandas()),\n",
    "                   sort_dataframe(async_results.to_table().to_pandas()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4.3 Deleting a job\n",
    "Once the job is finished and you have retrieved your results, you can delete the job and the results from the server. The results will be deleted automatically after a period of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job.delete()"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}