Adding cutout section

Author: kadrlica

04_Intro_to_Butler.ipynb

@@ -15,13 +15,14 @@
     "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>\n",
     "\n",
-    "This notebook provides an introduction to the use of the Gen-3 Butler to explore the structure and content of the DP0.1 data repository containing the LSST DESC DC2 data. The Butler is the LSST Science Pipelines interface for managing, reading, and writing datasets.  The Gen-3 Butler is still under development, so this notebook is expected to be updated after the official Gen-3 release.\n",
+    "This notebook provides an introduction to the use of the data Butler. 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.\n",
     "\n",
     "The goals of this notebook are:\n",
     "\n",
-    "1. Create a Butler\n",
+    "1. Create a Butler.\n",
     "2. Programmatically explore the content of the DP0.1 data repository\n",
-    "3. Grab some data!\n",
+    "3. Grab some DP0.1 data! \n",
+    "4. Grab a cutout of the DP0.1 coadds at a specific (RA,Dec).\n",
     "\n",
     "**Credit:** This notebook was originally developed by Alex Drlica-Wagner in the context of the LSST Stack Club.\n",
     "\n",
@@ -44,9 +45,7 @@
    "source": [
     "### 1. Create the Butler\n",
     "\n",
-    "In this section we are going to create a Gen-3 Butler by pointing it to a specific data set.\n",
-    "\n",
-    "We start with some general python package imports."
+    "In this section we are going to create an instance of the Butler. We start with some general python package imports."
    ]
   },
   {
@@ -65,7 +64,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We import LSST Science Pipelines packages to create the Butler and display some data. More details on data display can be found in [03_Intro_to_AFW_Display.ipynb](03_Intro_to_AFW_Display.ipynb)."
+    "We import several packages from the LSST Science Pipelines. The first import gives us access to the Butler, while the second provides tools for displaying data. More details on data display can be found in [03_Intro_to_AFW_Display.ipynb](03_Intro_to_AFW_Display.ipynb)."
    ]
   },
   {
@@ -84,7 +83,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "To create the Butler, we need give it a path to the data, which is called a \"repository\". Butler repositories can be remote (the data is on a server, across a network) or local (the data in on a local filesystem)."
+    "To create the Butler, we need to provide it with a path to the data set, which is called a \"data repository\". Butler repositories can be remote (i.e., pointing to an S3 bucket) or local (i.e., pointing to a directory on the local file system). In this case, we point to an S3 bucket."
    ]
   },
   {
@@ -103,7 +102,7 @@
    "source": [
     "### 2. Explore a data repository\n",
     "\n",
-    "With the butler created, we can now access the data `registry` (a database containing information about available data products). The registry will help us examine what collections of data products exist."
+    "Now that we've created an instance of the Butler, we can access the data `registry` (a database containing information about available data products). The registry will help us examine what collections of data products exist."
    ]
   },
   {
@@ -448,12 +447,103 @@
     "    print()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 5. Grabbing an Image Cutout\n",
+    "\n",
+    "Say we want to grab a cutout of the DP0.1 coadded images at a specific location. In order to do this, we need a few other packages from the LSST Science Pipelines. In particular, access to the geometry and coordinate packages."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import lsst.geom as geom\n",
+    "import lsst.afw.coord as afwCoord"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "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",
+    "        Servant 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",
+    "    radec = geom.SpherePoint(ra, dec, afwGeom.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, immediate=True, dataId=coaddId)\n",
+    "    \n",
+    "    return cutout_image"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ra,dec = 55.064, -29.783 # center of DC2 region\n",
+    "cutout_image = cutout_coadd(butler,ra,dec,datasetType='deepCoadd', cutoutSideLength=201)\n",
+    "print(\"The size of the cutout is: \",cutout_image.image.array.shape)"
+   ]
+  },
+  {
+   "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(cutout_image.image)\n",
+    "plt.gca().axis('off')"
+   ]
   }
  ],
  "metadata": {