DP02_01_Introduction_to_DP02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {
    "tags": []
   },
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<br><b>Introduction to Jupyter Notebooks for Data Preview 0.2</b> <br>\n",
    "Contact author: Melissa Graham <br>\n",
    "Last verified to run: 2024-05-01 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_16 <br>\n",
    "Container size: medium <br>\n",
    "Targeted learning level: beginner <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf7b00a0-a005-4a22-a4c7-a6ecb32c749a",
   "metadata": {},
   "source": [
    "**Description:** An introduction to using Jupyter Notebooks and Rubin Python packages to access LSST data products (images and catalogs)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ed8dff8-bcf7-4508-a2f2-cd2260321261",
   "metadata": {},
   "source": [
    "**Skills:** Execute Python code in a Jupyter Notebook. Use the TAP service to retrieve Object catalog data. Use the Butler to retrieve and display a deepCoadd image."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a748b3a3-00e8-4215-adfa-adcabff26f56",
   "metadata": {},
   "source": [
    "**LSST Data Products:** TAP dp02_dc2_catalogs.Object table. Butler deepCoadd image."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "387f46c8-d58f-400f-a6e6-f21c9ae7ce94",
   "metadata": {},
   "source": [
    "**Packages:** lsst.rsp.get_tap_service, lsst.rsp.retrieve_query, lsst.daf.butler, lsst.afw.display, lsst.geom, pandas, matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d622b84-83fa-4de9-94b6-d16af645a37e",
   "metadata": {},
   "source": [
    "**Credit:** Originally developed by Melissa Graham and the Rubin Community Science Team in the context of the Rubin DP0.1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "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": "cfc73be0",
   "metadata": {},
   "source": [
    "## 1.0. Introduction\n",
    "\n",
    "This Jupyter Notebook provides an introduction to how notebooks work. It demonstrates how to execute code and markdown text cells, how to import Python packages and learn about their modules, and provides links to further documentation.\n",
    "\n",
    "This Notebook also demonstrates the basic functionality of the Rubin Science Platform (RSP) installed at the Interim Data Facility (IDF; the Google Cloud), such as how to use the TAP service to query and retrieve catalog data; matplotlib to plot catalog data; the LSST Butler package to query and retrieve image data; and the LSST afwDisplay package to display images.\n",
    "\n",
    "This Notebook uses the Data Preview 0.2 (DP0.2) data set. This data set uses a subset of the DESC's Data Challenge 2 (DC2) simulated images, which have been *reprocessed* by Rubin Observatory using Version 23 of the LSST Science Pipelines. More information about the simulated data can be found in the <a href=\"https://ui.adsabs.harvard.edu/abs/2021ApJS..253...31L/abstract\">DESC's DC2 paper</a> and in the <a href=\"https://dp0-2.lsst.io\">DP0.2 data release documentation</a>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3507e32",
   "metadata": {},
   "source": [
    "### 1.1. How to use a Jupyter Notebook\n",
    "\n",
    "Jupyter Notebooks contain a mix of code, output, visualizations, and narrative text. The most comprehensive source for documentation about Jupyter Notebooks is https://jupyter-notebook.readthedocs.io, but there are many great beginner-level tutorials and demos out there. Usually a web search of a question, like \"how to make a table in markdown jupyter notebook\", will yield several good examples. Often the answers will be found in <a href=\"https://stackoverflow.com/\">StackOverflow</a>.\n",
    "\n",
    "A Jupyter Notebook is a series of cells. There are three types of cells: code, markdown, and raw. This text was generated from a markdown cell. Up in the menu bar you will find a drop-down menu to set the cell type."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3db41666-4e0a-4a7e-93a5-90b99286fedf",
   "metadata": {},
   "source": [
    "> **Warning:** All of the code cells in a notebook should be executed in the order that they appear."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fb67584-1bc1-4049-9f93-fa0f45d477c0",
   "metadata": {},
   "source": [
    "Click in the following code cell: with the cursor in the cell, simultaneously press \"shift\" and \"enter\" (or \"return\") to execute the cell code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba94de19",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# This is a code cell. Press shift-enter to execute.\n",
    "# The # makes these lines comments, not code. They are not executed.\n",
    "print('Hello, world!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06a6f7ed-284b-4a35-b0ba-d95a1300d463",
   "metadata": {},
   "source": [
    "<!---\n",
    "This is a markdown cell.\n",
    "Press shift-enter to execute, and see the formatted text reappear.\n",
    "-->\n",
    "\n",
    "Double click on THESE WORDS IN THIS MARKDOWN CELL to see the markdown source code."
   ]
  },
  {
   "cell_type": "raw",
   "id": "a3449766",
   "metadata": {},
   "source": [
    "# This is a raw cell. You can press shift-enter, but nothing will execute.\n",
    "# The # symbol does not mean anything in a raw cell.\n",
    "print('Hello, world!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a34975f6",
   "metadata": {},
   "source": [
    "#### 1.1.1. Jupyter Notebooks How-Tos\n",
    "\n",
    "**How to quickly execute all the cells:** \n",
    "Go to the top menu bar and select \"Kernel\", then \"Restart Kernel and Run All Cells\".\n",
    "\n",
    "**How to emergency-stop a notebook:** \n",
    "If a code cell is taking a long time to execute (e.g., if a process to retrieve an entire catalog was started by accident)\n",
    "kill it by going to \"Kernel\" in the top menu bar and selecting \"Restart Kernel and Clear All Outputs\".\n",
    "It might still a few tens of seconds, but it will stop the process and restart the kernel.\n",
    "\n",
    "**The kernel** is the computational engine for the notebook (the RSP uses a `python3` kernel),\n",
    "and can be thought of as a live compiler. \n",
    "Restarting the kernel and clearning all outputs means that all defined variables or functions are removed from memory, \n",
    "and all code cells revert to an \"unexecuted\" state.\n",
    "\n",
    "**How to view a table of contents for this notebook:** \n",
    "Click on the icon of a bullet list in the leftmost vertical menu bar, and an automatically-generated ToC will appear at left. \n",
    "Click on the icon of the file folder at the top of the leftmost vertical menu bar to return to a directory view.\n",
    "\n",
    "**How to know which version of the LSST Science Pipelines is running:** \n",
    "Look along the bottom bar of this browser window, and find the version of the LSST Science Pipelines that was selected as the \"image\". \n",
    "It is probably \"Recommended (Weekly yyyy_ww)\", and it should match the verified version listed in the notebook's header. \n",
    "Alternatively, uncomment the two lines in the following code cell and execute the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48e9d2d7-ebec-42d9-8de0-73c1475fd539",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# ! echo $IMAGE_DESCRIPTION\n",
    "# ! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "245d2838",
   "metadata": {},
   "source": [
    "### 1.2. Package Imports\n",
    "\n",
    "Most Jupyter Notebooks start out by importing all the packages they will need in the first code cell.\n",
    "\n",
    "Complete knowledge of these packages is not required in order to complete this tutorial, but here is a bit of basic information and some links for further learning.\n",
    "\n",
    "**numpy**: A fundamental package for scientific computing with arrays in Python. Its comprehensive documentation is available at <a href=\"https://numpy.org\">numpy.org</a>, and it includes quickstart beginner guides. (The numpy package is not used in this notebook, but is imported as a demonstration because it is a very commonly-used package.) <br>\n",
    "\n",
    "**matplotlib**: This package is a comprehensive library for creating static, animated, and interactive visualizations in Python. Its comprehensive documentation is at <a href=\"https://matplotlib.org/\">matplotlib.org</a>. The <a href=\"https://matplotlib.org/stable/gallery/index.html\">matplotlib gallery</a> is a great place to start and links to examples. <br>\n",
    " \n",
    "**pandas**: A package which allows users to deal efficiently with tabular data in dataframes. Learn more in the <a href=\"https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.html\">Pandas documentation</a>. <br>\n",
    "\n",
    "**astropy**: A Python package of useful astronomy tools. Learn more in the <a href=\"http://docs.astropy.org/en/stable/_modules/astropy/table/table.html\">astropy documentation</a>.\n",
    " \n",
    "**lsst**: These packages are all from the <a href=\"https://pipelines.lsst.io/\">LSST Science Pipelines</a>.\n",
    "The `lsst.rsp` package enables image and catalog access via the TAP service (see Section 2); \n",
    "the `lsst.daf.butler` package enables image and catalog access via the butler (see Section 3);\n",
    "and the `lsst.geom` has helper functions for image metadata and `lsst.afw.display` package enables image display (see Section 3)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8106313",
   "metadata": {},
   "source": [
    "Import the packages used in this notebook by executing the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddc1458",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas\n",
    "from lsst.rsp import get_tap_service, retrieve_query\n",
    "import lsst.daf.butler as dafButler\n",
    "import lsst.geom\n",
    "import lsst.afw.display as afwDisplay"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5ca54f6",
   "metadata": {},
   "source": [
    "#### 1.2.1. Learn more about the imported Python packages\n",
    "\n",
    "Print the version of numpy and matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ef89712-485d-4016-93df-4b3627bc3415",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print('numpy version: ', numpy.__version__)\n",
    "print('matplotlib version: ', matplotlib.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2649d146-4903-4c47-95f3-640970b55618",
   "metadata": {},
   "source": [
    "View a pop-up list of any package's modules by writing the package name, then a period, and then pressing tab. Use the up and down arrows to scroll through the pop-up list. This works whether or not the line is commented-out. In the cell below, `numpy.` is commented-out because that is not an executable code statement, and if the # were not there, this cell would fail to execute (try it -- remove the #, press shift-enter, and watch it fail)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ff354c2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# numpy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a4a559b",
   "metadata": {},
   "source": [
    "Use \"help\" function to view the help documentation for a package. Remove the # symbol to un-comment any one line, and execute the following cell. Help documentation can be really long. Re-comment the line by replacing the #, then re-execute the cell and the output will go away."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2d31bd7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# help(numpy)\n",
    "# help(matplotlib)\n",
    "# help(numpy.abs)\n",
    "# help(matplotlib.pyplot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec51ac0b",
   "metadata": {},
   "source": [
    "## 2.0. Catalog Data\n",
    "\n",
    "### 2.1. Table Access Protocol (TAP) service\n",
    "\n",
    "Table Access Procotol (TAP) provides standardized access to the catalog data for discovery, search, and retrieval. Full <a href=\"http://www.ivoa.net/documents/TAP\">documentation for TAP</a> is provided by the International Virtual Observatory Alliance (IVOA).\n",
    "\n",
    "The TAP service uses a query language similar to SQL (Structured Query Langage) called ADQL (Astronomical Data Query Language). The <a href=\"http://www.ivoa.net/documents/latest/ADQL.html\">documentation for ADQL</a> includes more information about syntax and keywords.\n",
    "\n",
    "> **Notice:** Not all ADQL functionality is supported by the RSP for Data Preview 0. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4936e520",
   "metadata": {},
   "source": [
    "Start the TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0a530b9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10804bb0",
   "metadata": {},
   "source": [
    "### 2.2. Exploring catalog tables and columns with TAP\n",
    "\n",
    "This example uses the DP0.2 Object catalog, which contains sources detected in the coadded images (also called stacked, combined, or deepCoadd images). \n",
    "\n",
    "Catalog contents can also be explored with the <a href=\"https://dm.lsst.org/sdm_schemas/browser/dp02.html\">DP0.2 schema browser</a>.\n",
    "\n",
    "Results from a TAP service search are best displayed using one of two functions:<br>\n",
    "`.to_table()`: convert results to an astropy table. <br>\n",
    "`.to_table().to_pandas()`: convert to an astropy table and then to a Pandas dataframe.\n",
    "\n",
    "> **Warning:** do not use the .to_table().show_in_notebook() method. This can cause issues in the RSP Jupyterlab environment that make your notebook hang indefinitely.\n",
    "\n",
    "The three optional exercises below teach different ways to explore using the TAP service. They show how to use the TAP service with ADQL statements to discover what catalogs exist and which columns they contain. Each cell uses a different method to display the TAP search results. Remove all of the # and execute each cell, and see that they create a lot of output -- add the # back to each line and re-execute the cell, and the output will go away."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beef6193",
   "metadata": {},
   "source": [
    "#### 2.2.1. Exercise 1\n",
    "Retrieve and display a list of all the table names and descriptions that are available via the TAP server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75ba75e6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# my_adql_query = \"SELECT description, table_name FROM TAP_SCHEMA.tables\"\n",
    "# results = service.search(my_adql_query)\n",
    "# results_table = results.to_table().to_pandas()\n",
    "# results_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "723cbe86",
   "metadata": {},
   "source": [
    "#### 2.2.2. Exercise 2\n",
    "Retrieve and display a list of the field names (columns names) in the DP0.2 Object catalog's TAP schema. Note that the results can be named anything else; here, 'res' is used instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88711ca7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# my_adql_query = \"SELECT * from TAP_SCHEMA.columns \"+\\\n",
    "#                 \"WHERE table_name = 'dp02_dc2_catalogs.Object'\"\n",
    "# res = service.search(my_adql_query)\n",
    "# print(res.fieldnames)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de30cca5",
   "metadata": {},
   "source": [
    "#### 2.2.3. Exercise 3\n",
    "Retrieve the names, data types, description, and units for all columns in the Object catalog. Display the number of columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5129c644",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# my_adql_query = \"SELECT column_name, datatype, description, unit \"+\\\n",
    "#                 \"FROM TAP_SCHEMA.columns \"+\\\n",
    "#                 \"WHERE table_name = 'dp02_dc2_catalogs.Object'\"\n",
    "# results = service.search(my_adql_query)\n",
    "# results_table = results.to_table().to_pandas()\n",
    "# print('Number of columns available in the Object catalog: ', len(results_table))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8f5cb6e-1844-4f0f-9329-cb20f20617be",
   "metadata": {},
   "source": [
    "Display all 991 column names and their information. It's so much output! Comment-out the line in the cell and re-execute the cell to make all that output disappear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2142c9a-3be3-4e61-b69e-47fd25d23178",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# results_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2577a97a-4a96-41bb-8e29-4cae62672b3d",
   "metadata": {},
   "source": [
    "Only display names and descriptions for columns that contain the string \"cModelFlux\". \n",
    "Try other strings like \"coord\", \"extendedness\", \"deblend\", or \"detect\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40bde1ad-15cd-409c-9bcb-162052ee7890",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# my_string = 'cModelFlux'\n",
    "# for col,des in zip(results_table['column_name'],results_table['description']):\n",
    "#     if col.find(my_string) > -1:\n",
    "#         print('%-40s %-200s' % (col,des))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ced15d7",
   "metadata": {},
   "source": [
    "### 2.3. Retrieving data with TAP\n",
    "\n",
    "A few tips about how to do efficient queries on the DP0.2 catalogs.\n",
    "\n",
    "**RA, Dec constraints yield faster queries:**\n",
    "LSST Query Services (Qserv) provides access to the LSST Database Catalogs.\n",
    "Users can query the catalogs using standard SQL query language with a few restrictions.\n",
    "Qserv stores catalog data sharded by coordinate (RA, Dec).\n",
    "ADQL query statements that include constraints by coordinate do not require a whole-catalog search, and are typically faster (and can be much faster) than ADQL query statements which only include constraints for other columns.\n",
    "\n",
    "**Retrieve a small sample of rows:**\n",
    "As demonstrated in Section 2.3.2, use `maxrec=10` or `SELECT TOP 10` when exploring data sets in order to return only a few rows to play with (this can also shorten query times for exploratory queries without WHERE statements).\n",
    "\n",
    "**Recommended constraint on `detect_isPrimary`:**\n",
    "When applicable, it is recommended to include `detect_isPrimary = True` in queries for the `Object`, `Source`, and `ForcedSource` catalogs.\n",
    "This parameter is `True` if a source has no children, is in the inner region of a coadd patch, is in the inner region of a coadd tract, and is not detected in a pseudo-filter.\n",
    "Including this constraint will remove any duplicates (i.e., will not include both a parent and its deblended children)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0370bb6e-6c3d-4cd2-978a-4b00fd485478",
   "metadata": {},
   "source": [
    "#### 2.3.1. Converting fluxes to magnitudes\n",
    "\n",
    "The object and source catalogs store only fluxes. There are hundreds of flux-related columns, and to store them also as magnitudes would be redundant, and a waste of space.\n",
    "\n",
    "All flux units are nanojanskys ($nJy$). The <a href=\"https://en.wikipedia.org/wiki/AB_magnitude\">AB Magnitudes Wikipedia page</a> provides a concise resource for users unfamiliar with AB magnitudes and jansky fluxes. To convert $nJy$ to AB magnitudes use: $m_{AB} = -2.5log( f_{nJy}) + 31.4$.\n",
    "\n",
    "As demonstrated in Section 2.3.2, to add columns of magnitudes after retrieving columns of flux, users can do this:<br>\n",
    "`results_table['r_calibMag'] = -2.50 * numpy.log10(results_table['r_calibFlux']) + 31.4`<br>\n",
    "`results_table['r_cModelMag'] = -2.50 * numpy.log10(results_table['r_cModelFlux']) + 31.4`\n",
    "\n",
    "As demonstrated in Section 2.3.3, to retrieve columns of fluxes *as magnitudes* in an ADQL query, users can do this:<br>\n",
    "`scisql_nanojanskyToAbMag(g_calibFlux) as g_calibMag`,\n",
    "and columns of magnitude errors can be retrieved with:<br>\n",
    "`scisql_nanojanskyToAbMagSigma(g_calibFlux, g_calibFluxErr) as g_calibMagErr`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5356f2e2-a369-42d2-8d46-ff155e221d7e",
   "metadata": {},
   "source": [
    "#### 2.3.2. Ten objects of any kind\n",
    "\n",
    "To quickly demonstrate how to retrieve data from the Object catalog,  use a cone search and request only 10 records be returned. Figure 2 of the <a href=\"https://ui.adsabs.harvard.edu/abs/2021ApJS..253...31L/abstract\">DESC's DC2 paper</a> shows the sky region covered by the DC2 simulation contains coordinates RA,Dec = 62,-37.\n",
    "\n",
    "This example uses `maxrec=10` in the `service.search()` function, but the same results could be achieved by replacing `SELECT` with `SELECT TOP 10` in the ADQL query.\n",
    "\n",
    "> **Warning:** The Object catalog contains hundreds of millions of rows. Searches that do not specify a region and/or a maximum number of records can take a long time, and return far too many rows to display in a notebook.\n",
    "\n",
    "Retrieve coordinates and g,r,i magnitudes for 10 objects within a radius 0.5 degrees of 62,-37."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fa7d255-584b-4be5-94a3-2b901147c4c3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "use_center_coords = \"62, -37\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c23a25b7-e861-4d74-a07d-cde752b10483",
   "metadata": {},
   "source": [
    "Below, `SELECT TOP 10` is used in the query statement to limit the returned data to 10 objects.\n",
    "An alternative is to use the `maxrec` keyword in the search statement: `service.search(my_adql_query, maxrec=10)`.\n",
    "However, use of `maxrec` might return a `DALOverflowWarning` to warn the user that partial results have been returned (even though partial results were desired)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "355a5cb2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_adql_query = \"SELECT TOP 10 \"+ \\\n",
    "                \"coord_ra, coord_dec, detect_isPrimary, \" + \\\n",
    "                \"r_calibFlux, r_cModelFlux, r_extendedness \" + \\\n",
    "                \"FROM dp02_dc2_catalogs.Object \" + \\\n",
    "                \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \" + \\\n",
    "                \"CIRCLE('ICRS', \" + use_center_coords + \", 0.01)) = 1 \"\n",
    "\n",
    "results = service.search(my_adql_query)\n",
    "results_table = results.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb72874d-8c81-44a4-a1f8-c37b117f734e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results_table['r_calibMag'] = -2.50 * numpy.log10(results_table['r_calibFlux']) + 31.4\n",
    "results_table['r_cModelMag'] = -2.50 * numpy.log10(results_table['r_cModelFlux']) + 31.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a863831-690c-4cb9-b48b-2e49256863b6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1f5b222",
   "metadata": {},
   "source": [
    "#### 2.3.2. Ten thousand point-like objects\n",
    "\n",
    "In addition to a cone search, impose query restrictions that detect_isPrimary is True (this will not return deblended \"child\" sources), that the calibrated flux is greater than 360 nJy (about 25th mag), and that the extendedness parameters are 0 (point-like sources).\n",
    "\n",
    "Retrieve g-, r- and i-band magnitudes for 10000 objects that are likely to be stars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9989da7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT TOP 10000 \"\n",
    "                         \"coord_ra, coord_dec, \"\n",
    "                         \"scisql_nanojanskyToAbMag(g_calibFlux) as g_calibMag, \"\n",
    "                         \"scisql_nanojanskyToAbMag(r_calibFlux) as r_calibMag, \"\n",
    "                         \"scisql_nanojanskyToAbMag(i_calibFlux) as i_calibMag, \"\n",
    "                         \"scisql_nanojanskyToAbMagSigma(g_calibFlux, g_calibFluxErr) as g_calibMagErr \"\n",
    "                         \"FROM dp02_dc2_catalogs.Object \"\n",
    "                         \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\n",
    "                         \"CIRCLE('ICRS', \"+use_center_coords+\", 1.0)) = 1 \"\n",
    "                         \"AND detect_isPrimary = 1 \"\n",
    "                         \"AND g_calibFlux > 360 \"\n",
    "                         \"AND r_calibFlux > 360 \"\n",
    "                         \"AND i_calibFlux > 360 \"\n",
    "                         \"AND g_extendedness = 0 \"\n",
    "                         \"AND r_extendedness = 0 \"\n",
    "                         \"AND i_extendedness = 0\")\n",
    "\n",
    "results_table = results.to_table()\n",
    "print(len(results_table))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d293260a-e7d3-4fdb-9a37-716f612c6220",
   "metadata": {},
   "source": [
    "The table display will automatically truncate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13eeaa07-ca19-4b1e-983f-3e29908adbc2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41b68c33",
   "metadata": {},
   "source": [
    "Put the results into a pandas dataframe for easy access to contents. This data is used to create a color-magnitude diagram in Section 2.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "798b0809",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data = results_table.to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e1dc71c",
   "metadata": {},
   "source": [
    "For users unfamiliar with Pandas, here are some optional lines of code that demonstrate how to print the column names, the 'ra' column info, or the 'ra' column values. Uncomment (remove #) and execute the cell to view the demo output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd0f247a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# data.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3dd33a7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# data['coord_ra']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a72e6ecd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# data['coord_ra'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd5b98bc",
   "metadata": {},
   "source": [
    "### 2.5 Make a color-magnitude diagram\n",
    "\n",
    "Use the plot task of the matplotlib.pyplot package (which was imported as plt). The plot task parameters are described in full at <a href=\"https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html#matplotlib.pyplot.plot\">this matplotlib website</a>, but briefly they are: x values, y values, symbol shape ('o' is circle), marker size (ms), and marker transparency (alpha)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98b2c106",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.plot(data['r_calibMag'].values - data['i_calibMag'].values,\n",
    "         data['g_calibMag'].values, 'o', ms=2, alpha=0.2)\n",
    "\n",
    "plt.xlabel('mag_r - mag_i', fontsize=16)\n",
    "plt.ylabel('mag_g', fontsize=16)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "\n",
    "plt.xlim([-0.5, 2.0])\n",
    "plt.ylim([25.5, 16.5])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c2c40f7",
   "metadata": {},
   "source": [
    "This plot generates many questions, such as \"Why are the colors quantized?\" and \"Are those all really stars?\". The answers are beyond the scope of this notebook, and are left as potential topics of scientific analysis that could be done with the DC2 data set."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7c27a75-8928-4ead-a958-c7142b0a9a6e",
   "metadata": {},
   "source": [
    "### 2.6 Optional: plot magnitude versus magnitude error\n",
    "\n",
    "To illustrate both the magnitudes and magnitude errors retrieved via the TAP query above,\n",
    "here is an option to plot the magnitude error as a function of magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9503291-d974-471a-9d53-fd880f57e6fd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# plt.plot(data['g_calibMag'].values, data['g_calibMagErr'].values, 'o', ms=2, alpha=0.2)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56baaa33",
   "metadata": {},
   "source": [
    "## 3.0. Image Data\n",
    "\n",
    "The two most common types of images that DP0 delegates will interact with are calexps and deepCoadds.\n",
    "\n",
    "**calexp**: A single image in a single filter.\n",
    "\n",
    "**deepCoadd**: A combination of single images into a deep stack or Coadd.\n",
    " \n",
    "The LSST Science Pipelines processes and stores images in tracts and patches.\n",
    "\n",
    "**tract**: A portion of sky within the LSST all-sky tessellation (sky map); divided into patches.\n",
    "\n",
    "**patch**: A quadrilateral sub-region of a tract, of a size that fits easily into memory on desktop computers.\n",
    " \n",
    "To retrieve and display an image at a desired coordinate, users have to specify their image type, tract, and patch."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6758b295",
   "metadata": {},
   "source": [
    "### 3.1. Create an instance of the butler\n",
    "\n",
    "The butler (<a href=\"https://pipelines.lsst.io/modules/lsst.daf.butler/index.html\">documentation</a>) is an LSST Science Pipelines software package to fetch LSST data without having to know its location or format. The butler can also be used to explore and discover what data exist. Other tutorials demonstrate the full butler functionality."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b640a69",
   "metadata": {},
   "source": [
    "Create an instance of the butler using the following DP0.2 configuration and collection.\n",
    "It will return an informative statement about credentials being found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2122a1fb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "butler = dafButler.Butler('dp02', collections='2.2i/runs/DP0.2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f25b7df-a901-4a69-81d0-c630cec28539",
   "metadata": {},
   "source": [
    "### 3.2. Identify and retrieve a deepCoadd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0519417",
   "metadata": {},
   "source": [
    "There is a cool-looking DC2 galaxy cluster at RA = 03h42m59.0s, Dec = -32d16m09s (in degrees, 55.745834, -32.269167).\n",
    "\n",
    "Use lsst.geom to define a SpherePoint for the cluster's coordinates (<a href=\"https://pipelines.lsst.io/modules/lsst.geom/index.html\">lsst.geom documentation</a>)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56ea8e01",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_ra_deg = 55.745834\n",
    "my_dec_deg = -32.269167\n",
    "\n",
    "my_spherePoint = lsst.geom.SpherePoint(my_ra_deg*lsst.geom.degrees,\n",
    "                                       my_dec_deg*lsst.geom.degrees)\n",
    "print(my_spherePoint)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "789fb4ff",
   "metadata": {},
   "source": [
    "Retrieve the DC2 sky map from the butler and use it to identify the tract and patch for the cluster's coordinates (<a href=\"http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/skymap.html\">skymap documentation</a>)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1796f237",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "skymap = butler.get('skyMap')\n",
    "\n",
    "tract = skymap.findTract(my_spherePoint)\n",
    "patch = tract.findPatch(my_spherePoint)\n",
    "\n",
    "my_tract = tract.tract_id\n",
    "my_patch = patch.getSequentialIndex()\n",
    "\n",
    "print('my_tract: ', my_tract)\n",
    "print('my_patch: ', my_patch)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63c0e237",
   "metadata": {},
   "source": [
    "Use the butler to retrieve the deep i-band Coadd for the tract and patch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "460e1df6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dataId = {'band': 'i', 'tract': my_tract, 'patch': my_patch}\n",
    "my_deepCoadd = butler.get('deepCoadd', dataId=dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4e2c4d1",
   "metadata": {},
   "source": [
    "### 3.3. Display the image with afwDisplay\n",
    "Image data retrieved with the butler can be displayed several different ways. A simple option is to use the LSST Science Pipelines package afwDisplay. There is some <a href=\"https://pipelines.lsst.io/modules/lsst.afw.display/index.html\">documentation for afwDisplay</a> available, and other DP0 tutorials go into more detail about all the display options (e.g., overlaying mask data to show bad pixels).\n",
    "\n",
    "Set the backend of afwDisplay to matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4ebe015",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f1ab1a8",
   "metadata": {},
   "source": [
    "Use afwDisplay to show the image data retrieved.\n",
    "\n",
    "The following code cell creates a matplotlib.pyplot figure; aliases `lsst.afw.display.Display` as `afw_display`; \n",
    "sets the scale for the pixel shading; displays the image data using `mtv`; and turns on the x and y axes labels (pixel coordinates)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd46970e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 8))\n",
    "afw_display = afwDisplay.Display(1)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(my_deepCoadd.image)\n",
    "plt.gca().axis('on')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a318232",
   "metadata": {},
   "source": [
    "To learn more about the afwDisplay package and its tasks, use the help function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e9a1ee0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# help(afw_display.scale)\n",
    "# help(afw_display.mtv)"
   ]
  }
 ],
 "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.7"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}