DP03_01_Introduction_to_DP03.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bc27f992-6cfd-4c2c-bb8b-5c8a7515f0db",
   "metadata": {},
   "source": [
    "# Introduction to the DP0.3\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "Contact author(s): Bob Abel, Douglas Tucker, and Melissa Graham<br>\n",
    "Last verified to run: 2024-12-17 <br>\n",
    "LSST Science Piplines version: Weekly 2024_50 <br>\n",
    "Container size: medium <br>\n",
    "Targeted learning level: beginner <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e57b4248-cf79-4945-956b-3186f3cc41c2",
   "metadata": {},
   "source": [
    "**Description:** An overview of the contents of the DP0.3 moving object catalogs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9fd4286-3bff-4345-affa-0d8e519b8d50",
   "metadata": {},
   "source": [
    "**Skills:** Use the TAP service and ADQL to access the DP0.3 tables."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e79e389-8abf-4dff-b503-55c0383084d1",
   "metadata": {},
   "source": [
    "**LSST Data Products:** TAP dp03_catalogs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07aae406-2b73-4995-a1c8-fc833a24cee2",
   "metadata": {},
   "source": [
    "**Packages:** `lsst.rsp.get_tap_service`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8740f32b-3d92-4c10-8d88-46e600a63a96",
   "metadata": {},
   "source": [
    "**Credit:**  The DP0.3 data set was generated by members of the Rubin Solar System Pipelines and Commissioning teams, with help from the LSST Solar System Science Collaboration, in particular: Pedro Bernardinelli, Jake Kurlander, Joachim Moeyens, Samuel Cornwall, Ari Heinze, Steph Merritt, Lynne Jones, Siegfried Eggl, Meg Schwamb, and Mario Jurić.\n",
    "\n",
    "Mario Jurić provided essential assistance in the initial stages of creating this notebook.  The notebook authors would also like to acknowledge help from Leanne Guy, Pedro Bernardinelli, Sarah Greenstreet, Megan Schwamb, Brian Rogers, Niall McElroy, and Jake Vanderplass, among others."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7b54235-6bab-4626-a603-13c15a3c93b6",
   "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",
   "id": "c84499ca-08c3-47a6-b1cc-6e289ec42fdf",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notebook demonstrates how to access the simulated Data Preview 0.3 (DP0.3) data set in the Rubin Science Platform.\n",
    "\n",
    "For the DP0.3 simulation, only moving objects were simulated, and only catalogs were created (there are no images).\n",
    "The DP0.3 simulation is *entirely independent of and separate from* the DP0.2 simulation.\n",
    "\n",
    "DP0.3 is a hybrid catalog that contains both real and simulated Solar System objects (asteroids, near-earth objects, Trojans, trans-Neptunian objects, comets, and even a simulated spaceship... but no major planets nor the Moon).\n",
    "See the <a href=\"https://dp0-3.lsst.io\">DP0.3 documentation</a> for more information about how the hybrid catalog was created.\n",
    "\n",
    "In Rubin Operations, these tables would be constantly changing, updated every day with the results of\n",
    "the previous night's observations.\n",
    "However, for DP0.3, both a static 1-year catalog and a static 10-year catalog have been simulated.\n",
    "\n",
    "> **Notice:** To re-iterate, for DP0.3, there are tables available for both the 1-year and 10-year catalogs.  For the remainder of this notebook, though, unless otherwise noted, we will just consider the tables for the 10-year catalog.  Users are encouraged to explore the tables for the 1-year catalog on their own.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fdd57af-83de-464a-863f-f0886a70d989",
   "metadata": {},
   "source": [
    "### 1.1. Package Imports\n",
    "\n",
    "Import general python packages and the Rubin TAP service utilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dab33fdd-8ff9-4f8f-b027-77ba5368de5d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colormaps\n",
    "import pandas as pd\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0484b1ba-eafd-465a-8691-5c1ba3bcc208",
   "metadata": {},
   "source": [
    "## 2. Create the Rubin SSO TAP Service client\n",
    "\n",
    "The DP0.3 data sets are available via the Table Access Protocol (TAP) service,\n",
    "and can be accessed with ADQL (Astronomical Data Query Language) statements.\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",
    "\n",
    "Get an instance of the SSO TAP service, and assert that it exists.\n",
    "\n",
    "> **Notice:** The DP0.3 TAP service is called `ssotap` (whereas `tap` is used for DP0.2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8daf49b-a897-479b-9524-f816bea54733",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"ssotap\")\n",
    "assert service is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d49c768a-e8c5-44a6-bde1-bfa456ac7f98",
   "metadata": {},
   "source": [
    "### 2.1. What `ssotap` schemas are available?\n",
    "\n",
    "Create an ADQL query to select all (\\*) available schemas, and use the TAP service to execute the query.\n",
    "\n",
    "Use `to_table` to convert to an `astropy` table, which will display the results in a user-friendly way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d62ab2ef-05df-4024-b514-dc02e4379041",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT * FROM tap_schema.schemas\"\n",
    "results = service.search(query).to_table()\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e56f59e6-28d7-49f9-9030-118f43bbf5ad",
   "metadata": {},
   "source": [
    "### 2.2 What DP0.3 tables are available?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16dbae17-6f9d-4356-941c-3e20498bc804",
   "metadata": {},
   "source": [
    "There are four tables for the 1-year and 10-year simulation: \n",
    "`MPCORB`, `SSObject`, `SSSource`, and `DiaSource`.\n",
    "\n",
    "Find descriptions of the tables and their schema, plus\n",
    "information and advice about accessing and querying the DP0.3 tables\n",
    "(including which columns are <a href=\"https://dp0-3.lsst.io/data-products-dp0-3/table-access-and-queries.html#unpopulated-columns\">currently unpopulated</a>),\n",
    "in the <a href=\"https://dp0-3.lsst.io/data-products-dp0-3/index.html#dp0-3-data-products-definition-document-dpdd\">DP0.3 data products definitions documentation</a>.\n",
    "\n",
    "This tutorial will explore the tables individually, with table joins to be\n",
    "demonstrated in future tutorials (see also <a href=\"https://dp0-3.lsst.io/data-products-dp0-3/table-access-and-queries.html#table-joins\">this advice on DP0.3 table joins</a>).\n",
    "\n",
    "Select all available tables and display information about them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79bbde6a-bd68-4de8-a02d-2160dfae6bc4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT * FROM tap_schema.tables \" \\\n",
    "        \"WHERE tap_schema.tables.schema_name LIKE 'dp03_catalogs%' \" \\\n",
    "        \"ORDER BY table_index ASC\"\n",
    "results = service.search(query).to_table()\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7302673-cab0-4003-a3dd-d22b1a77a97b",
   "metadata": {},
   "source": [
    "## 3. The `MPCORB` table\n",
    "\n",
    "During Rubin Operations, Solar System Processing will occur in the daytime, after a night of observing.\n",
    "It will link together the difference-image detections of moving objects and report discoveries\n",
    "to the <a href=\"https://minorplanetcenter.net\">Minor Planet Center</a> (MPC),\n",
    "as well as compute derived properties (magnitudes, phase-curve fits, coordinates in various systems).\n",
    "\n",
    "The MPC will calculate the orbital parameters and these results will be passed back to Rubin, and stored\n",
    "and made available to users as the `MPCORB` table \n",
    "(the other derived properties are stored in the other three tables explored below).\n",
    "The DP0.3 `MPCORB` table is a simulation of what this data product will be like after 10 years of LSST.\n",
    "\n",
    "> **Notice:** The MPC contains all reported moving objects in the Solar System, and is not limited to those detected by LSST. Thus, the `MPCORB` table will have more rows than the `SSObject` table.\n",
    "\n",
    "> **Notice:** For DP0.3, there was no fitting done by the MPC and the MPCORB table is the orbital elements used in the simulation (the `MPCORB` table is a truth table).\n",
    "\n",
    "For more information about Rubin's plans for Solar System Processing, see Section 3.2.2 of the \n",
    "<a href=\"https://docushare.lsstcorp.org/docushare/dsweb/Get/LSE-163/LSE-163_DataProductsDefinitionDocumentDPDD.pdf\">Data Products Definitions Document</a>.\n",
    "Note that there remain differences between Table 4 of the DPDD, which contain the anticipated schema \n",
    "for the moving object tables, and the DP0.3 table schemas."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd15dac9-45bf-4a0f-9f32-5d31bb48a4d5",
   "metadata": {},
   "source": [
    "### 3.1. Size\n",
    "\n",
    "Use the TAP service to count of the number of rows in the `MPCORB` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e965cf06-6e16-4117-a29b-97a6a8577fb2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT COUNT(*) FROM dp03_catalogs_10yr.MPCORB\")\n",
    "results.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dfddcd0-0a73-4f41-9324-2fc590bb0fbf",
   "metadata": {},
   "source": [
    "There are 14.5 million rows in the `MPCORB` table."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69152b41-0c0c-435c-884f-553ca4a682db",
   "metadata": {},
   "source": [
    "### 3.2. Columns\n",
    "\n",
    "Use the TAP service to query for the column information from `MPCORB`.\n",
    "\n",
    "Print the results as a `pandas` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d24359a4-6f1a-433b-84cf-b9184f4f673f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT column_name, datatype, description, \"\n",
    "                         \"       unit from TAP_SCHEMA.columns \"\n",
    "                         \"WHERE table_name = 'dp03_catalogs_10yr.MPCORB'\")\n",
    "results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06e52c7d-3e45-4f21-bddd-1b7e7d879e21",
   "metadata": {},
   "source": [
    "In some cases, the column descriptions cut off in the table above.\n",
    "\n",
    "Execute the following to see, for example, the full description for the `mpcDesignation` column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7743b75b-9c91-4756-892e-741e15bd3cef",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results['description'][10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50c4dcdf-9e70-4c56-aa2b-15c7bc6feb65",
   "metadata": {},
   "source": [
    "### 3.3. Retrieve a random subset\n",
    "\n",
    "To retrieve a random subset of rows, make use of the fact that `ssObjectId` is a randomly assigned 64-bit long unsigned integer. \n",
    "Since ADQL interprets a 64-bit long unsigned integer as a 63-bit _signed_ integer, \n",
    "these range from a very large negative integer value to a very large positive integer value.\n",
    "This will be fixed in the future so that all identifiers are positive numbers.\n",
    "\n",
    "> **Notice:** By using `ssObjectId`, the following query returns a random subset of `MPCORB` rows *that are associated with a row in the `SSObject` table*. In other words, this limits the query to only retrieve moving objects from the `MPCORB` table that have been detected by LSST.\n",
    "\n",
    "First, figure out the full range of `ssObjectId` values by using the ADQL `MIN` and `MAX` functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "daf15ffa-ac64-497b-8a08-2df29570fa36",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT min(ssObjectId), max(ssObjectId)\"\n",
    "                         \"FROM dp03_catalogs_10yr.MPCORB\")\n",
    "results.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "086af243-a5fe-4bc3-b7d9-57595cb0c637",
   "metadata": {},
   "source": [
    "Define a search range for `ssObjectId` that would return no more than 1% of all objects in `MPCORB`.  We will do this by estimating a new minimum `ssObjectId` that is 1% _below_ the maximum `ssObjectId` for the full range of `ssObjectId` values.\n",
    "\n",
    "> **Notice:** Since the _range_ of `ssObjectId`'s (-9223370383071521539 --> +9223370875126069107) is much larger than the number of _rows_ in the `MPCORB` table (14600302), we don't expect to get _exactly_ 1% of the rows from `MPCORB` via this method, but we should get approximately 1%, as long as the `ssObjectId` values are distributed reasonably uniformly over their large range.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad841f2b-80da-4ff4-9b2e-b56473cb6fee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "min_val = int(results[0].get('min1'))\n",
    "max_val = int(results[0].get('max2'))\n",
    "print('Full range: ', min_val, max_val)\n",
    "\n",
    "min_val = int(max_val - 0.01*(max_val-min_val))\n",
    "print('Search range: ', min_val, max_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ab5d7bb-6c76-4a06-9086-d6d92d458ce3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-27T22:39:44.565969Z",
     "iopub.status.busy": "2023-06-27T22:39:44.565076Z",
     "iopub.status.idle": "2023-06-27T22:39:44.570822Z",
     "shell.execute_reply": "2023-06-27T22:39:44.569984Z",
     "shell.execute_reply.started": "2023-06-27T22:39:44.565938Z"
    },
    "tags": []
   },
   "source": [
    "Execute the search, and retrieve all (\\*) columns from the `MPCORB` table."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfff1a2e-21ea-4ef4-9063-5f0bfac18b4f",
   "metadata": {},
   "source": [
    "Define query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94bde06e-a1eb-4c57-82ac-2a02a17d30e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"SELECT * FROM dp03_catalogs_10yr.MPCORB \n",
    "              WHERE ssObjectId BETWEEN 9038903462544093184 \n",
    "              AND 9223370875126069107\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bd96843-2309-4f06-9c12-47faad0305ea",
   "metadata": {},
   "source": [
    "Submit asynchronous query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1fc92d6-06ae-410b-8e95-6ad759d4c51e",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "results = job.fetch_result()\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e2fb6d0-dddf-4c35-aa1d-561f7fcf29a8",
   "metadata": {},
   "source": [
    "Store the results in `df` as a `pandas` dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f286f6a7-e5ec-4bf4-978a-26abbe71bf6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a94f02ec-84ae-4d9a-b81f-f335f27371fc",
   "metadata": {},
   "source": [
    "Display the resulting dataframe. Note that it is automatically truncated for both columns and rows \n",
    "(look for the \"...\" halfway down, and halfway across)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cdda7dc-977e-46f1-84ad-aad3b7ab087f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "528bd445-2510-4fbd-985e-f33c91f3e47f",
   "metadata": {},
   "source": [
    "As we see, 144,472 rows were returned, which -- _as expected_ -- is almost (but not exactly) 1% of the 14,462,388 rows in the `MPCORB` table.  \n",
    "\n",
    "> **Notice:** There are several columns that currently contain `NaN` (not a number) values.\n",
    "For the simulated DP0.3 data these columns might be replaced in the near future, \n",
    "and for real data releases there will not be all-`NaN` columns.  If desired, users can drop all-`NaN` columns with, e.g., `df.dropna(axis=1, how='all', inplace=True)`.\n",
    "However the better practice is to understand the columns and retrieve only what you are going to use."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaa28cea-3137-4a8a-a5a1-756a4943ddbb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-27T22:42:58.606568Z",
     "iopub.status.busy": "2023-06-27T22:42:58.605621Z",
     "iopub.status.idle": "2023-06-27T22:42:58.612817Z",
     "shell.execute_reply": "2023-06-27T22:42:58.611762Z",
     "shell.execute_reply.started": "2023-06-27T22:42:58.606523Z"
    },
    "tags": []
   },
   "source": [
    "**Optional:** use the `pandas` dataframe `info` function to learn more about the values in the retrieved columns.\n",
    "Uncomment the following code cell (remove the # symbol) and execute the cell to see the info on `df`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1752b220-1916-440b-8c3e-616c4e0623f7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4a18142-d555-4cba-a499-b80902ba4a2e",
   "metadata": {},
   "source": [
    "**Optional:** use the `pandas` fuction `describe` to display statistics for the \n",
    "numerical columns in the dataframe (count, mean, standard deviation, etc.)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cdf2982-1cc2-476c-a14b-9ba9d3561fb2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9721c47c-6302-4e0c-8d38-3825ff86266c",
   "metadata": {},
   "source": [
    "### 3.4. Plot histograms of selected columns\n",
    "\n",
    "Wikipedia provides a decent <a href=\"https://en.wikipedia.org/wiki/Orbital_elements\">beginner-level guide to orbital elements</a>.\n",
    "\n",
    "For a quick reference, distributions are shown below for five key orbital elements\n",
    "and the absolute $H$ magnitude (see Section 4.4 for a description of $H$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "478cba37-0c97-4132-83e5-6bde2083b6bc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 3, figsize=(10, 6), sharey=False)\n",
    "ax[0, 0].hist(df['e'], bins=100, log=True)\n",
    "ax[0, 0].set_xlabel('Eccentricity')\n",
    "ax[0, 0].set_ylabel('log(Number)')\n",
    "ax[0, 1].hist(df['incl'], bins=100, log=True)\n",
    "ax[0, 1].set_xlabel('Inclination (deg)')\n",
    "ax[0, 1].set_ylabel('log(Number)')\n",
    "ax[0, 2].hist(df['mpcH'], bins=100, log=True)\n",
    "ax[0, 2].set_xlabel('Absolute Magnitude, H (mag)')\n",
    "ax[0, 2].set_ylabel('log(Number)')\n",
    "ax[1, 0].hist(df['node'], bins=50)\n",
    "ax[1, 0].set_xlabel('Longitude of Ascending Node (deg)')\n",
    "ax[1, 0].set_ylabel('Number')\n",
    "ax[1, 0].set_ylim(0,3500)\n",
    "ax[1, 1].hist(df['peri'], bins=50)\n",
    "ax[1, 1].set_xlabel('Argument of Perihelion (deg)')\n",
    "ax[1, 1].set_ylabel('Number')\n",
    "ax[1, 1].set_ylim(0,3500)\n",
    "ax[1, 2].hist(df['q'], bins=100, log=True)\n",
    "ax[1, 2].set_xlabel('Perihelion Distance (au)')\n",
    "ax[1, 2].set_ylabel('log(Number)')\n",
    "fig.suptitle('Histograms for Key Orbital Elements')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "242b5cb8-0431-44e8-9aa1-2fb449845cc9",
   "metadata": {},
   "source": [
    "> **Figure 1**: Six histograms showing the log number of solar system objects in the `MPCORB` table as a function of their eccentricity, inclination, absolute H magnitude, longitude of the ascending node, argument of perihelion, and perihelion distance (from upper left to lower ight)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "463704a9-dc73-4f7c-abfd-43b2f9731539",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-17T15:47:28.771164Z",
     "iopub.status.busy": "2023-06-17T15:47:28.770077Z",
     "iopub.status.idle": "2023-06-17T15:47:28.774678Z",
     "shell.execute_reply": "2023-06-17T15:47:28.773860Z",
     "shell.execute_reply.started": "2023-06-17T15:47:28.771122Z"
    },
    "tags": []
   },
   "source": [
    "## 4. The `SSObject` table\n",
    "\n",
    "During Rubin Operations, Prompt Processing will occur during the night, detecting sources in \n",
    "difference images (`DiaSources`, see Section 6) and associating them into static-sky transients\n",
    "and variables (`DiaObjects`, not included in DP0.3).\n",
    "\n",
    "The Solar System Processing which occurs in the daytime, after a night of observing,\n",
    "links together the `DiaSources` for moving objects into `SSObjects`.\n",
    "Whereas the `MPCORB` table contains the orbital elements for these moving objects,\n",
    "the `SSObjects` contains the Rubin-measured properties such as phase curve fits and absolute magnitudes.\n",
    "\n",
    "> **Notice:** no artifacts or spurious difference-image sources have been injected into the DP0.3 catalogs.\n",
    "\n",
    "> **Notice:** although there are columns for them, no _u-_ or _Y_-band data were simulated for DP0.3.  Unless noted, we will ignore _u-_ or _Y_-band data for the remainder of this notebook.\n",
    "\n",
    "### 4.1. Size\n",
    "\n",
    "Use the ADQL count function to return the size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd943cfe-7c32-49da-b2eb-af12cf3b78b2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT COUNT(*) from dp03_catalogs_10yr.SSObject\")\n",
    "results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56f8dad2-5f9a-4393-960c-f652440ed3b1",
   "metadata": {},
   "source": [
    "The DP0.3 data set contains 4.4 million solar system objects detected by Rubin.\n",
    "\n",
    "This is less than the 14.5 million objects in the `MPCORB` catalog.\n",
    "It is left as an exercise for the learner in Section 7 to determine the characteristics of those \n",
    "objects from the `MPCORB` table are missing from the `SSObject` table."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb8901f4-f8d2-43a1-b2c4-cc020a31a30f",
   "metadata": {},
   "source": [
    "### 4.2. Columns\n",
    "\n",
    "**Option:** display a list of column names, data types, descriptions, and units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19816767-ea9b-4b23-a38d-c16a05450604",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT column_name, datatype, description, \"\n",
    "#                          \"unit from TAP_SCHEMA.columns \"\n",
    "#                          \"WHERE table_name = 'dp03_catalogs_10yr.SSObject'\")\n",
    "# results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4c7043b-4fc4-4c46-a819-9c3bdfed379a",
   "metadata": {},
   "source": [
    "### 4.3. Retrieve a random subset\n",
    "\n",
    "Use essentially the same query as was used for the `MPCORB` table, above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d27d03d6-91c5-4ec2-991f-8477a4e70fba",
   "metadata": {
    "tags": []
   },
   "source": [
    "Define query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4623d230-b4f5-49ba-9bdf-db283a6a19f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"SELECT * FROM dp03_catalogs_10yr.ssObject \n",
    "                WHERE ssObjectId BETWEEN 9038903462544093184 \n",
    "                AND 9223370875126069107\"\"\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "574d6c23-54f0-4731-ad8a-b940701143eb",
   "metadata": {},
   "source": [
    "Submit asynchronous query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9eed8a78-e6f9-4432-92aa-cd40672b253d",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "results = job.fetch_result()\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "387cd83a-a08f-4047-ae61-2d3a0633b995",
   "metadata": {},
   "source": [
    "Store the results in df as a pandas dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcb474d8-2081-419f-9c79-6e1cfae45164",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af55bc58-d60d-41fa-a6ef-829f81aa6052",
   "metadata": {},
   "source": [
    "Display the resulting dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26149b6f-2b7d-4cfd-b7b0-bfe8f06b4889",
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cd4af29-a072-4361-b65b-5979d4812783",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-06-27T21:51:31.545543Z",
     "iopub.status.busy": "2023-06-27T21:51:31.545077Z",
     "iopub.status.idle": "2023-06-27T21:51:31.549274Z",
     "shell.execute_reply": "2023-06-27T21:51:31.548724Z",
     "shell.execute_reply.started": "2023-06-27T21:51:31.545521Z"
    },
    "tags": []
   },
   "source": [
    "There are 44611 rows."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c1b56c6-5874-4f5c-a7cd-d0a6daa69a84",
   "metadata": {},
   "source": [
    "**Options:** use the `info` or `describe` functions on the `pandas` dataframe to learn more about the retrieved results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeeb504b-3ae8-4483-ac8a-809ab113e213",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8571bad-0431-4040-a083-cab4364a7927",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5c93746-8820-4dc1-8102-4d67b393b542",
   "metadata": {},
   "source": [
    "### 4.4. Plot a color-color diagram\n",
    "\n",
    "In the displayed dataframe above, it appears that for many `SSObjects` the phase-curve fits to derive \n",
    "absolute magnitudes were not successful, as the `<band>_H` (where `<band>` is a band u, g, r, i, z, or y) are `NaN`.  (Recall in  particular that no _u-_ or _y_-band data were simulated for DP0.3.)\n",
    "\n",
    "Before calculating and plotting the colors, drop all of the rows for which the phase-curve fits were not successful for g, r, i, and/or z bands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "937fb5ee-ef7b-4ef3-98ea-8d7237fe2a95",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df.dropna(subset=['g_H', 'r_H', 'i_H', 'z_H'], inplace=True)\n",
    "df.reset_index(inplace=True)\n",
    "print('Number of rows after dropping rows: %d' % len(df))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00180672-d7d6-40df-9d85-c51bb298dce6",
   "metadata": {},
   "source": [
    "For Solar System objects, absolute magnitudes are defined to be for an object 1 au from the Sun and 1 au \n",
    "from the observer, and at a phase angle (the angle Sun-object-Earth) of 0 degrees.\n",
    "\n",
    "Absolute magnitudes are derived by correcting for distance, fitting a function to the relationship between \n",
    "absolute magnitude and phase, and evaluating the function at a phase of 0 deg.\n",
    "\n",
    "The process for fitting phase curves will be covered in another tutorial.\n",
    "\n",
    "Calculate colors in the Rubin filters for the `SSObjects` that have absolute magnitudes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16eeddc9-f0dc-470c-ab96-252b4cfc66f6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df['gr'] = df['g_H'] - df['r_H']\n",
    "df['ri'] = df['r_H'] - df['i_H']\n",
    "df['iz'] = df['i_H'] - df['z_H']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8296082b-a528-4abd-92ce-52707da8ccc4",
   "metadata": {},
   "source": [
    "Plot color-color diagrams as 2-dimensional histograms (heatmaps)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "649578b9-496c-4703-8938-b3084d37cfd2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(8, 3))\n",
    "\n",
    "h = ax[0].hist2d(df['gr'], df['ri'],\n",
    "                 bins=(np.arange(-0.5, 1.5, 0.01),\n",
    "                       np.arange(-0.5, 0.75, 0.01)),\n",
    "                 norm='log')\n",
    "ax[0].set_xlabel('g - r')\n",
    "ax[0].set_ylabel('r - i')\n",
    "ax[0].grid()\n",
    "plt.colorbar(h[3])\n",
    "\n",
    "h = ax[1].hist2d(df['ri'], df['iz'],\n",
    "                 bins=(np.arange(-0.5, 0.75, 0.01),\n",
    "                       np.arange(-0.75, 0.5, 0.01)),\n",
    "                 norm='log')\n",
    "ax[1].set_xlabel('r - i')\n",
    "ax[1].set_ylabel('i - z')\n",
    "ax[1].grid()\n",
    "plt.colorbar(h[3])\n",
    "\n",
    "fig.suptitle('Color-Color Plots for the SSObject Catalog')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "653f251d-a172-4cda-b3c5-0aeecf28d8e6",
   "metadata": {},
   "source": [
    "> **Figure 2:** Two panels showing color-color plots for the retrieved subset of the `SSObject` catalog. At left, the $r-i$ vs. $g-r$ color and at right, the $i-z$ versus $r-i$ color. In both panels, the color distribution is bimodal.\n",
    "\n",
    "As shown in the figure above, there are two colors used in the simulation - but this is not the case for real Solar System objects.\n",
    "These plots will look very different in the future, when they are made with real Rubin data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef9ba117-4a6c-4a4c-9b82-5b74a9c40001",
   "metadata": {},
   "source": [
    "## 5. The `SSSource` table\n",
    "\n",
    "As described above, Solar System Processing links together the `DiaSources` (detections in the\n",
    "individual difference images) from moving objects into `SSObjects`. \n",
    "\n",
    "The `SSSource` table contains the 2-d (sky) coordinates and 3-d distances and velocities \n",
    "for every `SSObject` at the time of every LSST observation of that `SSObject`.\n",
    "The `SSSource` and `DiaSource` tables are 1:1, as they each contain data *per observation*,\n",
    "whereas `SSObject` contains data *per object*."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dc55526-5aa8-47cd-8d2a-75f02b8b48b7",
   "metadata": {},
   "source": [
    "### 5.1. Size\n",
    "\n",
    "It can take up to a minute to retrieve the size of the `SSSource` catalog."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "000450cd-dcbf-4001-bd35-7d925d82468e",
   "metadata": {},
   "source": [
    "Define query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "840d8a1d-2777-47b4-a050-8d8d627f58d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT COUNT(*) from dp03_catalogs_10yr.SSSource\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8956663-94b2-411a-80ea-56bf54351bdb",
   "metadata": {},
   "source": [
    "Submit asynchonous query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf8f6fe1-0d50-4134-81b6-82f64f3559f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "results = job.fetch_result()\n",
    "assert job.phase == 'COMPLETED'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7e104da-e779-46e2-8c92-ef71f58cce1f",
   "metadata": {},
   "source": [
    "Display results as a pandas data frame. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4900874-b7b1-49a5-af5d-a99de6affc70",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47b9a7c5-4f38-42f3-96b8-732c72010b81",
   "metadata": {},
   "source": [
    "**_This table contains over 653 million sources!_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05883359-a838-40ce-8106-494262b5435a",
   "metadata": {},
   "source": [
    "### 5.2. Columns\n",
    "\n",
    "**Option:** print the column information for the `SSSource` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "450c92ce-cc7e-4c6e-a6bd-a423d713c8ee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT column_name, datatype, description, \"\n",
    "#                          \"unit from TAP_SCHEMA.columns \"\n",
    "#                          \"WHERE table_name = 'dp03_catalogs_10yr.SSSource'\")\n",
    "# results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55f7fd6f-a5c8-4e58-8781-5d7f14ebf530",
   "metadata": {},
   "source": [
    "### 5.3. Retrieve data for one `SSObject`\n",
    "\n",
    "It is possible to obtain the `SSSource` data for a set of `SSObjects`.\n",
    "For example, to retrieve all `SSSources` for the `SSObjects` retrieved in Section 4,\n",
    "use a query such as: \n",
    "\n",
    "```\n",
    "SELECT * FROM dp03_catalogs_10yr.SSSource\n",
    "WHERE ssObjectId BETWEEN 9038903462544093184 AND 9223370875126069107\n",
    "```\n",
    "\n",
    "However, the better way to demonstrate the data in the `SSSource` table is to look at just one `SSObject`,\n",
    "and the one with an `ssObjectId` = `6793512588511170680` is a fun choice.\n",
    "\n",
    "Retrieve the heliocentric (sun-centered) and topocentric (observatory-centered) X and Y coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9f81967-a036-4f71-af59-2e2254e37977",
   "metadata": {},
   "source": [
    "Define the query and run the asynchronous query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0c0b207-7c44-46d7-b776-221f357e6467",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"\"\"SELECT heliocentricX, heliocentricY, \n",
    "              topocentricX, topocentricY, ssObjectId \n",
    "              FROM dp03_catalogs_10yr.SSSource \n",
    "              WHERE ssObjectId = 6793512588511170680\"\"\"\n",
    "\n",
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "results = job.fetch_result()\n",
    "assert job.phase == 'COMPLETED'\n",
    "\n",
    "df_xy = results.to_table().to_pandas()\n",
    "print('Retrieved ', len(df_xy), ' rows.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1a6bd96-c573-4b43-9310-cc3cc9411b42",
   "metadata": {},
   "source": [
    "**Options** to display the table in full or use the `info` or `describe` functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "defd1bb3-9065-4ca8-9171-0b2ce125a462",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df_xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77771ca4-ca51-4617-a9ce-8939595272b7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df_xy.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3daf9b96-9444-49b6-a702-aded212fa5ca",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df_xy.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11c32da4-7632-4ca0-8e51-c3ae130bde77",
   "metadata": {},
   "source": [
    "### 5.3. Plot the locations of one `SSObject`\n",
    "\n",
    "Plot the locations of the selected `SSObject` at the time of every \n",
    "LSST observation using the X and Y heliocentric (Sun-centered; orange star) \n",
    "and topocentric (observatory-centered; blue circle) coordinates.\n",
    "This can be considered a projection of the orbit into the plane of the Solar System.\n",
    "\n",
    "Notice how the points are not regularly spaced.\n",
    "This is because there is one point per LSST observation of the object,\n",
    "and in some years it receives more or fewer observations.\n",
    "\n",
    "Notice how the points appear in an ellipse around the Sun with heliocentric coordinates (left).\n",
    "This is because the selected object is in the main asteroid belt and close enough to the Sun \n",
    "to complete at least on orbit during the 10-year LSST survey.\n",
    "Had the selected object been in the outer Solar System, or were this tutorial using\n",
    "the 1-year data set, the plot below would show an arc instead of an ellipse.\n",
    "\n",
    "The plot of the topocentric coordinates (right) does not appear elliptical because\n",
    "the motion of the Earth with respect to the object over the 10 years of the LSST\n",
    "is imprinted into the data. \n",
    "For the topocentric coordinates, the Earth's rotation also contributes, \n",
    "but it is a much smaller effect on the scale of these plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0464b3fa-299a-4d27-8f47-f19cd21062b6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax[0].grid()\n",
    "ax[0].plot(df_xy['heliocentricX'], df_xy['heliocentricY'],\n",
    "           'o', ms=4, mew=0, color='black')\n",
    "ax[0].plot(0, 0, '*', ms=15, color='darkorange')\n",
    "ax[0].set_xlabel('heliocentric X (au)')\n",
    "ax[0].set_ylabel('heliocentric Y (au)')\n",
    "ax[0].set_title('Heliocentric Coordinates')\n",
    "ax[1].grid()\n",
    "ax[1].plot(df_xy['topocentricX'], df_xy['topocentricY'],\n",
    "           'o', ms=4, mew=0, color='black')\n",
    "ax[1].plot(0, 0, 'o', ms=15, color='dodgerblue')\n",
    "ax[1].set_xlabel('topocentric X (au)')\n",
    "ax[1].set_ylabel('topocentric Y (au)')\n",
    "ax[1].set_title('Topocentric Coordinates')\n",
    "fig.suptitle('XY Path of 6793512588511170680')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79a48c31-d2ad-41ed-916a-aa46f60c1894",
   "metadata": {},
   "source": [
    "> **Figure 3:** The two-dimensional projection of the location of the selected solar system object into the XY plan in heliocentric (sun-centered; left) and topocentric (earth-centered; right) coordinates. At left, the Sun is represented by an orange star, and at right, the Earth is represented by a blue circle. In both plots the object is represented by small black circles."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c568937-4f20-489c-b113-711f32cfad34",
   "metadata": {},
   "source": [
    "## 6. The `DiaSource` catalog\n",
    "\n",
    "Last but definitely not least, the `DiaSource` table - which is actually the first to be generated\n",
    "by the Prompt Processing pipeline.\n",
    "This table will contain measurements for all sources detected with a signal-to-noise ratio of at least 5 \n",
    "in a difference image, including moving and static-sky sources.\n",
    "\n",
    "However, for simulated DP0.3 the `DiaSource` table contains *only moving objects*, and no static sky time-domain objects (and no detector artifacts).\n",
    "\n",
    "For DP0.3 the simulated `DiaSource` table contains only a subset of the columns that the\n",
    "real `DiaSource` table will have; see Table 1 of the Rubin \n",
    "<a href=\"https://docushare.lsstcorp.org/docushare/dsweb/Get/LSE-163/LSE-163_DataProductsDefinitionDocumentDPDD.pdf\">Data Products Definitions Document</a>.\n",
    "Furthermore, the `DiaSource` table contains a few extra truth columns, such as `raTrue`, `decTrue`, `magTrue`.\n",
    "\n",
    "For DP0.3, no photometric variability due to, e.g., the rotation of non-spherical bodies, collisions,\n",
    "or outgassing events were simulated. \n",
    "All evolution in apparent magnitudes with time are due to the phase curve, which is explored\n",
    "in another tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c11588b-6129-4b89-acc4-a14f13c09053",
   "metadata": {},
   "source": [
    "### 6.1. Size\n",
    "\n",
    "**Option** to retrieve the size of the `DiaSource` table (it is 1:1 with the `SSSource` table, so the same size)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97cc39ef-af55-4244-939f-13a79a0f944d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT COUNT(*) from dp03_catalogs_10yr.DiaSource\")\n",
    "# results.to_table().to_pandas()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abdb5907-467b-4c2d-9797-0dbf7659a796",
   "metadata": {},
   "source": [
    "### 6.2. Columns\n",
    "\n",
    "**Option** to print the column information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad2f5716-4310-4393-89fc-0299b176a169",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT column_name, datatype, description, \"\n",
    "#                          \"unit from TAP_SCHEMA.columns \"\n",
    "#                          \"WHERE table_name = 'dp03_catalogs_10yr.DiaSource'\")\n",
    "# results.to_table().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebb6fe98-f238-4c46-8964-090eff82bec1",
   "metadata": {},
   "source": [
    "### 6.3. Retrieve data for one `SSObject`\n",
    "\n",
    "Similar to what we did for `SSSource` data in Section 5, it is possible to obtain the `DiaSource` data for a set of `SSObjects`.  For example, to retrieve all `DiaSources` for the `SSObjects` retrieved in Section 4,\n",
    "use a query such as: \n",
    "\n",
    "```\n",
    "SELECT * FROM dp03_catalogs_10yr.DiaSource\n",
    "WHERE ssObjectId BETWEEN 9038903462544093184 AND 9223370875126069107\n",
    "```\n",
    "\n",
    "However, as we saw for the `SSSource` data in Section 5, the better way to demonstrate the data in the `DiaSource` table is to look at just one `SSObject`.\n",
    "Here, as in Section 5, we will use the one with an `ssObjectId` = `6793512588511170680`.\n",
    "\n",
    "Retrieve the right ascension, declination, and the <a href=\"https://en.wikipedia.org/wiki/International_Atomic_Time\">TAI</a> at the exposure's midpoint (MJD) of the observation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10af5ce9-89b6-4d6b-9a52-55c3e1c4dfa2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"\"\"SELECT ra, dec, midpointMjdTai \n",
    "              FROM dp03_catalogs_10yr.DiaSource \n",
    "              WHERE ssObjectId = 6793512588511170680\"\"\"\n",
    "\n",
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)\n",
    "results = job.fetch_result()\n",
    "assert job.phase == 'COMPLETED'\n",
    "\n",
    "df = results.to_table().to_pandas()\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8985789e-a504-454d-b920-45c474e93747",
   "metadata": {},
   "source": [
    "### 6.4. Plot time-domain data for one `SSObject`\n",
    "\n",
    "Here, we plot the equatorial coordinates for all observations, colored by the MJD of the observation.\n",
    "\n",
    "This plot shows the location of the object on the sky, as seen from Earth.\n",
    "Recall from the plots in Section 5 that this object is only about 2 au from the Sun,\n",
    "and so is at the inner edge of the Main Asteroid Belt.\n",
    "Parallax from the Earth's orbit contributes to the appearance of this plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "064324af-f6fa-46ec-a425-036ce96af31c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cmap = colormaps['viridis']\n",
    "fig, ax = plt.subplots(1, 1, figsize=(12, 3))\n",
    "im = ax.scatter(df['ra'], df['dec'], c=df['midpointMjdTai'], cmap=cmap)\n",
    "ax.set_xlabel('ra (deg)')\n",
    "ax.set_ylabel('dec (deg)')\n",
    "fig.suptitle(\"Path of one SSObject in equatorial coordinates over 10 years\")\n",
    "fig.subplots_adjust(right=0.8)\n",
    "cbar_ax = fig.add_axes([0.85, 0.15, 0.02, 0.7])\n",
    "fig.colorbar(im, cax=cbar_ax, label='TAI')\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e75a0bd3-022a-4277-934d-23e7de048602",
   "metadata": {},
   "source": [
    "> **Figure 4:** The declination (dec) versus right ascencion (ra) in decimal degrees of the difference-image detections of the selected solar system object. The points are colored by the date (TAI) of the observation to illustrate the object's path in the sky over time.\n",
    "\n",
    "Use of the magnitude and filter data retrieved from the `DiaSource` table is left as\n",
    "an exercise for the learner in Section 7."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68e9e475-96e2-4f47-8164-37d3ddbf67a5",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 7. Exercises for the learner\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27380304-2669-4e14-ba81-f57fd072ebd9",
   "metadata": {},
   "source": [
    "**1. Calculate the semi-major axis for a subset of objects.**\n",
    "\n",
    "The orbital element of semi-major axis ($a$) is not pre-computed in the `MPCORB` table\n",
    "because it can be derived from the orbit’s ellipticiy ($e$) and perihelion distance ($q$)\n",
    "with $a = q/(1-e)$.\n",
    "\n",
    "Copy the query from Section 3.3 and alter it to retrieve only $e$ and $q$ for a subset of objects.\n",
    "Then add a column to the results table for semi-major axis and plot eccentricity \n",
    "versus semi-major axis.\n",
    "\n",
    "Hints:\n",
    " 1. Restrict the query to objects with eccentricities between 0 and 1.\n",
    " 2. `df['a'] = df['q'] / (1.0 - df['e'])`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6822d11e-b160-472b-972a-cb5d16f9be5c",
   "metadata": {},
   "source": [
    "**2. How does an SSObject's magnitude change with time?**\n",
    "\n",
    "In Section 6.4, the magnitude and filter for all detections of the `SSObject` in the\n",
    "`DiaSource` table were not retrieved; alter the query to retrieve `mag` and `band`,\n",
    "and plot the magnitude as a function of time for the filter of your choice\n",
    "(recall: there will be no *u*- or *y*-band magnitudes).\n",
    "Note that the only simulated cause of photometric changes in DP0.3 objects is the distance from Earth and the phase curve.\n",
    "Correcting for distance, and the applications phase curves, are covered in another tutorial."
   ]
  }
 ],
 "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": 5
}