DP03_06_User_Uploaded_Tables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "# User Uploaded Tables\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: Yumi Choi <br>\n",
    "Last verified to run: 2024-12-18 <br>\n",
    "LSST Science Piplines version: Weekly 2024_50 <br>\n",
    "Container Size: small <br>\n",
    "Targeted learning level: Beginner <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** Use the TAP upload functionality for user-supplied tables and join them with DP0.3 catalogs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Use the TAP service to upload a table and join it to an LSST table with ADQL."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** TAP tables dp03_catalogs_10yr.SSObject, dp03_catalogs_10yr.MPCORB, dp03_catalogs_10yr.DiaSource"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:** lsst.rsp.get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credit:**\n",
    "Developed by Yumi Choi. This tutorial is based on <a href=\"https://dp0-3.lsst.io/tutorials-dp0-3/portal-dp0-3-5.html\">a Portal tutorial</a> by Christina Williams for using user-supplied tables in queries for DP0.3 and <a href=\"https://github.com/rubin-dp0/cst-dev/blob/main/MLG_sandbox/DP03/gaia_hack_day_Feb2024.ipynb\">a Jupyter Notebook</a> by Melissan Graham and Jake Kurlander for accessing Gaia data and matching with DP0.3 data for solar system objects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a400d1bb-e1ec-4160-bd41-bc7ca953a97c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-27T17:55:39.523317Z",
     "iopub.status.busy": "2023-07-27T17:55:39.522902Z",
     "iopub.status.idle": "2023-07-27T17:55:40.035942Z",
     "shell.execute_reply": "2023-07-27T17:55:40.034503Z",
     "shell.execute_reply.started": "2023-07-27T17:55:39.523290Z"
    },
    "tags": []
   },
   "source": [
    "**Get Support:** 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."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cfc73be0",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notebook illustrates the process of uploading user-provided tables via the TAP service and integrating them into queries for DP0.3. It focuses on two types of user-provided tables: (1) tables generated outside the RSP, and (2) tables created within the RSP but retrieved from an external database using [PyVO](https://github.com/astropy/pyvo)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bf722b4-8227-441f-9f3e-0f8140040fc1",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1.1. Import packages\n",
    "\n",
    "Import general python packages and the Rubin Table Access Protocol (TAP) service.\n",
    "\n",
    "[PyVO](https://pyvo.readthedocs.io/en/latest/) is a package providing access to remote data and services of the Virtual observatory (VO) using Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d4078a1-8232-4f11-97ec-2a67ca19308e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import getpass\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pyvo\n",
    "from astropy.table import Table\n",
    "import pandas as pd\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c217adff-25ed-4fce-95e7-8aa04630f6cc",
   "metadata": {},
   "source": [
    "### 1.2. Define parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a4a8b3b-91cf-4729-a707-99b7b11a5b8d",
   "metadata": {},
   "source": [
    "Set a few style parameters for the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c780a382-5c5e-4f20-9ffd-7aa517c37146",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.style.use('tableau-colorblind10') \n",
    "params = {'axes.labelsize': 12,\n",
    "          'font.size': 12,\n",
    "          'legend.fontsize': 10}\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c5690df-b8c2-4d3f-9e47-b18aabb661c1",
   "metadata": {},
   "source": [
    "For the LSST filter set ugrizy, always use symbols to represent the filters in addition to color."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b4dac0c-6db0-423e-bffc-b3fd8917d451",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_filter_colors_white_background = {\n",
    "    'u': '#0c71ff',\n",
    "    'g': '#49be61',\n",
    "    'r': '#c61c00',\n",
    "    'i': '#ffc200',\n",
    "    'z': '#f341a2',\n",
    "    'y': '#5d0000'\n",
    "}\n",
    "\n",
    "plot_symbols = {\n",
    "    'u': 'o',\n",
    "    'g': '^',\n",
    "    'r': 'v',\n",
    "    'i': 's',\n",
    "    'z': '*',\n",
    "    'y': 'p'\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9ab2c64-b43b-4311-88a4-810e5bafb7dc",
   "metadata": {},
   "source": [
    "Define the path to to input files, which are stored in a shared directory in the Notebook Aspect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e349bb36-a9f4-4207-9563-83882db7e0bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_path = '/project/cst_repos/tutorial-notebooks-data/data/'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12ab5a66-b37d-4a3b-b207-b95a68c7a7a7",
   "metadata": {},
   "source": [
    "Define the path to the home directory, where outputs generated by this tutorial will be saved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f60d078e-a8b1-4925-8368-fe41064a1fa4",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_home_dir = '/home/' + getpass.getuser() + '/'\n",
    "print(my_home_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00b268e9-b011-4c57-9604-91e82bca58ce",
   "metadata": {},
   "source": [
    "Start the TAP service and assert that it exists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b62bc6f9-62ea-46a8-a05d-f6a7efaa855c",
   "metadata": {},
   "outputs": [],
   "source": [
    "rsp_tap = get_tap_service(\"ssotap\")\n",
    "assert rsp_tap is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec51ac0b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Spatial and temporal cross-match to diaSources\n",
    "\n",
    "The scientific scenario used for this demonstration is to answer the question\n",
    "of whether LSST detected a moving object which, based on its known orbit, was\n",
    "expected to appear at given coordinates on nine different nights.\n",
    "\n",
    "The list of expected coordinates and dates is provided in a file with nine rows\n",
    "(i.e., nine expected locations and times), with the format\n",
    "expected for user-uploaded tables: column names in the first row with no # symbol at the start of the row.\n",
    "The list contains columns of identifier index, right ascension, declination, and modified julian date (id, ra, dec, and mjd).\n",
    "\n",
    "The [ADQL](https://www.ivoa.net/documents/ADQL/) query below uploads this list to the TAP service and cross-matches to\n",
    "the table of detected sources that results from difference image analysis: the `diaSource` table\n",
    "(learn more about [DP0.3 catalogs](https://dp0-3.lsst.io/data-products-dp0-3/index.html#dp0-3-data-products-definition-document-dpdd))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43f597c0-8ba4-40c7-aee9-9f250dac82dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnm1 = input_path + '/dp03_06_user_table_1.cat'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c62755f2-f3fc-4759-9f10-e6f72d245edb",
   "metadata": {},
   "source": [
    "**Option:** View the contents of the file with `more`, or return the word count (lines, words, characters) with `wc`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81b8cd59-1ba3-4eaa-846f-6478ed0c3cf5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# os.system('more ' + fnm1)\n",
    "# os.system('wc ' + fnm1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf968de5-ce59-41c9-88f0-83cd88bdc083",
   "metadata": {},
   "source": [
    "Read the file as user table 1, `ut1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "220c90e6-5399-46c4-a1dc-49219c7f46fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "ut1 = Table.read(fnm1, format='ascii.basic')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e99ebf49-6bf9-47bf-bcf6-5b768be279d3",
   "metadata": {},
   "source": [
    "Define a query to cross match the coordinates\n",
    "and dates in the file to detections in the DP0.3 `diaSource` catalog.\n",
    "\n",
    "This query is applying a spatial threshold of 10 arcseconds (0.00278 degrees), and a temporal threshold of half a day.\n",
    "\n",
    "Notice that the columns from the user-uploaded table are being renamed in the query, e.g., `ut1.ra AS ut1_ra`.\n",
    "If this renaming is not done, the TAP service will rename the columns from the user table as `ra2` and `dec2`\n",
    "to distinguish them from the `ra` and `dec` columns from the `DiaSource` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "184ee9e9-c038-4aad-bbe1-87367794e63d",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "        SELECT dias.ra, dias.dec, dias.midPointMjdTai, dias.ssObjectId,\n",
    "        ut1.ra AS ut1_ra, ut1.dec AS ut1_dec, ut1.mjd AS ut1_mjd, ut1.id AS ut1_id\n",
    "        FROM dp03_catalogs_10yr.DiaSource AS dias, TAP_UPLOAD.ut1 AS ut1\n",
    "        WHERE CONTAINS(POINT('ICRS', dias.ra, dias.dec),\n",
    "        CIRCLE('ICRS', ut1.ra, ut1.dec, 0.00278))=1\n",
    "        AND ABS(dias.midPointMjdTai - ut1.mjd) < 0.5\n",
    "        ORDER BY dias.ssObjectId\n",
    "        \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46e8419f-ca60-4fc4-abf0-fa255f6d13ad",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-03T03:44:00.969710Z",
     "iopub.status.busy": "2024-06-03T03:44:00.969017Z",
     "iopub.status.idle": "2024-06-03T03:44:01.592200Z",
     "shell.execute_reply": "2024-06-03T03:44:01.590939Z",
     "shell.execute_reply.started": "2024-06-03T03:44:00.969685Z"
    }
   },
   "source": [
    "Create the job by submitting the query and then run it asynchronously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12d0080b-13f2-4a50-8b0e-126aa8b02faa",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = rsp_tap.submit_job(query, uploads={\"ut1\": ut1})\n",
    "job.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3e7aaf0-2f7b-4185-9d1f-f5eb40ae06d3",
   "metadata": {},
   "source": [
    "Check that the job is completed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "801f6215-f541-408e-82c1-e147267c8765",
   "metadata": {},
   "outputs": [],
   "source": [
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5383c21-d6c5-43b9-b189-62e429d11298",
   "metadata": {},
   "source": [
    "Retrieve the results and display them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4164a28a-0786-4f0d-ada6-c7a5a371477e",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "476b0ef8-a1d8-4cfb-b408-715949807ae0",
   "metadata": {},
   "source": [
    "In the table above, notice that the `ssObjectId` is all the same.\n",
    "This is because the file was created to contain the detections of a single moving object across multiple nights.\n",
    "\n",
    "Notice also that whereas the user-uploaded table had 9 rows, there are 15 matches to the `diaSource` table.\n",
    "Six of the rows had two matches.\n",
    "This means that six times, this moving object was detected twice within 10 arcseconds in the same night.\n",
    "\n",
    "### 2.1. Plot the spatial and temporal offsets\n",
    "\n",
    "For the purposes of demonstrating some kind of plot out of the results of this cross-matching, calculate the spatial and temporal\n",
    "offsets (`s_off` and `t_off`) between the detected `diaSource`s and the user uploaded expectations,\n",
    "and then plot the spatial vs. the temporal offsets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c452a798-221a-4443-8061-69d26bc33f25",
   "metadata": {},
   "outputs": [],
   "source": [
    "cos_dec = np.cos(np.deg2rad(results['dec']))\n",
    "delta_ra = cos_dec * np.abs(results['ra'] - results['ut1_ra'])\n",
    "delta_dec = np.abs(results['dec'] - results['ut1_dec'])\n",
    "s_off = 3600.0 * np.sqrt(delta_ra**2 + delta_dec**2)\n",
    "t_off = 24.0 * 60.0 * np.abs(results['midPointMjdTai'] - results['ut1_mjd'])\n",
    "del cos_dec, delta_ra, delta_dec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3245ed75-b3a1-478f-be36-e4a7732c2936",
   "metadata": {},
   "source": [
    "Moving objects do not exhibit a constant rate of motion on the sky in terms of, e.g., arcsec per minute,\n",
    "but a correlation between spatial and temporal offset is expected.\n",
    "Keep in mind that nine of these points will be at zero (no offset; perfect match).\n",
    "But for the six additional, newly-discovered detections in\n",
    "the `diaSource` catalog which were not in the user-uploaded table,\n",
    "there should be larger spatial offsets for larger temporal offsets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3bede14-c8f1-41b2-987a-bb00343668a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "plt.plot(t_off, s_off, 'o', ms=10, alpha=0.3, mew=0)\n",
    "plt.xlabel('time difference (minutes)')\n",
    "plt.ylabel('sky distance (arcseconds)')\n",
    "plt.title('distance between same-night observations')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7cba04f-1538-417f-8260-80193b092b3a",
   "metadata": {},
   "source": [
    "> Figure 1: Above, the sky distance in arcseconds is plotted versus the time difference in minutes for the 15 detections cross-matched with the 9 anticipated coordinates and times in the user-uploaded table.\n",
    "\n",
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc70f4d4-9551-435f-aa1c-84da684414b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "del t_off, s_off\n",
    "del fnm1, ut1, query, job, results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94714005-d795-48ea-b80a-2885c09b21e3",
   "metadata": {},
   "source": [
    "## 3. Object identifier cross-match to diaSources\n",
    "\n",
    "This section demonstrates how to upload a user-supplied table and join it with a DP0.3 table.\n",
    "\n",
    "In this scenario, a list of identifiers (`ssObjectId`) for moving objects has been assembled\n",
    "by the user and stored in a file (one column, two rows of data).\n",
    "\n",
    "As in Section 2, define the file name, optionally print the contents or the word count, \n",
    "and then read in the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55688974-a2cc-496b-801a-e1dfe38223e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fnm2 = input_path + 'dp03_06_user_table_2.cat'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66188120-7185-4838-920a-644473c99845",
   "metadata": {},
   "outputs": [],
   "source": [
    "# os.system('more ' + fnm2)\n",
    "# os.system('wc ' + fnm2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c82b41af-37db-4c27-839f-bccf87e8d0cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "ut2 = Table.read(fnm2, format='ascii.basic')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d17f1f77-92b3-4a16-be5a-035115e0b1f2",
   "metadata": {},
   "source": [
    "Define a query to use an inner join on the user-uploaded table and the `diaSource` table\n",
    "to return only records (objects) that have matching values in both tables.\n",
    "\n",
    "**Inner join:** An inner join is the default type of join, and is what is done if simply `JOIN` is used in a query.\n",
    "Other types of joins include outer joins (e.g., left, right, full) which optionally return\n",
    "all records in one or the other table, with matches where available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0c62f66-0498-438c-b436-12ba79000d34",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "        SELECT ut2.ssObjectId_user, dias.ssObjectId, dias.ra, dias.dec\n",
    "        FROM TAP_UPLOAD.ut2 as ut2\n",
    "        INNER JOIN dp03_catalogs_10yr.diaSource as dias\n",
    "        ON ut2.ssObjectId_user = dias.ssObjectId\n",
    "        \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b459919e-c47c-43fa-a0f4-208b92ff8e17",
   "metadata": {},
   "source": [
    "Run the query asynchronously and then retrieve the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22229072-b16b-4ade-9f8c-e7489c52ebdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = rsp_tap.submit_job(query, uploads={\"ut2\": ut2})\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5024461b-4a4d-4961-b893-4654e7cf679c",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b2af35e-2fb6-403e-a9f0-2e607f7e551d",
   "metadata": {},
   "source": [
    "Print the unique `ssObjectId` and how many detections each had."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3bf5e0f-8269-426e-a38e-d80106b279ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "uniqueIds, counts = np.unique(results['ssObjectId'], return_counts=True)\n",
    "for uniqueId, count in zip(uniqueIds, counts):\n",
    "    print(\"ssObjectId:\", uniqueId, \"Number of detections:\", count)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d66aa55-aca3-4703-a6a1-376b9fdb5da8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-05-07T05:18:45.699109Z",
     "iopub.status.busy": "2024-05-07T05:18:45.698800Z",
     "iopub.status.idle": "2024-05-07T05:18:45.704297Z",
     "shell.execute_reply": "2024-05-07T05:18:45.703757Z",
     "shell.execute_reply.started": "2024-05-07T05:18:45.699088Z"
    }
   },
   "source": [
    "### 3.1. Plot the sky distribution of LSST detections\n",
    "\n",
    "Create a plot showing sky distribution of the detections of these two unique objects over 10 years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b6651e9-4f82-47d7-80a6-a863a9702fd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "for uniqueId in uniqueIds:\n",
    "    tx = results['ssObjectId'] == uniqueId\n",
    "    plt.plot(results['ra'][tx], results['dec'][tx],\n",
    "             'o', ms=5, alpha=0.3, mew=0, label=str(uniqueId))\n",
    "plt.xlabel('Right Ascension [deg]')\n",
    "plt.ylabel('Declination [deg]')\n",
    "plt.title('LSST detections of two moving objects')\n",
    "plt.legend(loc='lower left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1c27c6b-f22e-4993-93c4-10413484cebe",
   "metadata": {},
   "source": [
    "> Figure 2: Above, the coordinates of the LSST difference-image detections from the `diaSource` catalog for the two moving objects listed in the user-uploaded table.\n",
    "\n",
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cbbe7d7-156f-4c46-8657-cb578c452bdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "del fnm2, ut2, query, job, results, uniqueIds, counts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28a8f930-4423-4677-aa35-42faa2300019",
   "metadata": {},
   "source": [
    "## 4. Gaia data for DP0.3 asteroids\n",
    "\n",
    "This section demonstrates how to do an external TAP query to retrieve a Gaia data table using `PyVO` from the RSP,\n",
    "and then upload the table and use it in a joint query with the LSST DP0.3 catalogs.\n",
    "\n",
    "Gaia data is also accessible via the [Gaia Archive](https://www.cosmos.esa.int/web/gaia-users/archive)\n",
    "(see the Help tab for more information about Gaia catalogs and data access). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bbad6ce-afde-44a4-9c1d-77e5074d67b4",
   "metadata": {},
   "source": [
    "### 4.1. Retrieve Gaia data for main-belt asteroids\n",
    "\n",
    "Get an instance of the Gaia TAP service using `PyVO` and assert that it exists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7644df2-4ac8-466a-83f8-7a896258f4fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_tap_url = 'https://gea.esac.esa.int/tap-server/tap'\n",
    "gaia_tap = pyvo.dal.TAPService(gaia_tap_url)\n",
    "assert gaia_tap is not None\n",
    "assert gaia_tap.baseurl == gaia_tap_url"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1f90a3b-db1f-41c7-bd48-d53d5584068f",
   "metadata": {},
   "source": [
    "Query the Gaia database for main-belt asteroids (MBAs) following the population definition used by the\n",
    "[JPL Horizons small body database query tool](https://ssd.jpl.nasa.gov/tools/sbdb_query.html):\n",
    "2.0 < `a` < 3.25 au and `q` > 1.666 au.\n",
    "\n",
    "To expedite the query, only return objects with more than 200 observations (i.e., well-observed objects),\n",
    "and limit the number of objects returned to 1000 with `SELECT TOP 1000`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c7a64ae-e845-42c8-9f1c-d1c85b7cf99e",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_query = \"\"\"\n",
    "             SELECT TOP 1000 denomination, inclination,\n",
    "                 eccentricity, semi_major_axis\n",
    "             FROM gaiadr3.sso_orbits\n",
    "             WHERE num_observations > 200\n",
    "             AND semi_major_axis > 2.0\n",
    "             AND semi_major_axis < 3.25\n",
    "             AND semi_major_axis*(1-eccentricity) > 1.666\n",
    "             \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93bdd35f-602c-42e0-8556-79bbb32727f6",
   "metadata": {},
   "source": [
    "Run the query and retrieve the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dee88e3-3aa8-411a-bc97-eb57b293db1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_job = gaia_tap.submit_job(gaia_query)\n",
    "gaia_job.run()\n",
    "gaia_job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', gaia_job.phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e590d892-326f-4abf-a0fb-61b9e15102e2",
   "metadata": {},
   "source": [
    "Retrieve the results and convert them to an `astropy` table with `to_table`, and then\n",
    "to a `pandas` table with `to_pandas`, in order to enable some column reformatting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "680fb0e8-7feb-4abf-83b4-9aeb2a4e06e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_results = gaia_job.fetch_result().to_table().to_pandas()\n",
    "gaia_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13ea1769-d3ad-47e3-a8ee-2d50db607ec4",
   "metadata": {},
   "source": [
    "#### 4.1.1. Plot the orbital parameters\n",
    "\n",
    "To explore the retrieved data for 1000 MBAs, make plots showing eccentricity (`e`) and inclination (`i`) against semi-major axis (`a`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "786b5319-f458-4c9b-8b75-c0d0c22eeff8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 3))\n",
    "plt.subplot(121)\n",
    "plt.scatter(gaia_results['semi_major_axis'], gaia_results['eccentricity'],\n",
    "            s=3, alpha=0.3)\n",
    "plt.xlabel('a [au]')\n",
    "plt.ylabel('e')\n",
    "plt.subplot(122)\n",
    "plt.scatter(gaia_results['semi_major_axis'], np.sin(gaia_results['inclination']),\n",
    "            s=3, alpha=0.3)\n",
    "plt.xlabel('a [au]')\n",
    "plt.ylabel(r'sin($i$)')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37e5c09f-3333-49ad-b26e-acb3c25a6889",
   "metadata": {},
   "source": [
    "> Figure 3: Above, scatter plots of the orbital parameters for 1000 MBAs from Gaia. At left, eccentricity (e) versus the semi-major axis (a). At right, inclination versus semi-major axis.\n",
    "\n",
    "#### 4.1.2. Reformat the MPC designation column\n",
    "\n",
    "The `MPCORB` table in DP0.3 follows the standard format for MPC designation, which is the year, then a space,\n",
    "and then an identifier that is several characters, capital letters and numbers.\n",
    "However, the `denomination` in some Gaia tables (like `sso_orbit`, queried above) uses underscores instead of\n",
    "spaces and/or lower-case letters instead of capital letters.\n",
    "\n",
    "To ensure accurate cross-matching (table joining) between the Gaia data and the DP0.3 data in Section 4.2.2,\n",
    "the `denomination` column must first be converted to the standard MPC designation format (`mpc_desig1`).\n",
    "A [known issue](https://dp0-3.lsst.io/data-products-dp0-3/data-simulation-dp0-3.html#known-issues)\n",
    "with the `mpcDesignation` column in the `MPCORB` table can lead to incorrect cross-matches for objects with `mpcDesignation` values longer than 8 characters. To avoid this problem, the `fullDesignation`\n",
    "column from the `MPCORB` table will be used in Section 4.2.2, despite its own issue with an unnecessary prefix of \"2011\".\n",
    "\n",
    "Additionally, the `denomination` column needs to be converted to use spaces instead of underscores, while retaining\n",
    "the lowercase letters (`mpc_desig2`), to correctly identify a moving object by name in the Gaia `sso_observation`\n",
    "table in Section 4.2.3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78eaedd6-a954-494c-bf1b-b7417a5311ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_results['mpc_desig1'] = '2011 '+ gaia_results['denomination'].str.replace('_', ' ').str.upper()\n",
    "gaia_results['mpc_desig2'] = gaia_results['denomination'].str.replace('_', ' ')\n",
    "gaia_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74a10b29-be39-4de0-adae-eff499aab84c",
   "metadata": {},
   "source": [
    "Save the result as a CSV (comma separated values) file in the home directory to explore an option of using `pandas` table joins in Section 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce31ab29-ce68-429a-b391-1c23fc467e91",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_results.to_csv(my_home_dir+'gaia_mba.cat', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e97f10b2-caf6-4532-9a38-345e4634bed9",
   "metadata": {},
   "source": [
    "Clean up. Only `gaia_results` is needed for the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c28a0fa-5424-4edf-a119-4002e2ba3c41",
   "metadata": {},
   "outputs": [],
   "source": [
    "del gaia_tap_url, gaia_tap\n",
    "del gaia_query, gaia_job"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af945265-bfe6-40c0-86dc-b5232a103af3",
   "metadata": {},
   "source": [
    "### 4.2. Gaia object cross-match to MPCORB and diaSources\n",
    "\n",
    "The scientific scenario used for this demonstration is that for the 1000 Gaia moving objects,\n",
    "the question is whether LSST detected them and, if so, what their $g$-band difference-image magnitudes\n",
    "are for each detection.\n",
    "\n",
    "#### 4.2.1. Reformat the table\n",
    "\n",
    "As `gaia_results` is a `pandas` table, but the TAP service does not accept uploads in a pandas `DataFrame` format,\n",
    "convert it back to an `astropy` table using the imported `Table` method.\n",
    "Follow the convention established in the sections above and name the table `ut3`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3422b55-7f4d-44ea-b1d2-0cd78a2bb569",
   "metadata": {},
   "outputs": [],
   "source": [
    "ut3 = Table.from_pandas(gaia_results)\n",
    "ut3.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6053f80-5dab-457a-b89b-8e5132c67bdf",
   "metadata": {},
   "source": [
    "Notice that the `denomination` and `mpc_desig` columns are of data type (`dtype`) Unicode (`U`).\n",
    "Uploading a table which includes Unicode columns will not work.\n",
    "Convert the Unicode columns (`dtype.kind='U'`) to bytestring (`dtype.kind='S'`) using the\n",
    "`convert_unicde_to_bytestring`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9411d0f0-434f-410d-838e-110a83dd6575",
   "metadata": {},
   "outputs": [],
   "source": [
    "ut3.convert_unicode_to_bytestring()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f70651c-2138-45dd-9caa-ea5b4bfdcab8",
   "metadata": {},
   "source": [
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6692d544-a0a1-4590-9692-a3e8a142b88c",
   "metadata": {},
   "outputs": [],
   "source": [
    "del gaia_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6732902a-008d-4e82-bfc4-fa737c6ad8e3",
   "metadata": {},
   "source": [
    "#### 4.2.2. Execute the query\n",
    "\n",
    "Prepare a query to join three tables: the user-specified Gaia table data, the DP0.3 `MPCORB` table, and the DP0.3 `diaSource` table.\n",
    "\n",
    "This query will return the magnitudes for all LSST detections in $g$-band difference images from the DP0.3 `diaSource` table\n",
    "for all objects in the user-uploaded Gaia table data (where a match exists).\n",
    "\n",
    "In the query, the first join between the Gaia user-uploaded table and the DP0.3 `MPCORB` table obtains the internal LSST `ssObjectId`\n",
    "for a given moving object by joining on MPC designation.\n",
    "The second join between the `MPCORB` table and the `diaSource` table obtains all $g$-band detections associated\n",
    "with the `ssObjectId`.\n",
    "This triple join is necessary because the `diaSource` table does not contain the MPC designation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd1b8d0e-e598-4027-a88c-9dd237c9fd42",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "        SELECT mpc.fullDesignation, mpc.ssObjectId,\n",
    "        ut3.denomination, ut3.mpc_desig2,\n",
    "        dias.midPointMjdTai, dias.mag, dias.band\n",
    "        FROM dp03_catalogs_10yr.MPCORB AS mpc\n",
    "        INNER JOIN TAP_UPLOAD.ut3 as ut3\n",
    "        ON ut3.mpc_desig1 = mpc.fullDesignation\n",
    "        INNER JOIN dp03_catalogs_10yr.DiaSource as dias\n",
    "        ON dias.ssObjectId = mpc.ssObjectId\n",
    "        Where dias.band = 'g'\n",
    "        \"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff3daf9a-ec3c-4c06-90ca-3f81f1c24cbb",
   "metadata": {},
   "source": [
    "Execute the query, then retrieve the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5faa4914-5b54-42dd-ac00-5643378fc348",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = rsp_tap.submit_job(query, uploads={\"ut3\": ut3})\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f0d4655-0da1-4c23-bb1a-08b46c29afef",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table()\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e642b90-93a9-42bb-8fc1-e88f007ca7a0",
   "metadata": {},
   "source": [
    "Print the number of objects (out of the 1000 from Gaia) that were detected by LSST,\n",
    "and also the minimum, maximum, and mean number of LSST detections over the ten-year survey."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2658ae97-fd23-4700-a384-df81d3b86d73",
   "metadata": {},
   "outputs": [],
   "source": [
    "uniqueIds, counts = np.unique(results['ssObjectId'], return_counts=True)\n",
    "print('Number of LSST-detected objects: ', len(uniqueIds))\n",
    "print('Minimum, maximum, and mean number of detections per object: ',\n",
    "      np.min(counts), np.max(counts), int(np.mean(counts)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e451849b-224f-4d36-9a42-9d5d340961a3",
   "metadata": {},
   "source": [
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61961c95-1cb7-4b7a-9d74-69d9f08e4e7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "del ut3, query, job\n",
    "del uniqueIds, counts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecc33d2d-b7fc-46a0-86f8-4efa45f0f936",
   "metadata": {},
   "source": [
    "#### 4.2.3. Plot the Gaia+LSST light curve for one MBA\n",
    "\n",
    "Select one random MBA and save its data entry as a small table called `random_mba`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fb292f4-55aa-417f-a872-e6e3829277be",
   "metadata": {},
   "outputs": [],
   "source": [
    "uniqueObj = np.random.choice(results['ssObjectId'], 1)[0]\n",
    "print('ssObjectId: ', uniqueObj)\n",
    "random_mba = results[results['ssObjectId'] == uniqueObj]\n",
    "print('MPC designation: ', random_mba['mpc_desig2'][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b45e4f7-acf5-463c-8120-58955f0e2111",
   "metadata": {},
   "source": [
    "Retrieve Gaia photometry for the selected MBA. \n",
    "\n",
    "Occasionally, the following error message may occur: \n",
    "```DALServiceError: 401 Client Error: 401 for url: https://gea.esac.esa.int/tap-server/tap/sync```.\n",
    "This usually represents an intermittent service error.\n",
    "To avoid it, rerun the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75ac700b-ee9c-4523-803c-91df0442161d",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_tap_url = 'https://gea.esac.esa.int/tap-server/tap'\n",
    "gaia_tap = pyvo.dal.TAPService(gaia_tap_url)\n",
    "assert gaia_tap is not None\n",
    "assert gaia_tap.baseurl == gaia_tap_url"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26b889a3-c6a5-410f-abb0-3c325d28d769",
   "metadata": {},
   "source": [
    "Define and execute a query to retrieve the dates (`epoch_utc`) and $g$-band magnitudes (`g_mag`)\n",
    "of Gaia observations for the randomly selected MBA by only returning rows\n",
    "where the `denomination` is equal to `mpc_desig2`, which is the MPC\n",
    "Designation reformatted to the convention of the Gaia `sso_observation` table\n",
    "(see Section 4.1.2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5872a57e-c18e-4cb6-8e38-37703ca58d65",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_query = \"\"\"\n",
    "             SELECT epoch_utc, g_mag\n",
    "             FROM gaiadr3.sso_observation\n",
    "             WHERE denomination = '{}'\n",
    "             \"\"\".format(random_mba['mpc_desig2'][0])\n",
    "gaia_job = gaia_tap.submit_job(gaia_query)\n",
    "gaia_job.run()\n",
    "gaia_job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', gaia_job.phase)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9bfa2c7-b930-4a32-af20-4495b4cd44a9",
   "metadata": {},
   "source": [
    "Retrieve the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cfe4c6a-5c0b-43b4-83ed-76d86450977d",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_results = gaia_job.fetch_result().to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bdafd75-62d8-42c0-af2d-41b7e15909c2",
   "metadata": {},
   "source": [
    "Convert Universal Time Coordinated (UTC) observation date (`epoch_utc`) to a Modified Julian Date (MJD) observation date.\n",
    "\n",
    "From the [Gaia documentation](https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_sso_tables/ssec_dm_sso_observation.html), `epoch_utc` is the Gaiacentric epoch UTC, while the LSST records observing time in MJD. The conversion is MJD = UTC + 55197.5 (day)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9cbe7ba-d81c-4d0c-98a1-4496ed1009d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_results['epoch_mjd'] = gaia_results['epoch_utc'] + 55197.5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d80558ce-b77b-4840-9423-9243bcc251c0",
   "metadata": {},
   "source": [
    "Plot the $g$-band light curve for the selected asteroid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "794fc254-9063-49da-9cde-580bce4ca65b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "plt.plot(gaia_results['epoch_mjd'], gaia_results['g_mag'],\n",
    "         's', color='dodgerblue', label='Gaia G')\n",
    "plt.plot(random_mba['midPointMjdTai'], random_mba['mag'],\n",
    "         plot_symbols['g'],\n",
    "         color=plot_filter_colors_white_background['g'], label='LSST g')\n",
    "ymin = min(gaia_results['g_mag'].min(), random_mba['mag'].min())\n",
    "ymax = max(gaia_results['g_mag'].max(), random_mba['mag'].max())\n",
    "plt.ylim(ymax+0.2, ymin-0.2)\n",
    "plt.legend(loc='upper left')\n",
    "plt.xlabel('MJD [days]')\n",
    "plt.ylabel('magnitude')\n",
    "plt.title('MPC designation:' + random_mba['mpc_desig2'][0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16471be3-2837-4365-910f-ca155beca412",
   "metadata": {},
   "source": [
    "> Figure 4: Above, the Gaia G-band detections in cyan squares and the simulated DP0.3 LSST $g$-band magnitudes in green circles as a function of time. The variation in apparent magnitude over time is primarily due to the object's changing distance from Earth.\n",
    "\n",
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abe33ff4-283e-4aff-8ab9-b989cf9fa962",
   "metadata": {},
   "outputs": [],
   "source": [
    "del results, uniqueObj, random_mba\n",
    "del gaia_tap_url, gaia_tap, gaia_query, gaia_job, gaia_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d14acce7-08ed-4a2b-ab68-47d88f66356e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-07-09T22:18:01.759735Z",
     "iopub.status.busy": "2024-07-09T22:18:01.759090Z",
     "iopub.status.idle": "2024-07-09T22:18:01.764494Z",
     "shell.execute_reply": "2024-07-09T22:18:01.763838Z",
     "shell.execute_reply.started": "2024-07-09T22:18:01.759708Z"
    }
   },
   "source": [
    "## 5. Combine a user-created table and the LSST query result table using `Pandas`\n",
    "\n",
    "The scientific scenario used for this demonstration is for a user-created `pandas` table\n",
    "of moving objects (the saved 1,000 Gaia MBA data is recycled as an example), whether LSST detected\n",
    "them during a specific night, and if so how the LSST data look like.\n",
    "The LSST query result will be also saved as a `pandas` table, and the join with the user-created\n",
    "table will be done using the `merge` function provided by `pandas`. \n",
    "\n",
    "### 5.1. Read the Gaia MBA table as a pandas `DataFrame`\n",
    "\n",
    "Follow the convention established in the sections above and name the table ut4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f9b4fea-1844-4725-b9aa-596e4cfdfa7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ut4 = pd.read_csv(my_home_dir+'gaia_mba.cat')\n",
    "ut4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abbe9138-473a-4c8e-9d8a-ddf9a42c5945",
   "metadata": {},
   "source": [
    "### 5.2 Retrieve the DP0.3 data and save the result as a `pandas` table\n",
    "\n",
    "As opposed to Section 4.2.2, the query in this section joins only two tables:\n",
    "the DP0.3 `MPCORB` table and the DP0.3 `diaSource` table.\n",
    "\n",
    "This query will return all LSST detections from the DP0.3 `diaSource` table\n",
    "for MBAs (2.0 < `a` < 3.25 au and `q` > 1.666 au) observed at 62001 < `midPointMjdTai` < 62002."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c77defa5-b937-4eed-9c31-2bc7d3bf899e",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"\n",
    "        SELECT mpc.fullDesignation, mpc.ssObjectId, dias.ssObjectId,\n",
    "        mpc.e, mpc.incl, mpc.q/(1-mpc.e) as a,\n",
    "        dias.midPointMjdTai, dias.mag, dias.band\n",
    "        FROM dp03_catalogs_10yr.MPCORB AS mpc\n",
    "        INNER JOIN dp03_catalogs_10yr.DiaSource as dias\n",
    "        ON dias.ssObjectId = mpc.ssObjectId\n",
    "        Where midPointMjdTai > 62001 AND midPointMjdTai < 62002\n",
    "        AND mpc.q/(1-mpc.e) > 2.0 AND mpc.q/(1-mpc.e) < 3.25\n",
    "        AND mpc.q > 1.666\n",
    "        \"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adaf666f-2a1d-4569-aa67-12104b8e6138",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = rsp_tap.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,
   "id": "758c216e-3013-4597-b739-69e1e2a4c561",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result().to_table().to_pandas()\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9db24f58-fb9e-4ddf-83c4-65239b7e8705",
   "metadata": {},
   "outputs": [],
   "source": [
    "uniqObjId = np.unique(results['ssObjectId'])\n",
    "print(f'{len(uniqObjId)} unique MBAs have LSST detections between MJD=62001 and 62002.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64348f7c-4c28-4166-a698-84fb480fa60d",
   "metadata": {},
   "source": [
    "### 5.3 Join two pandas `DataFrames` using the `merge` function\n",
    "\n",
    "`Pandas` can combine two DataFrames with different column names for a common key. \n",
    "The common key is the designation for objects, with the corresponding column name\n",
    "being `mpc_desig1` in the `ut4` table and `fullDesignation` in the `results` table.\n",
    "When using the 'inner' method, the merge operation returns only the rows that have\n",
    "matching values in both DataFrames. This type of join is known as an \"inner join\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "712c1013-f07e-42f8-8590-682d28be54c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_df = pd.merge(ut4, results, left_on='mpc_desig1', right_on='fullDesignation', how='inner')\n",
    "print(f'{len(merged_df)} matches between the ut4 ({len(ut4)} MBAs) and '\n",
    "      f'the results ({len(results)} diaSources) tables.')\n",
    "merged_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99467bd2-d9fb-493e-836d-b53b64b5af1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "uniqObj, cnts = np.unique(merged_df['denomination'], return_counts=True)\n",
    "most_detected = np.where(cnts == cnts.max())\n",
    "print(f'{len(uniqObj)} Gaia MBAs out of {len(ut4)} are detected with Rubin '\n",
    "      f'between MJD=62001 and 62002.')\n",
    "print(f'{len(most_detected[0])} of them have the maximum number '\n",
    "      f'of the LSST detections of {cnts.max()}.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e21472d4-2a88-4788-9e05-d9fa21d7450c",
   "metadata": {},
   "source": [
    "#### 5.3.1 Compare the orbital parameters from Gaia vs. the LSST for the matched MBAs\n",
    "\n",
    "Plot the orbital parameters such as eccentricity (`e`), inclination (`i`), and semi-major axis (`a`)\n",
    "from the Gaia and the LSST $g$-band for matched MBAs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6080b89-2280-4399-88dd-b5c0be0e90cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(16, 3))\n",
    "fig.suptitle('For matched MBAs between the Gaia and the LSST g-band')\n",
    "\n",
    "axs[0].scatter(merged_df['eccentricity'], merged_df['e'])\n",
    "axs[0].set_xlabel('e from Gaia')\n",
    "axs[0].set_ylabel('e from Rubin')\n",
    "\n",
    "axs[1].scatter(np.rad2deg(merged_df['inclination']), merged_df['incl'])\n",
    "axs[1].set_xlabel('i [deg] from Gaia')\n",
    "axs[1].set_ylabel('i [deg] from Rubin')\n",
    "\n",
    "axs[2].scatter(merged_df['semi_major_axis'], merged_df['a'])\n",
    "axs[2].set_xlabel('a [au] from Gaia')\n",
    "axs[2].set_ylabel('a [au] from Rubin')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83c07fad-a1c6-4231-9acd-3221ffdbc13f",
   "metadata": {},
   "source": [
    "Clean up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aad84648-eae2-4b13-b6bf-a4d41cd90d5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "del ut4, results, uniqObjId, merged_df, uniqObj, cnts, most_detected"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f880c4d0-6ffc-49e9-86f2-95f1b698e279",
   "metadata": {},
   "source": [
    "## 6. Excercises for the learner\n",
    "\n",
    "1. In Section 4.1, the Gaia database was queried for main-belt asteroids (MBAs) following the\n",
    "population definition used by the [JPL Horizons small body database query tool](https://ssd.jpl.nasa.gov/tools/sbdb_query.html).\n",
    "At that link, click on \"Limit by Orbit Class\" and then hover the cursor over any \"information\" symbol\n",
    "to see any definition. Use the definition of TransNeptunian Object (TNO) to rewrite the query in Section 4.1\n",
    "and create a table of the nine Gaia objects with >200 observations and the semi-major axis of a TNO.\n",
    "Clear the kernel and re-execute the notebook in full, using Gaia TNOs instead of MBAs; there\n",
    "should be 4 detected by LSST.\n",
    "\n",
    "2. (Advanced) Generate a table of DP0.3 objects of interest and save it as a file (e.g., the MPC designations of TNOs identified in step 1). Using Sections 3 and 4 as a guide, use the user-upload functionality and do a tripe-table join with the DP0.3 `SSObject` catalog (match on MPC designation)\n",
    "and `SSSource` catalog (match on `ssObjectId`) to obtain and plot, e.g., heliocentric distances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8501c9f-def1-4395-b0be-62c671201211",
   "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"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}