api_from_noirlab_short.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "565a9a28-2998-4736-94c3-5b236bcddaa3",
   "metadata": {},
   "source": [
    "# API for the RSP from NOIRLab's Astro DataLab (Short Version)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76b7cf8d-90fb-4b95-8457-045853740421",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=190 style=\"padding: 10px\">\n",
    "<br>\n",
    "\n",
    "**Contact authors:** Leanne Guy, Melissa Graham <br>\n",
    "**Last verified:** Fri Dec 1 2023 <br>\n",
    "**Recommended platform:** NOIRLab's Astro DataLab  (https://datalab.noirlab.edu/)<br>\n",
    "**Kernel for recommended platform:** Python3 <br>\n",
    "**Rubin data release:** Data Preview 0.2 (DP0.2) <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b26d6096-080f-4d9a-a3c5-60f4a7bd4f15",
   "metadata": {},
   "source": [
    "**This is a short, advanced version of the `api_from_noirlab.ipynb` notebook.**\n",
    "\n",
    "**Description:**\n",
    "This tutorial demonstrates how to access LSST-like data via the RSP from the NOIRLab Astro DataLab.\n",
    "The scenario starts with data that would be in an LSST alert packet, then does a remote query\n",
    "for additional data.\n",
    "\n",
    "**Credit:** Sections 1 and 2 are based on the <a href=\"https://dp0-2.lsst.io/data-access-analysis-tools/api-intro.html\">Introduction to the RSP API Aspect</a> webpage which had major contributions from Douglas Tucker.\n",
    "\n",
    "**Requirements:** Accounts in the NOIRLab Astro DataLab and the Rubin Science Platform (RSP) at https://data.lsst.cloud/. Only individuals with <a href=\"https://docushare.lsst.org/docushare/dsweb/Get/RDO-013\">Rubin data rights</a> may have an RSP account. See the <a href=\"https://dp0-2.lsst.io/dp0-delegate-resources/index.html#delegate-homepage-getting-started-checklist\">getting started with DP0.2 checklist</a> for instructions about how to request an RSP account."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45f4d0be-1030-4314-90d9-eb6f53a14854",
   "metadata": {},
   "source": [
    "## 1. Set up RSP token in DataLab\n",
    "\n",
    "Start a JupyterLab session at NOIRLab's Astro DataLab (https://datalab.noirlab.edu/).\n",
    "\n",
    "**THE TOKEN IS A PASSWORD.**\n",
    "Keep it secret. Keep it safe.\n",
    "\n",
    "Follow steps 1 through 5 of <a href=\"https://nb.lsst.io/environment/tokens.html#using-a-token-outside-the-science-platform\">these instructions to obtain a token for an RSP account</a>.\n",
    "Be sure to select the box for \"read:image\".\n",
    "\n",
    "Create a hidden file `~/.rsp-tap.token`. \n",
    "Use `chmod 600 .rsp-tap.token` to ensure user read/write permissions only.\n",
    "Then put the token in the file."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eeee261-191f-4165-ae61-bb9b3c885389",
   "metadata": {},
   "source": [
    "## 2. Set up RSP TAP service\n",
    "\n",
    "### 2.1. Import packages\n",
    "\n",
    "Import <a href=\"https://pyvo.readthedocs.io/en/latest/\">PyVO</a>,\n",
    "<a href=\"https://docs.python.org/3/library/urllib.htmlhttps://docs.python.org/3/library/urllib.html\">urllib</a>,\n",
    "and a variety of other common packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef3fc107-5ee9-4ac4-9d70-d62ce79b2b21",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pyvo\n",
    "from pyvo.dal.adhoc import DatalinkResults\n",
    "from pyvo.dal.adhoc import SodaQuery\n",
    "import os, getpass\n",
    "import pandas\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from urllib.request import urlretrieve\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.io import fits\n",
    "from astropy.wcs import WCS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c8361fb-1850-4c3a-9aed-a8879898d9c9",
   "metadata": {},
   "source": [
    "### 2.2. Set up RSP TAP credentials\n",
    "\n",
    "Define the filename of the token file created in Section 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d729bb7f-4e9b-42ac-83dc-f75360ceea6e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_username = getpass.getuser()\n",
    "token_filename = os.getenv('HOME')+'/.rsp-tap.token'\n",
    "assert os.path.exists(token_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23066f32-364f-477a-8d27-45efaebeff12",
   "metadata": {},
   "source": [
    "Get the token from the token file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad0bd516-4273-4e27-8e85-8467e9d9987c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "with open(token_filename, 'r') as f:\n",
    "    token = f.readline()\n",
    "assert token is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd3f9d6e-a551-4c96-84d4-6eeedf9a8548",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-30T07:04:39.283792Z",
     "iopub.status.busy": "2023-11-30T07:04:39.282998Z",
     "iopub.status.idle": "2023-11-30T07:04:39.288790Z",
     "shell.execute_reply": "2023-11-30T07:04:39.287355Z",
     "shell.execute_reply.started": "2023-11-30T07:04:39.283735Z"
    },
    "tags": []
   },
   "source": [
    "**Do not** `print(token)`.\n",
    "\n",
    "Establish the RSP TAP service (`rsp_tap`) using `pyvo`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "807bc17f-bc8b-429b-b42f-ea60eee466a4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cred = pyvo.auth.CredentialStore()\n",
    "cred.set_password(\"x-oauth-basic\", token)\n",
    "credential = cred.get(\"ivo://ivoa.net/sso#BasicAA\")\n",
    "rsp_tap_url = 'https://data.lsst.cloud/api/tap'\n",
    "rsp_tap = pyvo.dal.TAPService(rsp_tap_url, credential)\n",
    "assert rsp_tap is not None\n",
    "assert rsp_tap.baseurl == rsp_tap_url"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b111c155-8a2e-469d-b503-51dd5356202f",
   "metadata": {},
   "source": [
    "### 2.3. Optional DP0.2 test query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8b3f7ef-52be-47b6-9930-34660ddae2ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# query = \"SELECT * FROM tap_schema.schemas\"\n",
    "# results = rsp_tap.run_sync(query).to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87dc33f0-139f-4bed-8492-00f840fb99de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d173c585-f094-474b-a94f-5bad1363ea2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# del results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e02dff7-8f29-4708-b8f1-666a9eb35ae6",
   "metadata": {},
   "source": [
    "### 2.4. Optional DP0.3 test query"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28288ca1-8473-486d-abe3-9565b9251ea7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# rsp_tap_url_sso = 'https://data.lsst.cloud/api/ssotap'\n",
    "# rsp_tap_sso = pyvo.dal.TAPService(rsp_tap_url_sso, credential)\n",
    "# assert rsp_tap_sso is not None\n",
    "# assert rsp_tap_sso.baseurl == rsp_tap_url_sso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bc6e95a-1df7-42d5-b4ad-3451ba952ecd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# query_sso = \"SELECT * FROM tap_schema.schemas\"\n",
    "# results_sso = rsp_tap_sso.run_sync(query_sso).to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a62b56a3-c1a2-4cc0-a056-9b40f14906ab",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results_sso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2e856a7-8789-4a0f-9645-2d5ed51d969a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# del results_sso, rsp_tap_url_sso, rsp_tap_sso, query_sso"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4f60b2b-7b17-4fa3-8d28-c1fb4b363836",
   "metadata": {},
   "source": [
    "## 3. Obtain potential host galaxy information\n",
    "\n",
    "Pretend the `diaObjectId`, and the `objectId` for the three nearest\n",
    "extended objects from the latest data release, are known for five\n",
    "transients of interest.\n",
    "\n",
    "**In the future, this data would be in an LSST alert packet.**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eea98e6-cdc4-471d-a49b-de4c3c67badb",
   "metadata": {},
   "source": [
    "### 3.1. Simulate nearby-object data from LSST alert packets\n",
    "\n",
    "Create a pandas dataframe containing:\n",
    "\n",
    " * `diaObjectId` : identifier in the DP0.2 `DiaObject` table for the candidate supernova\n",
    " * `diaObject_coord` : coordinates [RA, Dec] in decimal degrees for the candidate supernova\n",
    " * `gals_objId` : identifier in the DP0.2 `Object` table for the three nearest galaxies\n",
    " * `gals_2Ddist` : the 2D sky distance in arcseconds between `DiaObject` and nearby galaxy's center\n",
    "\n",
    "for five potential Type Ia supernovae in DP0.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5befacdb-bd6b-40b4-a613-b7b6578dcacf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "d = {'diaObjectId' : [1568026726510894110, 1569909090417642499, 1653700672547196623, \n",
    "                      1734140943235288573, 1825796232526695593],\n",
    "     'diaObject_coord' : [[63.6025914, -38.634654],\n",
    "                          [69.9257038, -38.1424959],\n",
    "                          [70.8210894, -35.9915118],\n",
    "                          [52.5432991, -34.9028848],\n",
    "                          [71.7356252, -34.2191764]],\n",
    "     'gals_objId' : [[1568026726510919266, 1568026726510919261, 1568026726510919497],\n",
    "                     [1569425305301455007, 1569425305301455003, 1569425305301455014],\n",
    "                     [1653700672547231391, 1653700672547231402, 1653700672547231397],\n",
    "                     [1734140943235326493, 1734140943235293084, 1734140943235326492],\n",
    "                     [1739084347513803559, 1739084347513803574, 1739084347513803571]],\n",
    "     'gals_2Ddist' : [[3.15, 4.62, 5.08],\n",
    "                      [0.02, 3.08, 3.98],\n",
    "                      [2.13, 2.58, 4.34],\n",
    "                      [4.7, 5.98, 6.11],\n",
    "                      [0.03, 4.64, 5.93]]}\n",
    "df = pandas.DataFrame(data=d)\n",
    "del d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51f27d88-45ca-4d81-a578-247d932688a0",
   "metadata": {},
   "source": [
    "Option to display the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffa2f230-34b3-4bb6-9b94-67096908fdad",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15038959-a554-4ff6-ae7a-056884d8f423",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-30T18:50:00.207579Z",
     "iopub.status.busy": "2023-11-30T18:50:00.206900Z",
     "iopub.status.idle": "2023-11-30T18:50:00.215228Z",
     "shell.execute_reply": "2023-11-30T18:50:00.213445Z",
     "shell.execute_reply.started": "2023-11-30T18:50:00.207525Z"
    }
   },
   "source": [
    "### 3.2. Retrieve object data from the RSP's DP0.2 catalog\n",
    "\n",
    "Choose to continue with the first `diaObject` in the dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e955c1c9-80bc-4755-a7bd-033255387693",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "diao_index = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "347fe5ce-5b76-48d8-890a-7473abe5dbb4",
   "metadata": {},
   "source": [
    "Create `list_objId`, a string containing a comma-separated list of the three `objectId` for\n",
    "the three nearest galaxies to the selected `DiaObject`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b71be88-3503-4927-88e4-5c3c5210885e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "temp = np.asarray(df['gals_objId'][diao_index], dtype='int')\n",
    "list_objId = \"(\" + ','.join(['%20i' % num for num in temp]) + \")\"\n",
    "del temp\n",
    "print(list_objId)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8a2e56d-c7db-42a0-b7c5-360326bfa520",
   "metadata": {},
   "source": [
    "Create a query to retreive object astrometry, shape, size, and photometry measurements from the DP0.2 `Object` catalog.\n",
    "See the <a href=\"https://dp0-2.lsst.io/data-products-dp0-2/index.html#dp0-2-data-products-definition-document-dpdd\">DP0.2 DPDD</a> and <a href=\"https://dm.lsst.org/sdm_schemas/browser/dp02.html\">DP0.2 schema browser</a> for more information about the columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c10377d8-4e96-4e7d-8447-d7eb691e8fd6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT objectId, coord_ra, coord_dec, refExtendedness, \"\\\n",
    "        \"shape_xx, shape_xy, shape_yy, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(g_cModelFlux) AS g_cModelMag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(r_cModelFlux) AS r_cModelMag, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(i_cModelFlux) AS i_cModelMag \"\\\n",
    "        \"FROM dp02_dc2_catalogs.Object \"\\\n",
    "        \"WHERE objectId IN \"+list_objId\n",
    "del list_objId"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db5e74b8-6737-406e-9507-c26a440bfe87",
   "metadata": {},
   "source": [
    "Execute the query using the `rsp_tap` service, and store the results in `galaxies` as a table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7bceb46-4cd1-45a9-b32e-0308a431cfb4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "galaxies = rsp_tap.search(query).to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53bd1f44-ef10-461d-960a-334901a22dd6",
   "metadata": {},
   "source": [
    "Option to view the remotely retrieved data table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12de5c03-9bc2-4aab-8a56-7bb75381ee30",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9641df26-5649-4e86-aa65-9957a953ef18",
   "metadata": {},
   "source": [
    "### 3.3. Calculate additional galaxy properties\n",
    "\n",
    "#### 3.3.1. Galaxy colors\n",
    "\n",
    "Calculate the colors of the galaxies in $g-r$ and $r-i$ magnitude, and add them to the `galaxies` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc57f2c7-f76d-4034-a538-24548624cf97",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "galaxies['gr_clr'] = galaxies['g_cModelMag'] - galaxies['r_cModelMag']\n",
    "galaxies['ri_clr'] = galaxies['r_cModelMag'] - galaxies['i_cModelMag']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a68a331-946a-47aa-ac7d-8d5ad0f2ebb1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78238ee8-9fe4-47cf-adf8-ab664b29541e",
   "metadata": {},
   "source": [
    "#### 3.3.2. Separation in elliptical radii\n",
    "\n",
    "**In the future, the LSST alert packet will contain separation distances in elliptical radii\n",
    "that are based on the second moments of the galaxy's luminosity profile.**\n",
    "\n",
    "Calculate them and add them to the `galaxies` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b012d44-6771-4068-8738-2b2772e88efe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "galaxies['ell_rad'] = np.zeros(3, dtype='float')\n",
    "galaxies['2Ddist'] = np.zeros(3, dtype='float')\n",
    "\n",
    "snra = df['diaObject_coord'][diao_index][0]\n",
    "sndec = df['diaObject_coord'][diao_index][1]\n",
    "sncoord = SkyCoord(snra, sndec, unit='deg')\n",
    "\n",
    "for i in range(3):\n",
    "    objra = galaxies['coord_ra'][i]\n",
    "    objdec = galaxies['coord_dec'][i]\n",
    "    objcoord = SkyCoord(objra, objdec, unit='deg')\n",
    "    del objra, objdec\n",
    "    \n",
    "    temp = objcoord.separation(sncoord)\n",
    "    galaxies['2Ddist'][i] = temp.arcsec\n",
    "    del temp\n",
    "    \n",
    "    temp = objcoord.spherical_offsets_to(sncoord)\n",
    "    xr = 3600.0 * temp[0].deg\n",
    "    yr = 3600.0 * temp[1].deg\n",
    "    del temp, objcoord\n",
    "    \n",
    "    Ixx = galaxies['shape_xx'][i]\n",
    "    Iyy = galaxies['shape_yy'][i]\n",
    "    Ixy = galaxies['shape_xy'][i]\n",
    "    Cxx = Iyy / ((Ixx * Iyy) - Ixy)\n",
    "    Cyy = Ixx / ((Ixx * Iyy) - Ixy)\n",
    "    Cxy = -2.0 * (Ixy) / ((Ixx * Iyy) - Ixy)\n",
    "    galaxies['ell_rad'][i] = np.sqrt((Cxx * xr**2) + (Cyy * yr**2) + (Cxy * xr * yr))\n",
    "\n",
    "    del Ixx, Iyy, Ixy, Cxx, Cyy, Cxy\n",
    "\n",
    "del snra, sndec, sncoord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2120bdf4-e342-4272-a413-2ed0400b5959",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "724016cf-cd1c-43e3-94b4-2e7a40ac0c8b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-30T19:49:49.824731Z",
     "iopub.status.busy": "2023-11-30T19:49:49.823859Z",
     "iopub.status.idle": "2023-11-30T19:49:50.316878Z",
     "shell.execute_reply": "2023-11-30T19:49:50.315464Z",
     "shell.execute_reply.started": "2023-11-30T19:49:49.824672Z"
    }
   },
   "source": [
    "### 3.4. Interpret derived data for nearby galaxies\n",
    "\n",
    "For the nearest galaxy by 2D sky separation, print the `objectId`, separations, and colors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "004cf4ab-1e37-4cdc-a3c9-448bd3c7ffe3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mx = np.argmin(galaxies['2Ddist'])\n",
    "print(galaxies['objectId'][mx], \n",
    "      galaxies['2Ddist'][mx], galaxies['ell_rad'][mx],\n",
    "      galaxies['gr_clr'][mx], galaxies['ri_clr'][mx])\n",
    "del mx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fc53e00-b8e5-4ec5-bd54-23ffa9942ea8",
   "metadata": {},
   "source": [
    "For the nearest galaxy by elliptical radii separation, print the `objectId`, separations, and colors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da925462-639a-48e8-a61c-31835e3af9ab",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mx = np.argmin(galaxies['ell_rad'])\n",
    "print(galaxies['objectId'][mx], \n",
    "      galaxies['2Ddist'][mx], galaxies['ell_rad'][mx],\n",
    "      galaxies['gr_clr'][mx], galaxies['ri_clr'][mx])\n",
    "del mx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e33b04ac-79c1-41d1-bb0f-c49931e50cd4",
   "metadata": {},
   "source": [
    "**Summary:** The galaxy that is nearest by 2D sky separation (3.15\") is not the best candidate host galaxy:\n",
    "the best candidate is the one with a larger 2D sky separation (4.62\") but a smaller offset in \n",
    "terms of elliptical radii ($R=1.02$), which takes into account the size of the galaxy.\n",
    "Furthermore, the best candidate also has redder colors, ($g-r$ and $r-i>0$), which is more typical \n",
    "for the host galaxies of Type Ia supernovae."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "856fe28d-09ef-42ca-ae58-5ff39dae8367",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "del diao_index, galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d131457-84f7-4002-b717-4c954333f5e3",
   "metadata": {},
   "source": [
    "## 4. Retrieve and display a light curve\n",
    "\n",
    "**In the future, the LSST alert packet will contain the light curves.**\n",
    "\n",
    "As a demo, retrieve the difference-image 5-sigma detections\n",
    "for one of the `DiaObjects` and plot the light curve.\n",
    "\n",
    "### 4.1. Set light curve plot parameters\n",
    "\n",
    "These color and symbol combinations are colorblind-friendly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d834969-a6d0-42f5-8360-92c9dad74625",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_filter_labels = ['u', 'g', 'r', 'i', 'z', 'y']\n",
    "plot_filter_colors = {'u': '#56b4e9', 'g': '#008060', 'r': '#ff4000',\n",
    "                      'i': '#850000', 'z': '#6600cc', 'y': '#000000'}\n",
    "plot_filter_symbols = {'u': 'o', 'g': '^', 'r': 'v', 'i': 's', 'z': '*', 'y': 'p'}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c82d9c24-6ba0-4020-be7a-dee1ae2b48dd",
   "metadata": {},
   "source": [
    "### 4.2. Retrieve the light curve for one `DiaObject`\n",
    "\n",
    "Descriptions and units for the columns of the `DiaSource` catalog are\n",
    "available in the <a href=\"https://dm.lsst.org/sdm_schemas/browser/dp02.html\">DP0.2 schema browser</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0595f622-e446-43b4-8379-023018bb32d0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "diao_index = 0\n",
    "temp_string = str(df['diaObjectId'][diao_index])\n",
    "diasources = rsp_tap.search(\"SELECT midPointTai, filterName, \"\n",
    "                            \"scisql_nanojanskyToAbMag(psFlux) AS psAbMag \"\n",
    "                            \"FROM dp02_dc2_catalogs.DiaSource  \"\n",
    "                            \"WHERE diaObjectId = \"+temp_string+\" \").to_table()\n",
    "del diao_index, temp_string"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30bc3ed1-6762-4b9b-84fa-b64dae81844e",
   "metadata": {},
   "source": [
    "Option to display the data retrieved from the Rubin Science Platform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e95d1543-dc67-4953-87c7-e7e88fdeff92",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# diasources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02f0cdaf-dccf-4a3f-b2d9-7c129197afd0",
   "metadata": {},
   "source": [
    "### 4.3. Display the light curve\n",
    "\n",
    "Plot the light curve with matplotlib. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77da5163-8a10-4fff-b079-bd4f88387ee6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "\n",
    "for f, filt in enumerate(plot_filter_labels):\n",
    "    fx = np.where(diasources['filterName'][:] == filt)[0]\n",
    "    if len(fx) > 0:\n",
    "        plt.plot(diasources['midPointTai'][fx]-60965, \n",
    "                 diasources['psAbMag'][fx],\n",
    "                 plot_filter_symbols[filt], \n",
    "                 ms=10, mew=0, alpha=0.5,\n",
    "                 color=plot_filter_colors[filt],\n",
    "                 label=plot_filter_labels[f])\n",
    "    del fx\n",
    "\n",
    "plt.xlabel('Days Ago (MJD-60965)')\n",
    "plt.ylabel('Apparent AB Magnitude')\n",
    "plt.gca().invert_yaxis()\n",
    "plt.title('5-Sigma Difference-Image Detections')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "038a7275-a2f9-477f-b9f5-228a33f774a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "del diasources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bc40ea6-f1b0-4b17-b1d6-2fad05238ef1",
   "metadata": {},
   "source": [
    "## 5. Retrieve and display a large, deep cutout image\n",
    "\n",
    "**In the future, the LSST alert packet will \n",
    "contain the difference-image and the reference-image stamps.**\n",
    "However, they are relatively small cutouts: no smaller than 30 x 30 pixels\n",
    "(6\" x 6\").\n",
    "\n",
    "Obtain larger version to explore the transient's environment.\n",
    "\n",
    "Select the first `DiaObject` on the list, and obtain its coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80ab2b2d-37ca-44bc-8160-52a8e843f59c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "diao_index = 0\n",
    "snra = df['diaObject_coord'][diao_index][0]\n",
    "sndec = df['diaObject_coord'][diao_index][1]\n",
    "print(snra, sndec)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0f23f56-f986-4334-ad73-a888e47af0d4",
   "metadata": {},
   "source": [
    "### 5.1. Query for the r-band deepCoadd image\n",
    "\n",
    "Create the query for all deeply coadded image (`deepCoadd_calexp`) that overlaps the `DiaObject` coordinates, and then select the r-band.\n",
    "\n",
    "It is recommended to always select all rows with `*` when querying the `ivoa.ObsCore`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4474913a-08b7-46f7-9a8a-b4fb54c32902",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = '''SELECT * FROM ivoa.ObsCore \n",
    "WHERE dataproduct_type = 'image' \n",
    "AND obs_collection = 'LSST.DP02' \n",
    "AND dataproduct_subtype = 'lsst.deepCoadd_calexp' \n",
    "AND CONTAINS(POINT('ICRS', {}, {}), s_region) = 1\n",
    "'''.format(snra, sndec)\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56c689d5-a72e-4baf-83cb-21896079d00f",
   "metadata": {},
   "source": [
    "Query the `ivoa.ObsCore`, which is the DP0.2 images available via the Rubin Science Platform deployed in the Google Cloud.\n",
    "\n",
    "Assert that the `results` table contains only six rows, one per filter.\n",
    "Adding `AND lsst_band = 'r'` to the query above would result in just one row, for the r-band filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e06a4c0-6b0d-45b7-b158-31ab664bf189",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = rsp_tap.search(query)\n",
    "assert len(results) == 6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac748b16-baf6-4391-b6ff-48dc071a50a9",
   "metadata": {},
   "source": [
    "Show the results in a table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d3dee06-2177-47df-8027-13f9f005ba8c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results.to_table().show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca488c69-0c77-4278-a0c5-72a1c9333b29",
   "metadata": {},
   "source": [
    "Store the results for r-band, which is the third row (index=2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c49e641-7f69-4f78-8150-809ab0a1ba8e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results_r = results[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e399fd67-54bd-4e51-9cce-bef063cf14da",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-12-01T01:27:27.214682Z",
     "iopub.status.busy": "2023-12-01T01:27:27.213941Z",
     "iopub.status.idle": "2023-12-01T01:27:27.222536Z",
     "shell.execute_reply": "2023-12-01T01:27:27.220897Z",
     "shell.execute_reply.started": "2023-12-01T01:27:27.214628Z"
    }
   },
   "source": [
    "The `access_url` column contains the URL for the retrievable image, \n",
    "and the `access_format` specifies the format (content type) of the data product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09f7d362-50cb-499a-a89d-430b862c59a3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(results_r.getdataurl())\n",
    "print(results_r.getdataformat())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "293102c2-fe1a-4a4f-a415-8db11f4b786b",
   "metadata": {},
   "source": [
    "The `access_format` indicates that the `access_url` is a <a href=\"https://www.ivoa.net/documents/DataLink/\">DataLink</a> service.\n",
    "DataLink is an IVOA data access protocol that provides a linking mechanism to metadata about a dataset, and the dataset itself."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0686660e-86c0-46d7-9839-7210d20d74e4",
   "metadata": {},
   "source": [
    "### 5.2. Use DataLink to retrieve the image information\n",
    "\n",
    "From `results_r`, get the data URL and print it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e9f27f1-591f-4a2d-87ff-a6558a4035d0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dl_url = results_r.getdataurl()\n",
    "print('Datalink link service URL: ', dl_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1949ab01-4f83-49af-bfe2-61cb2db34aef",
   "metadata": {},
   "source": [
    "Pass the Datalink ULR and the RSP credentials to `DatalinkResults`.\n",
    "\n",
    "Executing the following cell gets the infomation stored at the datalink URL for this image:\n",
    "a list matching records, each record containing a set of metadata describing the record.\n",
    "\n",
    "Display the results as a table in the notebook, and see that there\n",
    "are two rows: one contains the `access_url` for the exact storage\n",
    "location for the full image data file,\n",
    "and the other contains the `ID` for the cutout service (SODA) for this image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93f7fcff-521b-49e5-823e-6cd2f51f7d3e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dl_results = DatalinkResults.from_result_url(dl_url, session=credential)\n",
    "dl_results.to_table().show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "028439a6-7fc0-42c1-9dc1-52bdf7c2b470",
   "metadata": {},
   "source": [
    "### 5.3. Download and display the full image\n",
    "\n",
    "Download the full image, save it to file, and display it with matplotlib.\n",
    "\n",
    "First, get the Google-signed URL for the image, and store it in `image_url`. \n",
    "Note that this is temporary and will expire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f790eff-6410-425b-b957-13e6c8e12eda",
   "metadata": {},
   "outputs": [],
   "source": [
    "image_url = dl_results.getrecord(0).get('access_url')\n",
    "print(image_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a129c151-c6a8-4add-9142-f0d2cb4a2f85",
   "metadata": {},
   "source": [
    "Retrieve the image save it to the current working directory \n",
    "(`getcwd`) with the filename \"image.fits\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b10576bd-ea54-49f7-8b6a-ed705119f67a",
   "metadata": {},
   "outputs": [],
   "source": [
    "image_file = os.path.join(os.getcwd(), 'image.fits')\n",
    "urlretrieve(image_url, image_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa090d15-106f-42d7-9b88-c05e44e70397",
   "metadata": {},
   "source": [
    "Read the image header into `img_hdr`,\n",
    "store the WCS in `img_wcs`,\n",
    "and read the pixel data into `img_data`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36ce8169-4548-40f4-92ce-63aefccbebe4",
   "metadata": {},
   "outputs": [],
   "source": [
    "hdulist = fits.open(image_file)\n",
    "img_hdr = hdulist[1].header\n",
    "img_wcs = WCS(img_hdr)\n",
    "img_data = fits.getdata(image_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57a0c451-adc0-41f2-952e-7f78201f4467",
   "metadata": {},
   "source": [
    "Option to view the contents of the image header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aee6b12-890f-4ff2-9361-d0efd66352a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# img_hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3372be61-3d87-4d87-a864-a82c1876717d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-12-02T01:37:56.624848Z",
     "iopub.status.busy": "2023-12-02T01:37:56.624079Z",
     "iopub.status.idle": "2023-12-02T01:37:56.629245Z",
     "shell.execute_reply": "2023-12-02T01:37:56.628603Z",
     "shell.execute_reply.started": "2023-12-02T01:37:56.624823Z"
    }
   },
   "source": [
    "Get the pixel coordinates of the supernova."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a8e0a0d-90ba-459a-b3e7-8424155ad135",
   "metadata": {},
   "outputs": [],
   "source": [
    "sncoord = SkyCoord(snra, sndec, unit='deg')\n",
    "snxy = img_wcs.world_to_pixel(sncoord)\n",
    "print(snxy[0], snxy[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "378dbd3e-fca5-4b39-88c0-26c9189afad2",
   "metadata": {},
   "source": [
    "Display the image and mark the location of the supernova with a cyan cross.\n",
    "\n",
    "> **Warning:** A pink deprecation warning for matplotlib might appear below. It is ok to ignore in the context of this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aad729b1-d848-439e-8f78-8b38db1a9dc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(5, 5))\n",
    "plt.subplot(projection=img_wcs)\n",
    "plt.imshow(img_data, cmap='gray', \n",
    "           vmin=0, vmax=0.5, norm='linear',\n",
    "           origin='lower')\n",
    "plt.plot(snxy[0], snxy[1], 'o', ms=14, color='none', mec='cyan')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c22dec8c-1ca7-492b-80b0-2bbdb0cb6256",
   "metadata": {},
   "source": [
    "Clean up; only `dl_results` is needed for the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccd54679-6a75-432a-a0e3-163059289be3",
   "metadata": {},
   "outputs": [],
   "source": [
    "del query, results, results_r, dl_url, image_url\n",
    "del image_file, hdulist, img_hdr, img_wcs, img_data\n",
    "del sncoord, snxy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "594c858c-cc19-472d-85c9-a654316c2f59",
   "metadata": {},
   "source": [
    "### 5.3. Use SODA to make and retrieve a cutout image\n",
    "\n",
    "Use the PyVO `SodaQuery` function to make a small cutout \n",
    "centered on the supernova, and *then* \n",
    "retrieve, save, and display the small cutout,\n",
    "instead of retrieving the full image.\n",
    "\n",
    "Get a link to the cutout service (`cutoutService`) for the r-band image\n",
    "returned and stored in `dl_results`, above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3e2a457-d508-4331-af0e-3a748ecfe791",
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoutService = dl_results.get_adhocservice_by_id(\"cutout-sync\")\n",
    "cutoutService.params"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4074f633-256b-444c-b1b7-17dc02e1899c",
   "metadata": {},
   "source": [
    "Create a SODA query (`sodaQuery`) service from the Datalink results\n",
    "and the cutout service by passing the RSP credentials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d41ac21-6718-4d2e-9ba0-9bdbb91a77a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sodaQuery = SodaQuery.from_resource(dl_results, cutoutService, session=credential)\n",
    "assert sodaQuery is not None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0ca0af9-d763-4253-a943-df4db3728cd6",
   "metadata": {},
   "source": [
    "Define the center and radius for the cutout.\n",
    "Use the supernova's coordinates, and a radius\n",
    "of about 30 arcsec (0.008 deg)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f3edb3-a40f-4c6f-94cc-a61d32caa74d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sodaQuery.circle = (snra, sndec, 0.008)\n",
    "print(sodaQuery.circle)\n",
    "print(sodaQuery)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a453a7c-d0d6-4cda-9c6c-db833efaffbe",
   "metadata": {},
   "source": [
    "Define the location of the cutout image file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b859fbb-c0a4-4171-8d7a-f2884c87372a",
   "metadata": {},
   "outputs": [],
   "source": [
    "cutout_file = os.path.join(os.getcwd(), 'image-cutout.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d4dd80f-7f5d-4ab9-ac9b-7656c00abe9f",
   "metadata": {},
   "source": [
    "Execute the defined SODA query to make the cutout, \n",
    "read the results (the cutout image), and\n",
    "save the cutout to the `cutout_file`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb8b0652-a526-483a-8fe1-fb9df733174d",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(cutout_file, 'bw') as f:\n",
    "    f.write(sodaQuery.execute_stream().read())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83d743c2-1f43-4819-b1c8-e0a70299ff65",
   "metadata": {},
   "source": [
    "Obtain the cutout image's header, WCS, and data.\n",
    "\n",
    "Convert the supernova's coordinates to pixel coordinates for the cutout image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e0a3098-dfec-4434-8554-178095b3fcfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "hdulist = fits.open(cutout_file)\n",
    "coimg_hdr = hdulist[1].header\n",
    "coimg_wcs = WCS(coimg_hdr)\n",
    "coimg_data = fits.getdata(cutout_file)\n",
    "sncoord = SkyCoord(snra, sndec, unit='deg')\n",
    "snxy = coimg_wcs.world_to_pixel(sncoord)\n",
    "print(snxy[0], snxy[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a95787-1d2e-424f-93f6-f028a2eb81ee",
   "metadata": {},
   "source": [
    "Display the cutout image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "501327e4-95b4-488f-b533-0916878ded2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(5, 5))\n",
    "plt.subplot(projection=coimg_wcs)\n",
    "plt.imshow(coimg_data, cmap='gray', \n",
    "           vmin=-0.1, vmax=0.5, norm='linear',\n",
    "           origin='lower')\n",
    "plt.plot(snxy[0], snxy[1], 'o', ms=14, color='none', mec='cyan')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a8cd0aa-0d7a-44de-875d-4cabe8743986",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}