04_Intro_to_Butler.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<b>Introduction to the LSST data Butler</b> <br>\n",
    "Last verified to run on 2021-10-12 with LSST Science Pipelines release w_2021_40 <br>\n",
    "Contact authors: Alex Drlica-Wagner, Melissa Graham <br>\n",
    "Target audience: All DP0 delegates. <br>\n",
    "Container Size: medium <br>\n",
    "Questions welcome at <a href=\"https://community.lsst.org/c/support/dp0\">community.lsst.org/c/support/dp0</a>. <br>\n",
    "Find DP0 documentation and resources at <a href=\"https://dp0-1.lsst.io\">dp0-1.lsst.io</a>. <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Credit:** This tutorial was originally developed by Alex Drlica-Wagner in the context of the LSST Stack Club. Please consider acknowledging Alex Drlica-Wagner in any publications or software releases that make use of this notebook's contents."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "\n",
    "The Butler is the LSST Science Pipelines interface for managing, reading, and writing datasets. The Butler can be used to explore the contents of the DP0.1 data repository and access the DP0.1 data. The current version of the Butler (referred to as \"Gen-3\") is still under development, and this notebook may be modified in the future. Full Butler documentation can be found [here](https://pipelines.lsst.io/v/weekly/modules/lsst.daf.butler/index.html).\n",
    "\n",
    "This notebook demonstrates how to:<br>\n",
    "1. Create an instance of the Butler<br>\n",
    "2. Explore the DP0.1 data repository<br>\n",
    "3. Retrieve and plot an image with sources<br>\n",
    "4. Create an image cutout at a user-specified coordinate<br>\n",
    "5. Exploring and retrieving catalog data from the Butler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Up "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import general python packages and several packages from the LSST Science Pipelines, including the Butler package and AFW Display, which will be used to display images.\n",
    "More details and techniques regarding image display with `AFW Display` can be found in the `rubin-dp0` GitHub Organization's [tutorial-notebooks](https://github.com/rubin-dp0/tutorial-notebooks) repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generic python packages\n",
    "import warnings\n",
    "import numpy as np\n",
    "import pylab as plt\n",
    "\n",
    "# Set a standard figure size to use\n",
    "plt.rcParams['figure.figsize'] = (8.0, 8.0)\n",
    "\n",
    "# LSST Science Pipelines (Stack) packages\n",
    "import lsst.daf.butler as dafButler\n",
    "import lsst.afw.display as afwDisplay\n",
    "import lsst.geom as geom\n",
    "import lsst.afw.coord as afwCoord\n",
    "afwDisplay.setDefaultBackend('matplotlib')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prevent some helpful but ancillary warning messages from printing\n",
    "#   during some LSST DM Release calls\n",
    "warnings.simplefilter(\"ignore\", category=UserWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This should match the verified version listed at the start of the notebook\n",
    "! echo $IMAGE_DESCRIPTION\n",
    "! eups list -s lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Create an instance of the Butler\n",
    "\n",
    "To create the Butler, we need to provide it with a path to the data set, which is called a \"data repository\".\n",
    "Butler repositories have both a database component and a file-like storage component; the latter can can be remote (i.e., pointing to an S3 bucket) or local (i.e., pointing to a directory on the local file system), and it contains a configuration file (usually `butler.yaml`) that points to the right database\n",
    "\n",
    "S3 (Simple Storage Service) buckets are public cloud storage resources that are similar to file folders, store objects, and which consist of data and its descriptive metadata.\n",
    "\n",
    "For access to the DP0.1 data set, always point to this S3 bucket."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repo = 's3://butler-us-central1-dp01'\n",
    "butler = dafButler.Butler(repo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Explore the DP0.1 data repository"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1 Butler registry and collections\n",
    "\n",
    "The database side of a data repository is called a registry.\n",
    "The registry contains entries for all data products, and organizes them by _collection_, _dataset type_, and _data ID_.\n",
    "Use the registry to investigate a repository by listing all collections.\n",
    "\n",
    "Find more about collections [here](https://pipelines.lsst.io/v/weekly/modules/lsst.daf.butler/organizing.html#collections).\n",
    "\n",
    "A registry client is part of our butler object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "registry = butler.registry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Learn more about the registry by uncommenting the following line.\n",
    "# help(registry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for c in sorted(registry.queryCollections()):\n",
    "    print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "Here are some definitions to help delegates understand the contents of the full DP0.1 repository. \n",
    "\n",
    "* `2.2i` - refers to the processing run of the LSST DESC DC2 data (the `i` stands for `imSim`)\n",
    "* `calib` - refers to calibration products that are used for instrument signature removal\n",
    "* `runs` - refers to processed data products\n",
    "* `refcats` - refers to the reference catalogs used for astrometric and photometric calibration\n",
    "* `skymaps` - definitions for the _tract_ and _patch_ grids that coadds are built on\n",
    "\n",
    "Some collections are nested, and DP0 delegates can access all the data for DC2 Run 2.2i, which is the DP0.1 data set, by selecting the collection `2.2i/runs/DP0.1`.\n",
    "\n",
    "Expand the pointer recursively to show the full contents of the selected collection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "collection = \"2.2i/runs/DP0.1\"\n",
    "print(collection)\n",
    "for c in sorted(registry.queryCollections(collection, flattenChains=True)):\n",
    "    print(c, registry.getCollectionType(c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Create a new Butler instance that specifies the `2.2i/runs/DP0.1` collection, and a new registry for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butler = dafButler.Butler(repo, collections=collection)\n",
    "registry = butler.registry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2 Butler DatasetType\n",
    "\n",
    "The LSST Science Pipelines classify data products as `DatasetTypes`.\n",
    "To demonstrate how to see the available `DatasetTypes`, the following cell prints them all to screen.\n",
    "\n",
    "As individual `DatasetTypes` are defined globally and do not belong to a specific collection, the following query returns *all* that belong to the repository, not just in the collection of interest. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in sorted(registry.queryDatasetTypes()):\n",
    "    print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "Here are some definitions to help delegates understand the contents of the full DP0.1 repository. \n",
    "\n",
    "* `calexp` - a single CCD of a processed visit image (PVI; individual calibrated exposures)\n",
    "* `deepCoadd` - products related to the coadded (stacked) images, including catalogs of coadd sources\n",
    "* `src` - a catalog of sources\n",
    "* `skyMap` - geometric representations of the sky coverage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Which data sets are most appropriate for DP0.1?** <br>\n",
    "> Most DP0.1 delegates will only be interested in data sets with types `ExposureF` or `SourceCatalog`. \n",
    "> For images, stick to the `calexp` (processed visit images, or PVIs) and `deepCoadd` (stacked PVIs).\n",
    "> For catalogs, the `src` should be used with the `calexp` images, and the `deepCoadd_forced_src` are the most appropriate to be used with the `coadds`.\n",
    "> More information can be found in the DP0.1 Data Products Definitions Document (DPDD) at [dp0-1.lsst.io](http://dp0-1.lsst.io)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3 Butler data IDs\n",
    "\n",
    "The data ID is a dictionary-like identifier for a data product.\n",
    "Find more about the data IDs [here](https://pipelines.lsst.io/v/weekly/modules/lsst.daf.butler/dimensions.html#data-ids).\n",
    "\n",
    "Each `DatasetType` uses a different set of keys in its data ID.\n",
    "For example, in the `DatasetType` list printed to screen (above), next to `calexp` in curly brackets is listed the band, instrument, detector, physical_filter, visit_system, and visit.\n",
    "These are the keys of the data ID for a `calexp`, which are also called \"dimensions\".\n",
    "\n",
    "In the following cell, the `DatasetRef` is queried for `calexp` data in our collection of interest, and the full data IDs are printed to screen (for just a few examples).\n",
    "Data IDs can be represented in code as regular Python `dict` objects, but when returned from the `Butler` the `DataCoordinate` class is used instead.\n",
    "\n",
    "The data ID contains both *implied* and *required* keys.\n",
    "For example, the value of *band* would be *implied* by the *visit*, because a single visit refers to a single exposure at a single pointing in a single band. \n",
    "In the following cell, printing the data ID without specifying `.full` shows only the required keys.\n",
    "The value of a single key, in this case *band*, can also be printed by specifying the key name.\n",
    "\n",
    "The following cell will fail and return an error if the query is requesting a `DatasetRef` for data that does not exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasetRefs = registry.queryDatasets(datasetType='calexp', collections=collection)\n",
    "\n",
    "for i, ref in enumerate(datasetRefs):\n",
    "    print(ref.dataId.full)\n",
    "    print(ref.dataId)\n",
    "    print(ref.dataId['band'])\n",
    "    print(' ')\n",
    "    if i > 2:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Uniform Resource Identifier (URI) is the closest thing to a \"filepath to the data\" that the Butler provides. \n",
    "Note that this URI does not refer to a local path on the filesystem.\n",
    "There is no need to know exactly where the data live in order to access it - that's the power of the Butler!!\n",
    "\n",
    "The following cell prints the URI to screen as a demonstration of an alternate way to uniquely identify data in the Butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, ref in enumerate(datasetRefs):\n",
    "    print('File URI: ', butler.getURI(ref))\n",
    "    if i > 2:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4 Butler queryDatasets\n",
    "\n",
    "Above demonstrated a very simple use of `queryDatasets`, but additional query terms can also be used, such as band and visit.\n",
    "When a query term is an equality, it can be specified like `band='g'`. \n",
    "When a query term is an inequality, it can be specified with `where`.\n",
    "Details about Butler queries can be found [here](https://pipelines.lsst.io/v/weekly/modules/lsst.daf.butler/queries.html).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasetRefs = registry.queryDatasets(datasetType='calexp', band='g',\n",
    "                                     where='visit > 700000', collections=collection)\n",
    "\n",
    "for i, ref in enumerate(datasetRefs):\n",
    "    print(ref.dataId.full)\n",
    "    if i > 2:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Each data ID key-value pair is associated with a metadata row called a `DimensionRecord`.\n",
    "Like dataset types, these exist independent of any collection, but they are also identified by data IDs.\n",
    "\n",
    "The `queryDimensionsRecords` method provides a way to query for these records.\n",
    "Most of the arguments accepted by `queryDatasets` can be used here (including `where`).\n",
    "\n",
    "An example of this is provided below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for dim in ['exposure', 'visit', 'detector']:\n",
    "    print(list(registry.queryDimensionRecords(dim, where='visit = 971990 and detector=0'))[0])\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another query method, `queryDataIds`, can be used to query for data IDs independent of any dataset, but it's less useful for general data exploration.\n",
    "\n",
    "It is also possible to pass `datasets` and `collections` to both `queryDataIds` and `queryDimensionRecords` in order to return records whose data IDs match those of existing datasets.\n",
    "But this is quite a bit more subtle than searching directly for a dataset, and rarely wanted when exploring a data repository.\n",
    "\n",
    "More information on all of the query methods can be found [here](https://pipelines.lsst.io/v/weekly/middleware/faq.html#when-should-i-use-each-of-the-query-methods-commands)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5 Temporal and spatial queries\n",
    "\n",
    "The following examples show how to query for data sets that include a desired coordinate and observation date.\n",
    "\n",
    "##### Temporal queries\n",
    "\n",
    "Above, we can see that for visit 971990, the (RA,Dec) are (70.37770,-37.1757) and the observation date is 20251201.\n",
    "But these are just human-readable summaries of the more precise spatial and temporal information stored in the registry, which are represented in Python by `Timespan` and `Region` objects, respectively.\n",
    "`DimensionRecord` objects that represent spatial or temporal concepts (a `visit` is both) have these objects attached to them.\n",
    "\n",
    "Retrieve the `DimensionRecord` for a visit and show its timespan and region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(record,) = registry.queryDimensionRecords('visit', visit=971990)\n",
    "\n",
    "print(record.timespan)\n",
    "print(' ')\n",
    "print(record.region)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the timespan or spatial region that are being used as query constraints are already associated with a data ID in the database, the spatial and temporal overlap constraints are automatic.\n",
    "For example, if we query for `deepCoadd` datasets with a `visit`+`detector` data ID, we'll get just the ones that overlap that observation and have the same band (because a visit implies a band):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ref in registry.queryDatasets(\"deepCoadd\", visit=971990, detector=50):\n",
    "    print(ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To query for dimension records or datasets that overlap an arbitrary time range, we can use the `bind` argument to pass times through to `where`.\n",
    "Using `bind` to define an alias for a variable saves us from having to string-format the times into the `where` expression.\n",
    "Note that a `dafButler.Timespan` will accept a `begin` or `end` value that is equal to `None` if it is unbounded on that side.\n",
    "\n",
    "Use `bind` and `where`, along with [astropy.time](https://docs.astropy.org/en/stable/time/index.html), to look for visits within one minute of this one on either side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import astropy.time\n",
    "# minute = astropy.time.TimeDelta(60, format=\"sec\")\n",
    "# timespan = dafButler.Timespan(record.timespan.begin - minute, record.timespan.end + minute)\n",
    "\n",
    "# for visit in registry.queryDimensionRecords(\"visit\", where=\"visit.timespan OVERLAPS my_timespan\", \n",
    "#                                             bind={\"my_timespan\": timespan}):\n",
    "#     print(visit.id, visit.timespan, visit.physical_filter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.time\n",
    "minute = astropy.time.TimeDelta(60, format=\"sec\")\n",
    "timespan = dafButler.Timespan(record.timespan.begin - minute, record.timespan.end + minute)\n",
    "\n",
    "datasetRefs = registry.queryDatasets(\"calexp\", where=\"visit.timespan OVERLAPS my_timespan\",\n",
    "                                     bind={\"my_timespan\": timespan})\n",
    "\n",
    "for i, ref in enumerate(datasetRefs):\n",
    "    print(ref)\n",
    "    if i > 6:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Spatial queries\n",
    "\n",
    "Arbitrary spatial queries are not supported at this time, such as the \"POINT() IN (REGION)\" example found in this [Butler queries](https://pipelines.lsst.io/v/weekly/modules/lsst.daf.butler/queries.html) documentation.\n",
    "In other words, at this time it is only possible to do queries involving regions that are already \"in\" the data repository, either because they are HTM pixel regions or because they are tract/patch/visit/visit+detector regions.\n",
    "\n",
    "Thus, for this example we use the set of dimensions that correspond to different levels of the HTM (hierarchical triangular mesh) pixelization of the sky ([HTM primer](http://www.skyserver.org/htm/)).\n",
    "The process is to transform a region or point into one or more HTM identifiers (HTM IDs), and then create a query using the HTM ID as the spatial data ID.\n",
    "The `lsst.sphgeom` library supports region objects and HTM pixelization in the LSST Science Pipelines.\n",
    "\n",
    "Import the `lsst.sphgeom` package, initialize a sky pixelization to level 10 (the level at which one sky pixel is about five arcmin across), and find the HTM ID for a desired sky coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lsst.sphgeom\n",
    "\n",
    "pixelization = lsst.sphgeom.HtmPixelization(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "htm_id = pixelization.index(\n",
    "    lsst.sphgeom.UnitVector3d(\n",
    "        lsst.sphgeom.LonLat.fromDegrees(70.376995, -37.175736)\n",
    "    )\n",
    ")\n",
    "\n",
    "# Obtain and print the scale to provide a sense of the size of the sky pixelization being used\n",
    "scale = pixelization.triangle(htm_id).getBoundingCircle().getOpeningAngle().asDegrees()*3600\n",
    "print(f'HTM ID={htm_id} at level={pixelization.getLevel()} is a ~{scale:0.2}\" triangle.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datasetRefs = registry.queryDatasets(\"calexp\", htm20=htm_id,\n",
    "                                     where=\"visit.timespan OVERLAPS my_timespan\",\n",
    "                                     bind={\"my_timespan\": timespan})\n",
    "\n",
    "for i, ref in enumerate(datasetRefs):\n",
    "    print(ref)\n",
    "    if i > 6:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, with the above query, we have uniquely identified the visit and detector for our desired temporal and spatial constraints.\n",
    "\n",
    "Note that if a smaller HTM level is used (like 7), which is a larger sky pixel (~2200 arcseconds), the above query will return many more visits and detectors which overlap with that larger region. Try it and see!\n",
    "\n",
    "Note that queries using the HTM ID can also be used to, e.g., find the set of all i-band `src` catalog data products that overlap this point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, src_ref in enumerate(registry.queryDatasets(\"src\", htm20=htm_id, band=\"i\")):\n",
    "    print(src_ref)\n",
    "    if i > 2:\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is does that search take tens of seconds?\n",
    "The butler's spatial reasoning is designed to work well for regions the size of full data products, like detector- or patch-level images and catalogs, and it's a poor choice for object-scale searches.\n",
    "The above search is slow in part because `queryDatasets` searches for all `src` datasets that overlap a larger region and then filters the results down to the specified HTM ID pixel.\n",
    "\n",
    "Options for exploring and retrieving catalog data with the Butler is covered in more depth in Section 5."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Retrieve and plot an image with sources\n",
    "\n",
    "At this point, we have all the information we need to ask the Butler to get a specific data product.\n",
    "In the following example, the visit, detector, and band are used to define the `dataId` and the corresponding `calexp` is retrieved by the Butler.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId = {'visit': '703697', 'detector': 80, 'band': 'g'}\n",
    "calexp = butler.get('calexp', dataId=dataId)\n",
    "\n",
    "# This will print a warning related to the gen2 to gen3 Butler conversion.\n",
    "# It is ok to ignore this warning for DP0.1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the `calexp` is a calibrated CCD from a single exposure.\n",
    "Use the afwDisplay interface to show the pixel values and mask plane.\n",
    "\n",
    "The blue coloring and red streaks seen in the image below is set by the \"mask\" plane of the `calexp`.\n",
    "The mask encodes information such as bad or saturated pixels.\n",
    "In this case, blue indicates a detected source.\n",
    "See the Image Display and Manipulation tutorial for more about afwDisplay and the mask plane. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "display = afwDisplay.Display()\n",
    "display.scale('linear', 'zscale')\n",
    "display.mtv(calexp)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "To retrieve the catalog of sources detected in this `calexp`, pass the `dataId` to the Butler and request the `src` `datasetType`. \n",
    "This performs another query on the registry database.\n",
    "See Section 5 for more examples of querying the Butler for catalog data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = butler.get('src', dataId)\n",
    "\n",
    "# Define src as a copy in order to manipulate it.\n",
    "src = src.copy(True)\n",
    "\n",
    "# If desired, show src as an AstroPy table (nicely formatted).\n",
    "# src.asAstropy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the `calexp` with the `src` catalog overlaid.\n",
    "An investigation of this image is left as an exercise to the user! :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "afw_display = afwDisplay.Display(1)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(calexp)\n",
    "plt.gca().axis('off')\n",
    "\n",
    "# Use display buffering to avoid re-drawing the image after each source is plotted\n",
    "with afw_display.Buffering():\n",
    "    for s in src:\n",
    "        afw_display.dot('+', s.getX(), s.getY(), ctype=afwDisplay.RED)\n",
    "        afw_display.dot('o', s.getX(), s.getY(), size=20, ctype='orange')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Create an image cutout at a user-specified coordinate\n",
    "\n",
    "In this section, the LSST Science Pipelines geometry (recall \"import lsst.geom as geom\") and coordinate (\"import lsst.afw.coord as afwCoord\") packages are used to define a function which returns a cutout at a defined location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cutout_coadd(butler, ra, dec, band='r', datasetType='deepCoadd',\n",
    "                 skymap=None, cutoutSideLength=51, **kwargs):\n",
    "    \"\"\"\n",
    "    Produce a cutout from a coadd at the given afw SpherePoint radec position.\n",
    "\n",
    "    Adapted from DC2 tutorial notebook by Michael Wood-Vasey.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    butler: lsst.daf.persistence.Butler\n",
    "        Client providing access to a data repository\n",
    "    ra: float\n",
    "        Right ascension of the center of the cutout, degrees\n",
    "    dec: float\n",
    "        Declination of the center of the cutout, degrees\n",
    "    filter: string\n",
    "        Filter of the image to load\n",
    "    datasetType: string ['deepCoadd']\n",
    "        Which type of coadd to load.  Doesn't support 'calexp'\n",
    "    skymap: lsst.afw.skyMap.SkyMap [optional]\n",
    "        Pass in to avoid the Butler read.  Useful if you have lots of them.\n",
    "    cutoutSideLength: float [optional]\n",
    "        Side of the cutout region in pixels.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    MaskedImage\n",
    "    \"\"\"\n",
    "\n",
    "    radec = geom.SpherePoint(ra, dec, geom.degrees)\n",
    "    cutoutSize = geom.ExtentI(cutoutSideLength, cutoutSideLength)\n",
    "\n",
    "    if skymap is None:\n",
    "        skymap = butler.get(\"skyMap\")\n",
    "\n",
    "    # Look up the tract, patch for the RA, Dec\n",
    "    tractInfo = skymap.findTract(radec)\n",
    "    patchInfo = tractInfo.findPatch(radec)\n",
    "    xy = geom.PointI(tractInfo.getWcs().skyToPixel(radec))\n",
    "    bbox = geom.BoxI(xy - cutoutSize//2, cutoutSize)\n",
    "    patch = tractInfo.getSequentialPatchIndex(patchInfo)\n",
    "\n",
    "    coaddId = {'tract': tractInfo.getId(), 'patch': patch, 'band': band}\n",
    "    parameters = {'bbox':bbox}\n",
    "\n",
    "    cutout_image = butler.get(datasetType, parameters=parameters,\n",
    "                              immediate=True, dataId=coaddId)\n",
    "\n",
    "    return cutout_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the center of the DC2 region as an example\n",
    "ra, dec = 55.064, -29.783\n",
    "cutout_image = cutout_coadd(butler, ra, dec, datasetType='deepCoadd', \n",
    "                            cutoutSideLength=201)\n",
    "print(\"The size of the cutout is: \", cutout_image.image.array.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the image cutout\n",
    "fig = plt.figure()\n",
    "afw_display = afwDisplay.Display(1)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(cutout_image.image)\n",
    "plt.gca().axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Exploring and retrieving catalog data from the Butler\n",
    "\n",
    "The TAP service is the recommended way to retrieve DP0.1 catalog data for a notebook, and there are several other [tutuorials](https://github.com/rubin-dp0/tutorial-notebooks) that demonstrate how to use the TAP service.\n",
    "\n",
    "However, as we saw above, the Butler can also be used to access to catalog data. We can investigate the table schema for a specific source catalog by  Butler appending `_schema` to the `datasetType`. Note that this does not require you to specify the ``dataId``. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schema_coadd_src = butler.get('deepCoadd_forced_src_schema')\n",
    "schema_coadd_src.asAstropy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The table `schema` stores information about the columns stored in the table. Each of the following lines will print the schema to the screen in different ways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# schema_coadd_src.schema\n",
    "# schema_coadd_src.schema.getNames()\n",
    "# schema_coadd_src.schema.getOrderedNames()\n",
    "print('Number of columns in this table = ', len(schema_coadd_src.schema.getNames()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps you want to search for all schema elements that contain the term 'psf'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define an array that is all of the column names\n",
    "all_names = schema_coadd_src.schema.getOrderedNames()\n",
    "\n",
    "# Loop over the names and look for the term 'psf'\n",
    "for i, name in enumerate(all_names):\n",
    "    if name.find('psf') >= 0:\n",
    "        print(i, name)\n",
    "del all_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Probably you will want to know more about the values in these columns. You can do that by printing the documentation string in the schema."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn the schema into a python dictionary, to be able to call a column by name\n",
    "schema_dict = schema_coadd_src.schema.extract('*')\n",
    "\n",
    "# Print the associated docstring for each of the named columns of interest\n",
    "for name in ['base_SdssShape_psf_xx', 'base_SdssShape_psf_yy', 'base_SdssShape_psf_xy']:\n",
    "    doc = schema_dict[name].getField().getDoc()\n",
    "    units = schema_dict[name].getField().getUnits()\n",
    "    print(name, '[%s]'%units, ' = ', doc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Refer to the DP0.1 Data Products Definitions Document (DPDD) at [dp0-1.lsst.io](http://dp0-1.lsst.io) to find out more about the columns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "The full catalogs are very large and it is not feasible to try and retrieve them in their entirety.\n",
    "Instead, in this example we identify the tract and patch of interest and retrieve only catalog data for a small region of sky.\n",
    "We use the same ra and dec coordinates as above to find the patch and tract."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "radec = geom.SpherePoint(ra, dec, geom.degrees)\n",
    "\n",
    "skymap = butler.get(\"skyMap\")\n",
    "\n",
    "tractInfo = skymap.findTract(radec)\n",
    "tract = tractInfo.getId()\n",
    "\n",
    "patchInfo = tractInfo.findPatch(radec)\n",
    "patch = tractInfo.getSequentialPatchIndex(patchInfo)\n",
    "\n",
    "print(tract, patch)\n",
    "\n",
    "coaddId = {'tract': tract, 'patch': patch, 'band': 'i'}\n",
    "\n",
    "coadd_src = butler.get('deepCoadd_forced_src', coaddId)\n",
    "coadd_src = coadd_src.copy(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the table contents if desired\n",
    "# coadd_src.asAstropy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Convert to a Pandas dataframe (see the first tutorial) for easy interaction.\n",
    "The following cells offer options for printing the column names or the data values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = coadd_src.asAstropy().to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(data.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for col in data.columns:\n",
    "#     print(col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data['coord_ra'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the locations of sources in the patch as well as the ra,dec that we requested. Note that the `coord_ra` and `coord_dec` are in radians, so we need to convert them to degrees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.plot(np.degrees(data['coord_ra'].values),\n",
    "         np.degrees(data['coord_dec'].values),\n",
    "         'o', ms=2, alpha=0.5)\n",
    "plt.plot(ra, dec, '*', ms=25, mec='k')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Dec (deg)')\n",
    "plt.title('Butler coadd_forced_src objects in tract 4638 patch 43')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we noted, the `coord_ra` and `coord_dec` columns have units of _radians_. As an exercise, you could use what you've learned from above to confirm this by accessing the table schema. (Also note that you can scroll up and find the answer in the outputs from a cell you already executed.) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}