DP02_12a_PSF_Data_Products.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "# Introduction to Point Spread Function (PSF) Data Products\n",
    "\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<br>\n",
    "Contact author(s): Andrés A. Plazas Malagón<br>\n",
    "Last verified to run: 2024-09-16 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_37 <br>\n",
    "Container size: medium <br>\n",
    "Targeted learning level: intermediate <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** A demonstration of how to access `calexp` and `deepCoadd` PSF properties."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Use of single-epoch and coadded PSF models."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** DP0.2 `calexp` and `deepCoadd` images."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:**  lsst.daf.butler, lsst.geom, lsst.afw.display"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credits:** Developed by Andrés A. Plazas Malagón in collaboration with Melissa Graham, Jeff Carlin, Ryan Lau, and the Rubin Community Science Team for DP0.2. The functions for studying the PSF image profile are built upon the rapid-analysis code originally created by Merlin Fisher-Levine to characterize the Point Spread Function (PSF) of the Rubin Auxiliary Telescope LSST Atmospheric Transmission and Slitless Spectrograph (LATISS) images.\n",
    "This notebook incorporates suggestions from the [Accessible Authoring Checklist](https://iota-school.github.io/accessibility_hackathon/hack#checklist) and utilizes NASA's Astrophysics Data System Bibliographic Services."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28e91cbf-ab7f-4e26-9276-b00299d6065e",
   "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": {
    "tags": []
   },
   "source": [
    "## 1. Introduction\n",
    "\n",
    "The Point Spread Function (PSF) can be understood as a function that describes how a bundle of rays, initially converging on a single point, spreads out spatially around that point (in simplest terms, the PSF can be thought of as how light from an astronomical point source is \"spread out\" in an image).\n",
    "\n",
    "To accurately characterize the PSF, it is common to analyze the observed images of stars, which act as point sources before being distorted by the PSF.\n",
    "This analysis helps determine the convolution kernel, representing the size and shape of the blurring effect caused by the PSF.\n",
    "Accurate PSF modeling is crucial because any inaccuracies in understanding these effects can lead to erroneous conclusions about fundamental aspects of the universe (e.g., the properties of dark matter and dark energy).\n",
    "Therefore, understanding and characterizing the PSF is essential to properly interpret and extract reliable information from astronomical observations.\n",
    "\n",
    "This tutorial demonstrates how the PSF is measured by the LSST Science Pipelines for all images,\n",
    "with a focus on `calexp` and `deepCoadd` images.\n",
    "It provides an overview of the PSF-related information and data products that are available for LSST images,\n",
    "and shows how to calculate PSF profiles and contours, along with other properties such as size.\n",
    "\n",
    "**How are PSFs measured by the LSST Science Pipelines?**\n",
    "\n",
    "1. The initial stage in processing LSST science observations involves Instrument Signature Removal (ISR).\n",
    "This encompasses fundamental detrending operations like flat-fielding, bias subtraction, fringe correction,\n",
    "and rectification of flawed and oversaturated pixels. \n",
    "\n",
    "2. Subsequently, the focus shifts to characterizing single-epoch direct images, constructing models that \n",
    "depict the observational system and its conversion of the true celestial scene into the observed image.\n",
    "This process encompasses background subtraction, PSF modeling, addressing cosmic ray effects, \n",
    "applying aperture corrections, and source measurement. \n",
    "\n",
    "3. After segregating stars and galaxies, the PSF model in single-epoch images is assembled using data from a star catalog. \n",
    "\n",
    "The LSST Science Pipelines have partially transitioned to using the [PIFF (PSF in the Full-Field of View)](https://rmjarvis.github.io/Piff/_build/html/index.html) code for PSF modeling in the full field of view. However, the DP0.2 PSF models were generated using the Point Spread Function Extractor, [PSFEx](https://ui.adsabs.harvard.edu/abs/2013ascl.soft01001B/abstract), software.\n",
    "Understanding or using these codes is not necessary for completing this tutorial, nor for doing PSF\n",
    "analysis with the LSST images."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1.1. Import packages\n",
    "\n",
    "The [matplotlib](https://matplotlib.org/) (and especially sublibrary `matplotlib.pyplot`), \n",
    "[numpy](http://www.numpy.org/), and [scipy](https://scipy.org/) libraries are widely used \n",
    "Python libraries for plotting, scientific computing, and scientific analysis.\n",
    "The `getpass` package includes functionality to programmatically retrieve user account information,\n",
    "such as username.\n",
    "\n",
    "\n",
    "The `lsst.daf.butler` library is used to access data products via the butler, \n",
    "and the `lsst.afw.display` library provides access to image visualization routines.\n",
    "The `lsst.geom` library provides the representation of a 2D coordinate `Point2D`.\n",
    "\n",
    "<!--- and [astropy](http://www.astropy.org/) --->\n",
    "\n",
    "<!--- The `lsst.rsp` package provides access to the Table Access Protocol (TAP) service for queries to the DP0 catalogs. --->\n",
    "\n",
    "<!--- The `lsst.afw.image` provide visualization tools. --->\n",
    "\n",
    "<!--- The `lsst.afw.analysis.tools` library provides acces to analysis tools.--->\n",
    "\n",
    "<!--- The `lsst.pipe.base` library provides access to the `lsst.pipe.base.Struct` container. --->\n",
    "\n",
    "<!--- [treecorr](https://rmjarvis.github.io/TreeCorr/_build/html/index.html) by [Jarvis et al. 2004](https://ui.adsabs.harvard.edu/abs/2004MNRAS.352..338J/abstract) is a widely used tool for fast correlations measurements based on a ball tree method (similar to a [k-d tree](https://en.wikipedia.org/wiki/K-d_tree)). In particular, it is widely used in state-of-the-art weak lensing cosmological analyses. --->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d185358a-4c15-4b04-8e3d-13083191ee72",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.stats import skew\n",
    "import getpass\n",
    "\n",
    "import lsst.daf.butler as dafButler\n",
    "import lsst.afw.display as afwDisplay\n",
    "from lsst.geom import Point2D, radToDeg, SpherePoint, degrees"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa240aaf-d675-4d72-bb68-d736200f4980",
   "metadata": {},
   "source": [
    "### 1.2. Define parameters and functions\n",
    "\n",
    "The following functions and parameters are defined once, and used throughout this notebook.\n",
    "\n",
    "> **Notice: It is not necessary to read and understand all functions!**\n",
    "\n",
    "It is OK to simply execute all of the cells in Section 1.2 and move on to Section 2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb203afe-4448-4a6f-add3-193c1e824302",
   "metadata": {},
   "source": [
    "Define the Butler instance for DP0.2 in order to access images with this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb1614d4-3280-4580-8900-ade50aea6132",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "butler = dafButler.Butler('dp02', collections='2.2i/runs/DP0.2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8dd7f04-f3d2-4f92-9927-acb9f47fe772",
   "metadata": {},
   "source": [
    "The following cell will set a standard figure size and `afwDisplay` backend to use throughout the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59f35be0-b960-46d7-bc1e-b875f641fd7a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8.0, 8.0)\n",
    "afwDisplay.setDefaultBackend('matplotlib')\n",
    "plt.style.use('tableau-colorblind10')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0781578e-6b82-4820-9fbb-bbee1c432e5e",
   "metadata": {},
   "source": [
    "**Constant**: `SIGMA_TO_FWHM` \n",
    "\n",
    "Defines the convertion factor between the standard deviation `sigma` of a Gaussian function a the full-withd at half-maximum (FWHM): $2\\sqrt{2\\ln(2)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d654417-a7d2-4a2c-9846-de3ee6bcbafb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "SIGMA_TO_FWHM = 2.0*np.sqrt(2.0*np.log(2.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a58cc22a-01ff-4b01-9b62-e2fb5ee08618",
   "metadata": {},
   "source": [
    "**Function**: `gauss`\n",
    "\n",
    "Defines a one-dimensional Gaussian profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f3386e1-b0e7-49a1-b3e3-a90fafd95163",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def gauss(x, a, x0, sigma):\n",
    "    return a*np.exp(-(x-x0)**2/(2*sigma**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "182df6b2-16ef-4bed-b91f-202334a73f15",
   "metadata": {},
   "source": [
    "**Function**: `getPsfProperties`\n",
    "\n",
    "Given a PSF model (`lsst.meas.extensions.psfex.PsfexPsf`) and a coordinate where the model is being evaluated\n",
    "(`lsst.geom.Point2D`) this functions returns the PSF FWHM, flux from aperture photometry, \n",
    "peak value of the normalized PSF (the PSF is normalized to a sum of 1), and size of the PSF postage stamp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25957b72-15ba-4459-a160-27c77b68d483",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def getPsfProperties(psf, point):\n",
    "    \"\"\"Function to obtain PSF properties.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    psf : `lsst.meas.extensions.psfex.PsfexPsf`\n",
    "        PSF object.\n",
    "    point : `lsst.geom.Point2D`\n",
    "        Coordinate where the PSF is being evaluated.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    fwhm : `float`\n",
    "        Full-width at half maximum: PSF determinant radius\n",
    "        from SDSS adaptive moments matrix (sigma) times\n",
    "        SIGMA_TO_FWHM.\n",
    "    ap_flux : `float`\n",
    "        PSF flux from aperture photometry weighted\n",
    "        by a sinc function.\n",
    "    peak : `float`\n",
    "        Peak PSF value.\n",
    "    dims : `lsst.geom.ExtendI`\n",
    "        PSF postage stamp dimensions.\n",
    "    \"\"\"\n",
    "    sigma = psf.computeShape(point).getDeterminantRadius()\n",
    "    fwhm = sigma * SIGMA_TO_FWHM\n",
    "    ap_flux = psf.computeApertureFlux(radius=sigma, position=point)\n",
    "    peak = psf.computePeak(position=point)\n",
    "    dims = psf.computeImage(point).getDimensions()\n",
    "\n",
    "    print(f\"PSF FWHM: {fwhm:.4} pix \\n\"\n",
    "          f\"PSF flux from aperture photometry: {ap_flux:.4} \\n\"\n",
    "          f\"Peak PSF value: {peak:.4} \\n\"\n",
    "          f\"PSF postage stamp dimensions: {dims} \\n\")\n",
    "\n",
    "    return (sigma, ap_flux, peak, dims)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9376ef28-0cb2-44aa-ae2c-d569c44f1966",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-16T23:31:51.836818Z",
     "iopub.status.busy": "2023-09-16T23:31:51.836052Z",
     "iopub.status.idle": "2023-09-16T23:31:52.051703Z",
     "shell.execute_reply": "2023-09-16T23:31:52.050845Z",
     "shell.execute_reply.started": "2023-09-16T23:31:51.836792Z"
    }
   },
   "source": [
    "**Function**: `plotRadialAverage`\n",
    "\n",
    "Given an image of the PSF (`lsst.afw.image.ExposureF`), this function \n",
    "plots the azimuthally-averaged radial profile of the PSF, and fits a Gaussian function to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd51deb0-4d72-405c-a561-090754b2ccbc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plotRadialAverage(exp, ax=None):\n",
    "    \"\"\"\n",
    "    Function to plot the radial average of a point spread function (PSF) image.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    exp : `lsst.afw.image.ExposureF`\n",
    "        The PSF image exposure.\n",
    "    ax : `matplotlib.axes.Axes`, optional\n",
    "        If provided, the plot will be drawn on this axis (default is None).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    ax : `matplotlib.axes.Axes`\n",
    "        The axis on which the plot was drawn.\n",
    "    \"\"\"\n",
    "    data = exp.array\n",
    "    xlen, ylen = data.shape\n",
    "    center = np.array([xlen / 2, ylen / 2])\n",
    "    distances = []\n",
    "    values = []\n",
    "\n",
    "    for i in range(xlen):\n",
    "        for j in range(ylen):\n",
    "            value = data[i, j]\n",
    "            dist = np.linalg.norm(np.array([i, j]) - center)\n",
    "            if dist <= xlen / 2:\n",
    "                values.append(value)\n",
    "                distances.append(dist)\n",
    "\n",
    "    peakPos = 0\n",
    "    amplitude = np.max(values)\n",
    "    width = 10\n",
    "\n",
    "    bounds = ((0, 0, 0), (np.inf, np.inf, np.inf))\n",
    "\n",
    "    try:\n",
    "        pars, pCov = curve_fit(gauss, distances, values,\n",
    "                               [amplitude, peakPos, width], bounds=bounds)\n",
    "        pars[0] = np.abs(pars[0])\n",
    "        pars[2] = np.abs(pars[2])\n",
    "    except RuntimeError:\n",
    "        pars = None\n",
    "\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    ax.plot(distances, values, 'x', label='Radial average')\n",
    "\n",
    "    if pars is not None:\n",
    "        fitAmp = pars[0]\n",
    "        fitGausMean = pars[1]\n",
    "        fitFwhm = pars[2] * SIGMA_TO_FWHM\n",
    "        fitline = gauss(distances, *pars)\n",
    "\n",
    "        x_fit = np.linspace(-np.max(distances), np.max(distances), 100)\n",
    "        y_fit = gauss(x_fit, *pars)\n",
    "        skewness = skew(y_fit, axis=0, bias=True)\n",
    "\n",
    "        ax.plot(distances, fitline,\n",
    "                label=f\" Gaussian Fit \\n Amp: {fitAmp:.3f}\"\n",
    "                      f\"\\n Position: {fitGausMean:.2f}\"\n",
    "                      f\"\\n FWHM: {fitFwhm:.2f}\"\n",
    "                      f\"\\n \\n Skewness: {skewness:.2f}\")\n",
    "\n",
    "    ax.set_ylabel('Flux (ADU)')\n",
    "    ax.set_xlabel('Radius (pix)')\n",
    "    ax.set_aspect(1.0 / ax.get_data_ratio(), adjustable='box')\n",
    "    ax.legend()\n",
    "    ax.set_title(\"Azimuthally-averaged radial profile\")\n",
    "\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fde6adc9-3a90-45c2-9b12-da811d75953a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-18T14:08:48.613653Z",
     "iopub.status.busy": "2023-09-18T14:08:48.612698Z",
     "iopub.status.idle": "2023-09-18T14:08:49.548203Z",
     "shell.execute_reply": "2023-09-18T14:08:49.546510Z",
     "shell.execute_reply.started": "2023-09-18T14:08:48.613621Z"
    }
   },
   "source": [
    "**Function**: `plotCurveOfGrowth`\n",
    "\n",
    "Given an image of the PSF (`lsst.afw.image.ExposureF`), this function\n",
    "plots the encircled energy of the PSF as a function of radius (curve of growth)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17d482c9-b4c5-4151-bf1b-3c7162f6ca32",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plotCurveOfGrowth(exp, ax=None):\n",
    "    \"\"\"\n",
    "    Function to plot the curve of growth of a point\n",
    "    spread function (PSF) image.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    exp : `lsst.afw.image.ExposureF`\n",
    "        The PSF image exposure.\n",
    "    ax : `matplotlib.axes.Axes`, optional\n",
    "        If provided, the plot will be drawn on this axis (default is None).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    ax : `matplotlib.axes.Axes`\n",
    "        The axis on which the plot was drawn.\n",
    "    \"\"\"\n",
    "    data = exp.array\n",
    "    xlen, ylen = data.shape\n",
    "    center = np.array([xlen / 2, ylen / 2])\n",
    "    distances = []\n",
    "    values = []\n",
    "\n",
    "    for i in range(xlen):\n",
    "        for j in range(ylen):\n",
    "            value = data[i, j]\n",
    "            dist = np.linalg.norm(np.array([i, j]) - center)\n",
    "            if dist <= xlen / 2:\n",
    "                values.append(value)\n",
    "                distances.append(dist)\n",
    "\n",
    "    d = np.array([(r, v) for (r, v) in sorted(zip(distances, values))])\n",
    "    radii = d[:, 0]\n",
    "    values = d[:, 1]\n",
    "    cum_fluxes = np.cumsum(values)\n",
    "    cum_fluxes_norm = cum_fluxes / np.max(cum_fluxes)\n",
    "\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    ax.plot(radii, cum_fluxes_norm, markersize=10)\n",
    "    ax.set_ylabel('Encircled flux (%)')\n",
    "    ax.set_xlabel('Radius (pix)')\n",
    "    ax.set_title(\"Encircled flux\")\n",
    "\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8e6b637-76ca-4079-a43c-35f6def23dbd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-16T22:28:53.117282Z",
     "iopub.status.busy": "2023-09-16T22:28:53.116494Z",
     "iopub.status.idle": "2023-09-16T22:28:53.666615Z",
     "shell.execute_reply": "2023-09-16T22:28:53.665738Z",
     "shell.execute_reply.started": "2023-09-16T22:28:53.117257Z"
    }
   },
   "source": [
    "**Function**: `plotContours` \n",
    "\n",
    "Given an image of the PSF (`lsst.afw.image.ExposureF`) and, optionally,\n",
    "the plot axis and number of contours (`nContours`), this function\n",
    "plots the two-dimensional contous of the PSF postage stamp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f73d9bea-80f2-442d-99ac-7e23783304fd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plotContours(exp, ax=None, nContours=10):\n",
    "    \"\"\"\n",
    "    Function to plot contour lines of a point spread function (PSF) image.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    exp : `lsst.afw.image.ExposureF`\n",
    "        The PSF image exposure.\n",
    "    ax : `matplotlib.axes.Axes`, optional\n",
    "        If provided, the plot will be drawn on this axis (default is None).\n",
    "    nContours : int, optional\n",
    "        The number of contour lines to plot (default is 10).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    ax : `matplotlib.axes.Axes`\n",
    "        The axis on which the plot was drawn.\n",
    "    \"\"\"\n",
    "    data = exp.array\n",
    "    xlen, ylen = data.shape\n",
    "\n",
    "    vmin = np.percentile(data, 0.1)\n",
    "    vmax = np.percentile(data, 99.9)\n",
    "    lvls = np.linspace(vmin, vmax, nContours)\n",
    "\n",
    "    xx, yy = np.meshgrid(np.linspace(-xlen / 2, xlen / 2, xlen),\n",
    "                         np.linspace(-ylen / 2, ylen / 2, ylen))\n",
    "\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "\n",
    "    ax.contour(xx, yy, data, levels=lvls)\n",
    "    ax.tick_params(which=\"both\", direction=\"in\", top=True,\n",
    "                   right=True, labelsize=8)\n",
    "    ax.set_aspect(\"equal\")\n",
    "    ax.set_xlabel('x (pix)')\n",
    "    ax.set_ylabel('y (pix)')\n",
    "    ax.set_title(\"Contour plot\")\n",
    "    ax.set_xlim([-5, 5])\n",
    "    ax.set_ylim([-5, 5])\n",
    "\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af207cd0-3eb8-49d0-be1f-85f078ca9380",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T16:27:20.994415Z",
     "iopub.status.busy": "2023-07-05T16:27:20.993404Z",
     "iopub.status.idle": "2023-07-05T16:27:20.997680Z",
     "shell.execute_reply": "2023-07-05T16:27:20.996991Z",
     "shell.execute_reply.started": "2023-07-05T16:27:20.994378Z"
    },
    "tags": []
   },
   "source": [
    "## 2. Explore the PSF for a `calexp` image\n",
    "\n",
    "Recall that a `calexp` is the image from a single detector for a single processed visit image.\n",
    "\n",
    "### 2.1. Obtain a `calexp`\n",
    "\n",
    "Retrieve a particular single `calexp` image from the butler using a known `visit` and `detector`\n",
    "(this is the the same image as used in section 2 of DP0.2 Notebook Tutorial 3a)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39cf8d6d-4df2-45e6-9ceb-cf95887e527b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "datasetType = 'calexp'\n",
    "dataId = {'visit': 192350, 'detector': 175}\n",
    "calexp = butler.get(datasetType, dataId=dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bbe0819-eba3-4006-9e1b-4f8d823525e2",
   "metadata": {},
   "source": [
    "Display the `calexp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a774d0f6-cd45-40bb-86ff-9a19ab1f64aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5, 5))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(calexp.image)\n",
    "plt.title(f'Calexp dataID: {dataId}')\n",
    "plt.xlabel('x (pix)')\n",
    "plt.ylabel('y (pix)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "095bcbd5-9b84-46a5-a370-8ddba208f279",
   "metadata": {},
   "source": [
    "> **Figure 1:** *The plot above displays a 2D image of the `calexp` with coordinates ranging from 0 to 4000 pixels in both axes,\n",
    "and with a contrast bar ranging from -300 to 400 digital units.\n",
    "Point and extended sources are scattered around the image.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6cb02b8-4297-45cc-82a3-ad2b37570335",
   "metadata": {},
   "source": [
    "### 2.2. Extract PSF information\n",
    "\n",
    "Extract the PSF model from the exposure information, which is obtained from the `calexp` object itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75cbfc59-9fb1-4169-99aa-63fff065506d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "info_calexp = calexp.getInfo()\n",
    "psf_calexp = info_calexp.getPsf()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e1277a3-deb5-4574-bd32-74ef8930d138",
   "metadata": {},
   "source": [
    "**Option:** explore the methods available for `info_calexp` by uncommenting the following cell (removing the `#`),\n",
    "putting the cursor after the period, pressing the tab key, and scrolling through the list of options.\n",
    "Find that `getPsf` is just one option."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0811160a-99a2-43e5-97a7-83416d9cd82e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# info_calexp."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27af42b2-6a2e-4038-99f2-1a0932090682",
   "metadata": {},
   "source": [
    "**Option:** explore the methods available for `psf_calexp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dcbe9aa-ae63-458c-afbc-acf8745f17ab",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# psf_calexp."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0514d1f6-46e7-4c29-967c-c8ad3bbbb365",
   "metadata": {},
   "source": [
    "The methods available for `psf_calexp` will return PSF information for the PSF at a given point in the image.\n",
    "\n",
    "Pick the specific point in the image with coordinates `(x, y) = (2000, 3500)` and use\n",
    "it to define the `point_tuple` variable.\n",
    "\n",
    "Convert the `point_tuple` into a `Point2D` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d50cbe7-4729-4045-866c-ad9ca041a34e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "point_tuple = (2000, 3500)\n",
    "point_image = Point2D(point_tuple)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dafbc992-f17f-469e-8376-8b5c99fa52b8",
   "metadata": {},
   "source": [
    "Use the methods available for `psf_calexp` to retrieve and display a few\n",
    "properties of the PSF at the defined location in the image.\n",
    "\n",
    "Explore properties available for the the shape of the PSF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd50a13e-5674-4fa9-bd49-0d3953dc7e5e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "psf_shape = psf_calexp.computeShape(point_image)\n",
    "psf_shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e8b1fe0-587d-42f4-aed2-45051c4ef0c1",
   "metadata": {},
   "source": [
    "**Option:** explore the methods available for `psf_shape`.\n",
    "Commonly-used shape parameters can be accessed with `getArea`, `getDeterminantRadius`, `getIxx`, `getIyy`, and `getIxy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78ed53d0-1c4a-4f52-9113-e87b24d2b59e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# psf_shape."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a27b992-e80b-401c-aea4-586cda0e166b",
   "metadata": {},
   "source": [
    "Obtain the size of the PSF, as a $\\sigma$ and full-width half-max (FWHM)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db719c33-a0d4-4c38-87c2-f4f84e24fc2b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "psf_sigma = psf_shape.getDeterminantRadius()\n",
    "psf_fwhm = psf_sigma * SIGMA_TO_FWHM\n",
    "print(psf_sigma, psf_fwhm)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44920b55-eb25-4428-9f4a-e384007b9dfa",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-20T21:31:37.516462Z",
     "iopub.status.busy": "2023-09-20T21:31:37.515527Z",
     "iopub.status.idle": "2023-09-20T21:31:37.521373Z",
     "shell.execute_reply": "2023-09-20T21:31:37.520515Z",
     "shell.execute_reply.started": "2023-09-20T21:31:37.516424Z"
    },
    "tags": []
   },
   "source": [
    "Print the aperture flux within 1, 2, and 3$\\sigma$.\n",
    "Recall that the total flux of the PSF kernel has been normalized to an integrated flux of 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "462daa1a-c5ba-4f9a-a84f-ba94a34392a6",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "psf_apflux_1s = psf_calexp.computeApertureFlux(radius=1.0*psf_sigma, position=point_image)\n",
    "psf_apflux_2s = psf_calexp.computeApertureFlux(radius=2.0*psf_sigma, position=point_image)\n",
    "psf_apflux_3s = psf_calexp.computeApertureFlux(radius=3.0*psf_sigma, position=point_image)\n",
    "print(psf_apflux_1s, psf_apflux_2s, psf_apflux_3s)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88835313-967f-4df9-a12a-d361bb519a64",
   "metadata": {},
   "source": [
    "Print the peak flux of the PSF kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63518042-9438-4e4f-b998-e05ea38f3e5b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "psf_peak = psf_calexp.computePeak(position=point_image)\n",
    "psf_peak"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff71a24d-98d4-4918-b91d-7dffa9f2bc8f",
   "metadata": {},
   "source": [
    "Use the `getPsfProperties` function to return several commonly-used PSF properties.\n",
    "\n",
    "Pass the `psf_calexp` retrieved from the image, and the defined `point_image` coordinates,\n",
    "to the `getPsfProperties` helper function in order to obtain PSF properties at that location in the `calexp` image.\n",
    "\n",
    "The returned PSF properties are size (FWHM in pixels), \n",
    "aperture photometry flux (within a 1-sigma radius), \n",
    "peak flux value (recall that the PSF has been normalized to a sum of 1), \n",
    "and the dimensions of the PSF postage stamp (stamp display is demonstrated below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b0c9e9d-c6f3-4ffa-a720-76abdd78db24",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "props = getPsfProperties(psf_calexp, point_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc90e98b-b36f-40a2-a3f9-11242f6ed4ff",
   "metadata": {},
   "source": [
    "For comparison, return the PSF properties at another location in the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0b91ffa-04be-49d8-8955-ea150b724a99",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "point_tuple_2 = (100, 100)\n",
    "point_image_2 = Point2D(point_tuple_2)\n",
    "props_2 = getPsfProperties(psf_calexp, point_image_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aeab067-2ea8-4e71-bdc8-7be52780de62",
   "metadata": {},
   "source": [
    "### 2.3. Create and display PSF kernel images\n",
    "\n",
    "Use the extracted PSF information to create and display postage stamps of the PSF using\n",
    "two functions: `computeKernelImage`, which uses central pixels `(0, 0)`,\n",
    "and `computeImage`, which accepts user-specified central pixels.\n",
    "\n",
    "\n",
    "#### 2.3.1. Use `computeKernelImage`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "994ee3b5-7cf6-48ea-a4b1-e10067b93f4b",
   "metadata": {},
   "source": [
    "Use the `computeKernelImage` method to retrieve the PSF kernel at the location of the `point_image` defined above\n",
    "as a postage stamp with a center pixel of (0, 0).\n",
    "\n",
    "First compute the kernel image, and then convert it to `exposureF` format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b1b46f5-6fa8-4059-88fc-115113e84ed4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "psf_calexp_kernel = psf_calexp.computeKernelImage(point_image)\n",
    "first_psf_image_calexp = psf_calexp_kernel.convertF()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aa19867-59a3-4bc2-969d-a09a3fa1631d",
   "metadata": {},
   "source": [
    "Display the PSF postage stamp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ab8e49b-fb44-4d20-ae72-0e1c816cb49f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(3, 3))\n",
    "afw_display = afwDisplay.Display(frame=fig)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(first_psf_image_calexp)\n",
    "plt.xlabel('x (pix)')\n",
    "plt.ylabel('y (pix)')\n",
    "plt.gca().axis('on')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47a8571c-c84c-48c3-841c-a157391d4330",
   "metadata": {},
   "source": [
    "> **Figure 2:** *The plot above depicts the postage stamp (or cut-out) of the PSF model extracted from \n",
    "the `calexp` image using the `computeKernelImage` function.\n",
    "The PSF is centered on the origin, with central pixel coordintes (0, 0).\n",
    "The stamp is a square of 40 pixels per side that spans\n",
    "from negative 20 to positive 20 pixels on each side.\n",
    "The total flux in the PSF kernel has been normalized to sum to 1, and the color contrast bar\n",
    "spans a range from approximately negative 0.0002 to about 0.0003.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c991693c-b950-4b03-a396-649a7527983b",
   "metadata": {},
   "source": [
    "Note that with `computeKernelImage`, the postage stamp's coordinates are centered at the origin of the image.\n",
    "The coordinates of this origin point are (0, 0), resulting in negative coordinates for the lower left point.\n",
    "\n",
    "Print the coordinates for the lower left point of this postage stamp image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "161c6f7e-ead2-44d1-8a52-5fbb8fba61a4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(first_psf_image_calexp.getXY0())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4343915c-88e7-4980-b00b-3e77f67a7655",
   "metadata": {},
   "source": [
    "#### 2.3.2. Use `computeImage`\n",
    "\n",
    "Use the `computeImage` method to retrieve the PSF kernel at the location of the `point_image` defined above\n",
    "as a postage stamp with a center pixel defined by `point_image`.\n",
    "\n",
    "(I.e., instead of a PSF kernel postage stamp centered on `0, 0`, as returned by `computeKernelImage`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dad459b-4063-42c4-80f2-bc5f25e9a794",
   "metadata": {},
   "source": [
    "The `computeImage` function requires the astrometric solution for the `calexp`: in other words, \n",
    "the \"World Coordinate System\" (WCS) for the image is needed.\n",
    "The WCS maps pixel coordinates to sky coordinates and is obtained from the exposure information\n",
    "in a manner similar to how the PSF model was obtained.\n",
    "\n",
    "Obtain the WCS from the `calexp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "639aaf61-ed71-411b-ae76-42f755652b18",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wcs_calexp = info_calexp.getWcs()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a514f58-47e9-49af-893c-bbeffce05022",
   "metadata": {},
   "source": [
    "Compute the kernel image with `computeImage` and convert it to `exposureF` format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e309acb1-4213-4817-b64c-649557f35dbf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "second_psf_image_calexp = psf_calexp.computeImage(point_image).convertF()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9a3decc-a8e3-4f14-aabc-06f7a03bbe72",
   "metadata": {},
   "source": [
    "Display the PSF postage stamp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23265383-a8df-4afe-add5-fc5499dd6c48",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(3, 3))\n",
    "afw_display = afwDisplay.Display(frame=fig)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(second_psf_image_calexp)\n",
    "plt.xlabel('x (pix)')\n",
    "plt.ylabel('y (pix)')\n",
    "plt.gca().axis('on')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "246a4689-eb36-4e37-a86a-ed984fe02107",
   "metadata": {},
   "source": [
    "> **Figure 3:** *The plot above shows the postage stamp (or cut-out) of the PSF model \n",
    "from the `calexp` image obtained using `computeKernelImage`.\n",
    "The PSF is centered on the user-specified pixel coordinates (2000, 3500).\n",
    "The stamp is a square of 40 pixels per side that spans from 1980 to 2020 pixels on the `x` axis \n",
    "and from about 3480 to 3520 on the `y` axis.\n",
    "The total flux in the PSF kernel has been normalized to sum to 1, and the color contrast bar spans a range from approximately negative 0.0002 to about 0.0003.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64bc5b3a-b2e2-4642-9bad-29fbb16a6259",
   "metadata": {},
   "source": [
    "Print the coordinates for the lower left point of this postage stamp image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ca44e9b-1024-466e-82e5-35077156f11c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(second_psf_image_calexp.getXY0())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db08d1e2-115b-4d94-a052-928fd85cc853",
   "metadata": {},
   "source": [
    "#### 2.3.3. Save PSF kernel image to file\n",
    "\n",
    "Save copies of the PSF kernel images in the folder `/home/USER_NAME/DATA/.`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d588e76-4069-4e55-b502-da106c4a7cbf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_username = getpass.getuser()\n",
    "path_home = '/home/'+my_username+'/DATA/' \n",
    "print(path_home)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "851af851-0a34-49f7-96ff-4a0e9e7d5c4e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "first_psf_image_calexp.writeFits(path_home + 'first_psf_image.fits')\n",
    "second_psf_image_calexp.writeFits(path_home + 'second_psf_image.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4dbd5c3-34eb-4341-b092-343bc890d693",
   "metadata": {},
   "source": [
    "Use the file navigator in the left sidebar to locate the FITS files in the `DATA/` folder.\n",
    "Right-click on the filename and select Download. \n",
    "The downloaded files can then be opened with, e.g., [SAOImage ds9](https://sites.google.com/cfa.harvard.edu/saoimageds9) or other FITS file viewers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c731c48f-3126-4ca9-9bb5-b9549b0b2cb3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-05T18:22:30.461406Z",
     "iopub.status.busy": "2023-07-05T18:22:30.460243Z",
     "iopub.status.idle": "2023-07-05T18:22:30.469955Z",
     "shell.execute_reply": "2023-07-05T18:22:30.469108Z",
     "shell.execute_reply.started": "2023-07-05T18:22:30.461368Z"
    },
    "tags": []
   },
   "source": [
    "## 3. Explore the PSF for a `deepCoadd` image\n",
    "\n",
    "### 3.1. Context for coadded image PSFs\n",
    "\n",
    "When conducting multi-epoch surveys for static-sky science, the traditional method involves creating coadds.\n",
    "This process entails resampling images from different observations onto a common grid and combining them \n",
    "to generate a single deeper image (called a `deepCoadd` within the context of the LSST Science Pipelines).\n",
    "As part of the coaddition process, a coadded point spread function (PSF) model is established.\n",
    "Handling PSF coadding with care is crucial to ensure a well-defined PSF.\n",
    "\n",
    "One of the challenges in coadding PSFs arises from slight variations observed between PSFs in different visits.\n",
    "Even minor differences make it practically impossible to accurately model the effective PSF of the coadd using\n",
    "coadded star images.\n",
    "Small positional shifts, known as dithers, that are used to fill the gaps between charge-coupled devices \n",
    "(CCDs) introduce disruptions in the effective PSF of the coadd.\n",
    "As the number of dithers increases, the areas within the coadd with a continuous effective PSF become \n",
    "smaller, making it increasingly unlikely to find stars suitable for PSF modeling in each region.\n",
    "\n",
    "For example, for the Hyper-Suprime Camera survey, \n",
    "[Bosch et al 2018](https://ui.adsabs.harvard.edu/abs/2018PASJ...70S...5B/abstract) \n",
    "use an approach that involves resampling and combining existing PSF models from the input \n",
    "images using the same coordinate transformations and weights applied to the image data.\n",
    "The Hyper-Suprime Camera survey analysis pipeline is based on the LSST Science Pipelines code."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "222c76e3-37ee-4b18-9a3a-c6b1d9ea23f7",
   "metadata": {},
   "source": [
    "### 3.2. Obtain a `deepCoadd`\n",
    "\n",
    "Obtain the `deepCoadd` which covers the same point in the sky in the same filter as\n",
    "was used in Section 2, in order to better appreciate the differences between the \n",
    "`deepCoadd` and an individual `calexp` PSF model.\n",
    "\n",
    "Retrieve the right ascension and declination coordinates of the point used in the previous section \n",
    "(which was stored in the variable `point_tuple`), and print the coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d0a0a0f-d055-4d49-9035-56e33f76205d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "x, y = point_tuple\n",
    "ra, dec = wcs_calexp.pixelToSky(x, y)\n",
    "my_ra_deg = radToDeg(ra)\n",
    "my_dec_deg = radToDeg(dec)\n",
    "print(my_ra_deg, my_dec_deg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "164225f0-8589-4d18-a61b-448837171464",
   "metadata": {},
   "source": [
    "In order to retrieve a `deepCoadd` image, the butler will need the filter (`band`) to be specified in the `dataId`.\n",
    "\n",
    "Obtain the band of the `calexp` used above, and print it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c906ff57-a604-4984-89d9-72a509712727",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_band = info_calexp.getFilter().bandLabel\n",
    "print(my_band)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e76e390a-544c-4e22-974e-94b611d7cdc2",
   "metadata": {},
   "source": [
    "Follow the steps outlined in Section 3.2 of the DP0.2 tutorial notebook 01 (Introduction to DP0.2) \n",
    "to determine the `tract` and `patch` of the `deepCoadd`, based on the sky coordinates.\n",
    "Print the `tract` and `patch`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aba3928-01c2-4e95-8004-cc6f3191e47a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_spherePoint = SpherePoint(my_ra_deg*degrees,\n",
    "                             my_dec_deg*degrees)\n",
    "skymap = butler.get('skyMap')\n",
    "tract = skymap.findTract(my_spherePoint)\n",
    "patch = tract.findPatch(my_spherePoint)\n",
    "my_tract = tract.tract_id\n",
    "my_patch = patch.getSequentialIndex()\n",
    "print('my_tract: ', my_tract)\n",
    "print('my_patch: ', my_patch)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "729d070e-56ca-420b-855b-fe26ad7abcc9",
   "metadata": {},
   "source": [
    "Use the `tract`, `patch`, and `band` information to define the `dataId` for the desired `deepCoadd`.\n",
    "\n",
    "Retrieve the `deepCoadd` from the butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a3806f9-de0a-42a0-bcfd-c8b122ce326c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "datasetType = 'deepCoadd'\n",
    "dataId = {'tract': my_tract, 'patch': my_patch, 'band': my_band}\n",
    "coadd = butler.get(datasetType, dataId=dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d886e6f8-d98d-4726-bc09-8b533a8d9d8f",
   "metadata": {},
   "source": [
    "Display the `deepCoadd`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a88abef-9142-44f0-a22d-ba65098e6ffe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5, 5))\n",
    "display = afwDisplay.Display(frame=fig)\n",
    "display.scale('asinh', 'zscale')\n",
    "display.mtv(coadd.image)\n",
    "plt.title(f'deepCoadd. dataID: {dataId}')\n",
    "plt.xlabel('x (pix)')\n",
    "plt.ylabel('y (pix)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "614e7c42-3bb1-4cdd-8e11-99a543c60b51",
   "metadata": {},
   "source": [
    "> **Figure 4:** *The plot above displays the 2D image of the `deepCoadd` from the previous butler query.\n",
    "Coordinates range from 12000 to 16000 pixels in the horizontal axis and 0 to 4000 pixels in the vertical axis.\n",
    "The contrast bar ranges from negative 0.2 to 0.3 digital units. Point and extended sources are scattered around the image.\n",
    "Many more sources can be seen in this `deepCoadd` than in the `calexp` at the same sky location, shown in Figure 1.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3ab3ac8-e588-4fe5-9956-c580adda7e52",
   "metadata": {},
   "source": [
    "### 3.3. Extract PSF information\n",
    "\n",
    "Extract the PSF model from the `deepCoadd` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "485babcf-e6a5-47ce-8dd5-ae7ffb35453e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "info_coadd = coadd.getInfo()\n",
    "psf_coadd = info_coadd.getPsf()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ac87233-8dba-4f19-be6d-23fb1f51eeca",
   "metadata": {},
   "source": [
    "Obtain the pixel coordinates in the `deepCoadd` that correspond to the same\n",
    "sky coordinates (RA, Dec) used above \n",
    "(i.e., the sky coordinates corresponding to `calexp` pixel coordinates `(x, y) = (2000, 3500)`).\n",
    "\n",
    "Retrieve the WCS for the `deepCoadd`, then convert the RA and Dec stored as `my_spherePoint` into\n",
    "pixel coordinates for the `deepCoadd`, and store the pixel coordinates as `point_image`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46310f91-8936-481d-b35c-e4b99782067c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wcs_coadd = info_coadd.getWcs()\n",
    "point_image = wcs_coadd.skyToPixel(my_spherePoint)\n",
    "point_image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bc2c5a6-5527-437d-aa48-b874e3129ed6",
   "metadata": {},
   "source": [
    "Pass the `deepCoadd` PSF model and desired pixel coordinates to the `getPsfProperties` helper function \n",
    "to obtain the PSF properties of size, aperture flux, peak flux, and postage stamp dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "484d6b01-c1e0-471b-a905-e9c2af4252eb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "props_psf_coadd = getPsfProperties(psf_coadd, point_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24943e66-6a6f-475e-a28f-4b77a9105e56",
   "metadata": {},
   "source": [
    "> **Notice:** The PSF FWHM for the `deepCoadd` is 4.1 pixels -- larger than the PSF FWHM for the `calexp` which was 3.8 pixels.\n",
    "This indicates that we chose a `calexp` of better-than-average image quality, \n",
    "as the width of the `deepCoadd` FWHM\n",
    "is approximately indicative of the average image quality."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b44f6b9-480f-4e22-8e82-c69a298fd778",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-18T16:49:05.541062Z",
     "iopub.status.busy": "2023-09-18T16:49:05.540231Z",
     "iopub.status.idle": "2023-09-18T16:49:05.703038Z",
     "shell.execute_reply": "2023-09-18T16:49:05.702332Z",
     "shell.execute_reply.started": "2023-09-18T16:49:05.541038Z"
    },
    "tags": []
   },
   "source": [
    "### 3.4. Create and display PSF kernel images\n",
    "\n",
    "#### 3.4.1. Use `computeKernelImage`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7648e581-4bcd-4a76-847f-71477f3a9e9a",
   "metadata": {},
   "source": [
    "Use the `computeKernelImage` method to display the PSF kernel as a postage stamp with a center pixel of (0, 0).\n",
    "\n",
    "First compute the kernel image, and then convert it to `exposureF` format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9455bbb-1fa9-4ff7-860c-cc3223dbfefc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "psf_kernel_coadd = psf_coadd.computeKernelImage(point_image)\n",
    "first_psf_image_coadd = psf_kernel_coadd.convertF()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10997be2-331a-4ad2-ba1e-3bd3ecbcdda2",
   "metadata": {},
   "source": [
    "Display the PSF postage stamp for the `deepCoadd`.\n",
    "\n",
    "Note the higher signal-to-noise ratio for the `deepCoadd`, compared to the single `calexp` PSF image in Section 2.3.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab1d1c0c-d55f-47cb-acbe-109828b6a7f3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(3, 3))\n",
    "afw_display = afwDisplay.Display(frame=fig)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(first_psf_image_coadd)\n",
    "plt.xlabel('x (pix)')\n",
    "plt.ylabel('y (pix)')\n",
    "plt.gca().axis('on')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "819644ad-2cc7-40c2-adce-5080ee8cb264",
   "metadata": {},
   "source": [
    "> **Figure 5:** *The plot above displays the postage stamp of the PSF model from the `deepCoadd` image \n",
    "obtained using `computeKernelImage`. \n",
    "The PSF is centered at pixel coordinates (0, 0).\n",
    "The stamp is a square with sides of about 50 pixels, ranging from negative 25 to positive 25 pixels on both axes.\n",
    "The contrast bar ranges from -1e-5 to 7e-5.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d18805e7-b26c-45c4-92c5-b9d1443003cc",
   "metadata": {},
   "source": [
    "Print the coordinates of the lower-left corner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45d599ac-2c47-420b-8c9a-411fa34a4e06",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(first_psf_image_coadd.getXY0())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cf24788-f26e-4b1c-8752-514bca1a3687",
   "metadata": {},
   "source": [
    "#### 3.4.2. Use `computeImage`\n",
    "\n",
    "Use the `computeKernelImage` method to display the PSF as a postage stamp with a center pixel of \n",
    "(12937.5, 2325.0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a116e4c-8108-445e-a5c3-fc07efe06b09",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "second_psf_image_coadd = psf_coadd.computeImage(point_image).convertF()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1bcd692-3b79-4793-8ff2-f1b0daa71771",
   "metadata": {},
   "source": [
    "Display the PSF postage stamp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e17192f-b5d4-4817-8e11-ec0c2532d2fa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(3, 3))\n",
    "afw_display = afwDisplay.Display(frame=fig)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(second_psf_image_coadd)\n",
    "plt.xlabel('x (pix)')\n",
    "plt.ylabel('y (pix)')\n",
    "plt.gca().axis('on')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "635cef8f-5712-401a-b0d2-f0b31df6bd81",
   "metadata": {},
   "source": [
    "> **Figure 6:** *The plot above displays the postage stamp of the PSF model from the `deepCoadd` image \n",
    "obtained using `computeImage`. \n",
    "The PSF is centered at coordinates 12937.5, 2325.0 pixels.\n",
    "Te stamp is a square with sides of about 50 pixels, ranging from 12910 to 12970 pixels on the horizontal axis\n",
    "and 2300 to 2350 pixels on the vertical axis.\n",
    "The contrast bar ranges from pixels values of negative 2e-5 to positive 7e-5.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f5ea2b8-fd89-4ff2-ae9c-71b7090afd82",
   "metadata": {},
   "source": [
    "Print the coordinates of the lower-left corner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f4f53d3-80d2-4d0b-9d76-766b2124d4ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "second_psf_image_coadd.getXY0()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92385ff5-cb7a-4670-9c04-d3df99b947cf",
   "metadata": {},
   "source": [
    "## 4. Create PSF analysis plots\n",
    "\n",
    "Use the functions `plotRadialAverage`, `plotCurveOfGrowth`, and `plotContours`, \n",
    "which are based on image quality code currently employed for the \n",
    "[swift PSF analysis of the LATISS images](https://roundtable.lsst.codes/rubintv/summit/auxtel/im_current) \n",
    "obtained by Rubin's AuxTel at Cerro Pachón, Chile."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3fec3cc-7152-45b2-a05d-74753d265e86",
   "metadata": {},
   "source": [
    "### 4.1. Radial profile plots\n",
    "\n",
    "Use the `plotRadialAverage` function to plot an azimuthally-averaged one dimensional PSF profile \n",
    "and fit it with a Gaussian function, for both the `calexp` and `deepCoadd` PSF images. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce37032c-277c-47b2-ad9c-d49b8a24ddc5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax1 = plotRadialAverage(first_psf_image_calexp, ax=axes[0])\n",
    "ax1.set_title(\"Radial Average (calexp)\")\n",
    "ax2 = plotRadialAverage(first_psf_image_coadd, ax=axes[1])\n",
    "ax2.set_title(\"Radial Average (deepCoadd)\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16d0427-d191-4e88-b9b9-e9e20235e6ba",
   "metadata": {},
   "source": [
    "> **Figure 7:** *The plots above display the PSF kernel flux in ADU as a function of radius in pixels\n",
    "from the center of the PSF (blue crosses), and the best-fit Gaussian function (orange line), \n",
    "for the PSF from the `calexp` (left) and the `deepCoadd` (right).\n",
    "The legend contains the amplitude (peak flux), position (peak radial coordinate), \n",
    "and FWHM (in pixels), along with the [skewness](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skew.html) of the Gaussian fit.*\n",
    "\n",
    "As noted in Section 3.3, the FWHM of the `deepCoadd` PSF is larger than that for the selected `calexp`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9c30173-712a-4d04-b2d8-aa3b963ea679",
   "metadata": {},
   "source": [
    "### 4.2. 2D contour plots\n",
    "\n",
    "Use the `plotContours` function to create contour lines on the PSF image.\n",
    "These contour lines depict regions of equal brightness. \n",
    "As above, show the contour plots for the `calexp` and `deepCoadd` PSFs side-by-side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bfb5ed3-7bae-434a-8d32-1884c55c5ae2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax1 = plotContours(first_psf_image_calexp, ax=axes[0])\n",
    "ax1.set_title(\"Contours (calexp)\")\n",
    "ax1.set_xlim([-5, 5])\n",
    "ax1.set_ylim([-5, 5])\n",
    "ax2 = plotContours(first_psf_image_coadd, ax=axes[1])\n",
    "ax2.set_title(\"Contours (deepCoadd)\")\n",
    "ax2.set_xlim([-5, 5])\n",
    "ax2.set_ylim([-5, 5])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be322acf-b0ea-41b1-b31d-8e44102ce5e1",
   "metadata": {},
   "source": [
    "> **Figure 8:** *The plot above shows contour lines for the PSF kernel flux for the `calexp` (left) and \n",
    "the `deepCoadd` (right).\n",
    "Recall that the PSF kernel fluxes have been normalized to an integrated flux of 1.\n",
    "In both plots the inner and outer contour lines represent the 0.1 (yellow) and 99.9 (purple) percentile of the peak flux.*\n",
    "\n",
    "As for the radial profiles, notice that the PSF kernel is broader for the `deepCoadd`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f67dc1e-0f83-4484-b3e9-309c3cdcc9da",
   "metadata": {},
   "source": [
    "### 4.3. Curve of growth (encircled energy)\n",
    "\n",
    "Use the `plotCurveOfGrowth` function to calculate the cumulative amount of light enclosed by successively \n",
    "larger radii from the center of the PSF. \n",
    "This cumulative light measurement is normalized to the total light, resulting in a curve that reveals \n",
    "how light accumulates with increasing distance from the PSF center. \n",
    "As above, show the plots for the `calexp` and the `deepCoadd` PSFs side-by-side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6317c5ba-519f-4887-b22a-fbc9bc53a6cd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax1 = plotCurveOfGrowth(first_psf_image_calexp, ax=axes[0])\n",
    "ax1.set_title(\"calexp\")\n",
    "ax2 = plotCurveOfGrowth(first_psf_image_coadd, ax=axes[1])\n",
    "ax2.set_title(\"deepCoadd\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84c43b95-faa2-49b3-95cb-687c7736f55a",
   "metadata": {},
   "source": [
    "> **Figure 9:** *The plot above shows the curve of growth -- the fraction of the cumulative flux\n",
    "as a function of radial offset from the center of the PSF -- for the `calexp` (left) and the \n",
    "`deepCoadd` (right).*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "157656d2-fdaf-4a35-aecc-d6f3000db2c8",
   "metadata": {},
   "source": [
    "## 5. Exercises for the learner"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22b3ff6f-07cb-41f5-b884-f4ca9d623c24",
   "metadata": {},
   "source": [
    "**1. Co-plot `calexp` and `deepCoadd` PSF properties.**\n",
    "\n",
    "In Section 4, the PSF properties are plotted side-by-side.\n",
    "Use the code in the functions defined in Section 1.2 to co-plot the PSF properties for the\n",
    "`calexp` and the `deepCoadd` in the same panel (same axis). \n",
    "\n",
    "HINT: This is relatively easy for the curve of growth in Section 4.3.\n",
    "The code below works.\n",
    "So does a similar approach for the radial profiles.\n",
    "But co-plotting the contours is more challenging.\n",
    "\n",
    "```\n",
    "fig, axes = plt.subplots(1, 1, figsize=(4, 4))\n",
    "ax1 = plotCurveOfGrowth(first_psf_image_calexp, ax=axes)\n",
    "ax1.set_title(\"calexp\")\n",
    "ax1 = plotCurveOfGrowth(first_psf_image_coadd, ax=axes)\n",
    "ax1.set_title(\"deepCoadd\")\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "```\n",
    "\n",
    "For an even greater challenge, write new functions that *extract and return* \n",
    "the data being plotted (e.g., a new function that *returns* the `distances` and `values` \n",
    "arrays in `plotRadialAverage`, instead of plotting them and returning the axis),\n",
    "then create your own plots with the returned arrays."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59ee62ca-3725-4f10-acf9-124dcbfa78b8",
   "metadata": {},
   "source": [
    "**2. Fit the PSF radial profiles with a non-Gaussian function.**\n",
    "\n",
    "The function `plotRadialAverage` uses a Gaussian function to fit the averaged PSF profile.\n",
    "Write a new function that uses a different (non-Gaussian) functional form for fitting the one-dimensional PSF radial profiles, \n",
    "such as the Moffat profile (as described in Section 3.1 of \n",
    "[Jarvis et al. 2021](https://ui.adsabs.harvard.edu/abs/2021MNRAS.501.1282J/abstract))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7d2f5d5-aadb-4473-a60d-2f04f3ea6a39",
   "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.9"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}