Visualization/LsstCameraGeometry.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LSST Camera Geometry\n",
    "\n",
    "<br>Owner(s): **Alex Drlica-Wagner** ([@kadrlica](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@kadrlica))\n",
    "<br>Last Verified to Run: **2021-09-11**\n",
    "<br>Verified Stack Release: **w_2021_33**\n",
    "\n",
    "This notebook demonstrates how to interact with the Stack representation of the LSST camera (`lsstCam`). We use it to get the geometry of the various nested components of the camere focal plane -- i.e., amps, detectors, and rafts. We then produce a labeled figure of the LSST focal plane geometry.\n",
    "\n",
    "### Learning Objectives:\n",
    "\n",
    "After working through this tutorial you should be able to:\n",
    "\n",
    "1. Access the geometry of specific LSST amps, detectors, and rafts.\n",
    "2. Plot the geometry of these camera components.\n",
    "3. Create a labeled plot of the LSST focal plane geometry.\n",
    "\n",
    "### Logistics\n",
    "This notebook is intended to be run at `lsst-lsp-stable.ncsa.illinois.edu` or `data.lsst.cloud` from a local git clone of the [StackClub](https://github.com/LSSTScienceCollaborations/StackClub) repo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "You can find the Stack version by using `eups list -s` on the terminal command line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Site, host, and stack version\n",
    "! echo $EXTERNAL_INSTANCE_URL\n",
    "! echo $HOSTNAME\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from matplotlib.patches import Polygon\n",
    "from matplotlib.collections import PatchCollection\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from lsst.daf.butler import Butler \n",
    "from lsst.afw.cameraGeom import utils as cgUtils\n",
    "from lsst.obs.lsst.cameraTransforms import LsstCameraTransforms\n",
    "from lsst.obs.lsst.cameraTransforms import ampPixelToCcdPixel\n",
    "\n",
    "import lsst.afw.cameraGeom as cameraGeom\n",
    "import lsst.afw.geom as afwGeom\n",
    "\n",
    "from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Existing Tools\n",
    "\n",
    "The stack provides a pre-wrapped visualization of the camera geometry. However, this figure has a few limitations: the detector labels are very small, rafts are only indicated in the sensor names, and amplifier geometry is not shown. A more visually apealing version of the LSST camera geometry can be found [here](https://confluence.lsstcorp.org/display/LSWUG/Representation+of+a+Camera?preview=/4129064/10190878/LSST_FocalPlane.png#/); however, it has a different orientation than is used by the DM Stack."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "URL = os.getenv('EXTERNAL_INSTANCE_URL')\n",
    "if URL.endswith('data.lsst.cloud'): # IDF\n",
    "    repo = \"s3://butler-us-central1-dp01\"\n",
    "elif URL.endswith('ncsa.illinois.edu'): # NCSA\n",
    "    repo = \"/repo/dc2\"\n",
    "else:\n",
    "    raise Exception(f\"Unrecognized URL: {URL}\")\n",
    "\n",
    "dataset='DC2'\n",
    "collection='2.2i/runs/DP0.1'\n",
    "\n",
    "butler = Butler(repo,collections=collection)\n",
    "camera = butler.get('camera',instrument='LSSTCam-imSim')\n",
    "cgUtils.plotFocalPlane(camera)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The source code for `cgUtils.plotFocalPlane` found [here](https://github.com/lsst/afw/blob/01e196b6519ef91d51c61435065e16477972b897/python/lsst/afw/cameraGeom/utils.py#L87) gives us a very good place to start our investigation of the LSST camera geometry."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot CCDs on Focal Plane\n",
    "\n",
    "To get started, it's helpful to understand how to reproduce the figure above. The basic idea is to build an instance of the LSST camera object, and then loop through the detectors that it contains. For each detector, we plot the detector location and label the detector. We make a slight change to the code above by labeling with the CCD ID (a unique integer for each detector) rather than the `name` of the detector (which is built from the raft/sensor information)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig= plt.figure(figsize=(10,10))\n",
    "ax = plt.gca()\n",
    "\n",
    "xvals,yvals = [],[]\n",
    "colors,patches = [],[]\n",
    "for det in camera:\n",
    "    corners = [(c.getX(), c.getY()) for c in det.getCorners(FOCAL_PLANE)]\n",
    "    for corner in corners:\n",
    "        xvals.append(corner[0])\n",
    "        yvals.append(corner[1])\n",
    "    colors.append('b')\n",
    "    patches.append(Polygon(corners, True))\n",
    "    center = det.getOrientation().getFpPosition()\n",
    "    text = det.getId()\n",
    "    rot = 90 if text in (195,196,199,200) else 0\n",
    "    ax.text(center.getX(), center.getY(), text,\n",
    "            ha='center', va='center', size=12, rotation=rot)\n",
    "\n",
    "patchCollection = PatchCollection(patches, alpha=0.6, facecolor=colors)\n",
    "ax.add_collection(patchCollection)\n",
    "\n",
    "ax.set_xlim(min(xvals) - abs(0.1*min(xvals)),\n",
    "            max(xvals) + abs(0.1*max(xvals)))\n",
    "ax.set_ylim(min(yvals) - abs(0.1*min(yvals)),\n",
    "            max(yvals) + abs(0.1*max(yvals)))\n",
    "ax.set_xlabel('Focal Plane X (mm)')\n",
    "ax.set_ylabel('Focal Plane Y (mm)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Amps on a CCD\n",
    "\n",
    "We'll start by grabbing the amps from a specific detector (the central detector R22,S11). We can then plot the extents of the amps in pixel coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the camera transformation\n",
    "transform = LsstCameraTransforms(camera)\n",
    "# Get the central detector\n",
    "det = transform.getDetector('R22_S11')\n",
    "\n",
    "# Get the amps for this detector\n",
    "amps = det.getAmplifiers()\n",
    "\n",
    "fig,ax = plt.subplots(1,figsize=(10,10))\n",
    "patches,colors = [],[]\n",
    "for amp in amps:\n",
    "    corners = [(c.getX(), c.getY()) for c in amp.getBBox().getCorners()]\n",
    "    patches.append(Polygon(corners, True))\n",
    "    colors.append('b')\n",
    "    center = amp.getBBox().getCenter()\n",
    "    text = amp.getName()\n",
    "    ax.text(center.getX(), center.getY(), text, color='k',\n",
    "            ha='center', va='center', size=14)\n",
    "\n",
    "# Add the patch collection\n",
    "patchCollection = PatchCollection(patches, alpha=0.6, facecolor=colors,edgecolor='k')\n",
    "ax.add_collection(patchCollection)\n",
    "\n",
    "# Set some labels and extent\n",
    "ax.set_xlim(-200,4250)\n",
    "ax.set_ylim(-200,4250)\n",
    "ax.set_xlabel('CCD X (pix)')\n",
    "ax.set_ylabel('CCD Y (pix)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so this is all well and good, but what if we want to plot the physical positions of the amps in focal plane coordinates? We should be able to do this with a transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotFancyDetector(camera, detectorName='R22_S11', figsize=(10.,10.)):\n",
    "    # Create the camera transformation\n",
    "    transform = LsstCameraTransforms(camera)\n",
    "    # Get the central detector\n",
    "    det = transform.getDetector(detectorName)\n",
    "\n",
    "    # Get the amps for this detector\n",
    "    amps = det.getAmplifiers()\n",
    "\n",
    "    fig,ax = plt.subplots(1,figsize=figsize)\n",
    "    patches,colors = [],[]\n",
    "    xvals,yvals = [],[]\n",
    "    for amp in amps:\n",
    "        points = [transform.ccdPixelToFocalMm(c.getX(), c.getY(), det.getName()) for c in amp.getBBox().getCorners()]\n",
    "        corners = [(p.getX(),p.getY()) for p in points]\n",
    "        for corner in corners:\n",
    "            xvals.append(corner[0])\n",
    "            yvals.append(corner[1])\n",
    "        patches.append(Polygon(corners, True))\n",
    "        colors.append('skyblue')\n",
    "        # Center in pixels\n",
    "        center_pix = amp.getBBox().getCenter()\n",
    "        # center in mm\n",
    "        center = transform.ccdPixelToFocalMm(center_pix.getX(),center_pix.getY(),det.getName())\n",
    "        text = amp.getName()\n",
    "        ax.text(center.getX(), center.getY(), text, color='k',\n",
    "                ha='center', va='center', size=14)\n",
    "\n",
    "    # Add the patch collection\n",
    "    patchCollection = PatchCollection(patches, alpha=0.6, facecolor=colors,edgecolor='k')\n",
    "    ax.add_collection(patchCollection)\n",
    "\n",
    "    ax.set_xlim(min(xvals) - abs(0.02*min(xvals)),\n",
    "                max(xvals) + abs(0.02*max(xvals)))\n",
    "    ax.set_ylim(min(yvals) - abs(0.02*min(yvals)),\n",
    "                max(yvals) + abs(0.02*max(yvals)))\n",
    "    ax.set_xlabel('Focal Plane X (mm)')\n",
    "    ax.set_ylabel('Focal Plane Y (mm)')\n",
    "    ax.set_title(det.getName())\n",
    "\n",
    "# Set some labels and extent\n",
    "plotFancyDetector(camera,'R11_S11')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot CCDs on a Raft\n",
    "\n",
    "It looks like the LSST Camera object doesn't have the concept of a raft (just a list of detectors). If we want to assemble a raft, we can do so directly from the CCDs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotFancyRaft(camera, raftName='R22', figsize=(10.,10.)):\n",
    "    colorMap = {0: 'skyblue', 1: 'y', 2: 'g', 3: 'r'}\n",
    "    transform = LsstCameraTransforms(camera)\n",
    "\n",
    "    plt.figure(figsize=figsize)\n",
    "    ax = plt.gca() \n",
    "\n",
    "    patches, colors = [],[]\n",
    "    xvals, yvals = [],[]\n",
    "    for det in camera:\n",
    "        if not det.getName().startswith(raftName): continue\n",
    "        corners = [(c.getX(), c.getY()) for c in det.getCorners(FOCAL_PLANE)]\n",
    "        for corner in corners:\n",
    "            xvals.append(corner[0])\n",
    "            yvals.append(corner[1])\n",
    "        colors.append(colorMap[int(det.getType())])\n",
    "        patches.append(Polygon(corners, True))\n",
    "        center = det.getOrientation().getFpPosition()\n",
    "        name = det.getName()\n",
    "\n",
    "        # Label central raft\n",
    "        text = '(%s,%s)'%tuple(name.split('_')[1].strip('S'))\n",
    "        ax.text(center.getX(), center.getY(), text, color='0.3',\n",
    "                ha='center', va='center',size=18)      \n",
    "             \n",
    "    patchCollection = PatchCollection(patches, alpha=0.6, facecolor=colors)\n",
    "    ax.add_collection(patchCollection)\n",
    "    \n",
    "    ax.set_xlim(min(xvals) - abs(0.1*min(xvals)),\n",
    "                max(xvals) + abs(0.1*max(xvals)))\n",
    "    ax.set_ylim(min(yvals) - abs(0.1*min(yvals)),\n",
    "                max(yvals) + abs(0.1*max(yvals)))\n",
    "    ax.set_xlabel('Focal Plane X (mm)')\n",
    "    ax.set_ylabel('Focal Plane Y (mm)')\n",
    "    ax.set_title(raftName)\n",
    "    \n",
    "plotFancyRaft(camera)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Focal Plane\n",
    "\n",
    "Now we put together what we've learned to try to replicate the image [here](https://confluence.lsstcorp.org/display/LSWUG/Representation+of+a+Camera?preview=/4129064/10190878/LSST_FocalPlane.png#/). Note that the corner rafts (containing the focus and guiding CCDs) are not yet included in the `lsstCamera` model. Also, since the concept of a raft doesn't really exist, we don't plot outlines around the rafts (we could hack this if we wanted)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotFancyFocalPlane(camera, figsize=(10., 10.), showFig=True, savePath=None):\n",
    "    \"\"\"Make a plot of the focal plane along with a set points that sample\n",
    "    the field of view.\n",
    " \n",
    "    Parameters\n",
    "    ----------\n",
    "    camera : `lsst.afw.cameraGeom.Camera`\n",
    "        A camera object\n",
    "    figsize : `tuple` containing two `float`\n",
    "        Matplotlib style tuple indicating the size of the figure in inches\n",
    "    showFig : `bool`\n",
    "        Display the figure on the screen?\n",
    "    savePath : `str` or `None`\n",
    "        If not `None`, save a copy of the figure to this name.\n",
    "    \"\"\"\n",
    "    try:\n",
    "        from matplotlib.patches import Polygon\n",
    "        from matplotlib.collections import PatchCollection\n",
    "        import matplotlib.pyplot as plt\n",
    "    except ImportError:\n",
    "        raise ImportError(\n",
    "            \"Can't run plotFocalPlane: matplotlib has not been set up\")\n",
    "\n",
    "    colorMap = {0: 'skyblue', 1: 'y', 2: 'g', 3: 'r'}\n",
    "    transform = LsstCameraTransforms(camera)\n",
    "\n",
    "    plt.figure(figsize=figsize)\n",
    "    ax = plt.gca()\n",
    "\n",
    "    patches, colors = [],[]\n",
    "    xvals, yvals = [],[]\n",
    "    for det in camera:\n",
    "        corners = [(c.getX(), c.getY()) for c in det.getCorners(FOCAL_PLANE)]\n",
    "        for corner in corners:\n",
    "            xvals.append(corner[0])\n",
    "            yvals.append(corner[1])\n",
    "        colors.append(colorMap[int(det.getType())])\n",
    "        patches.append(Polygon(corners, True))\n",
    "        center = det.getOrientation().getFpPosition()\n",
    "        name = det.getName()\n",
    "\n",
    "        # Label central raft\n",
    "        if name.startswith('R22'):\n",
    "            if not name.endswith('S11'):\n",
    "                # Label CCDs for central raft\n",
    "                text = '(%s,%s)'%tuple(name.split('_')[-1].strip('S'))\n",
    "                ax.text(center.getX(), center.getY(), text,\n",
    "                        ha='center', va='center', size=10)\n",
    "            else:\n",
    "                # Draw the amps for the central CCD\n",
    "                amp_patches = []\n",
    "                for amp in det.getAmplifiers():\n",
    "                    points = [transform.ccdPixelToFocalMm(c.getX(), c.getY(), det.getName()) for c in amp.getBBox().getCorners()]\n",
    "                    corners = [(p.getX(),p.getY()) for p in points]\n",
    "                    amp_patches.append(Polygon(corners, True))\n",
    "                # Add the amp patch collection\n",
    "                patchCollection = PatchCollection(amp_patches, alpha=0.6, facecolor='none',edgecolor='k')\n",
    "                ax.add_collection(patchCollection)\n",
    "        elif name.endswith('S11'):\n",
    "            text = '(%s,%s)'%tuple(name.split('_')[0].strip('R'))\n",
    "            ax.text(center.getX(), center.getY(), text, color='0.3',\n",
    "                    ha='center', va='center',size=22)\n",
    "    for raft in ('R00', 'R40', 'R44', 'R04'):\n",
    "        # These rafts don't have an S11 sensor, so need to figure out the raft center from the other rafts around them\n",
    "        _, y, x = list(raft)\n",
    "        column = camera[f'R1{x}_S11'].getOrientation().getFpPosition()  # Just needs to be in the column, could have used R2{x}\n",
    "        row = camera[f'R{y}1_S11'].getOrientation().getFpPosition()  # Same for rows\n",
    "        text = f'({y},{x})'\n",
    "        ax.text(column.getX(), row.getY(), text, color='0.3',\n",
    "                ha='center', va='center',size=22)\n",
    "        \n",
    "             \n",
    "    patchCollection = PatchCollection(patches, alpha=0.6, facecolor=colors)\n",
    "    ax.add_collection(patchCollection)\n",
    "    \n",
    "    ax.set_xlim(min(xvals) - abs(0.1*min(xvals)),\n",
    "                max(xvals) + abs(0.1*max(xvals)))\n",
    "    ax.set_ylim(min(yvals) - abs(0.1*min(yvals)),\n",
    "                max(yvals) + abs(0.1*max(yvals)))\n",
    "    \n",
    "    ax.set_xlabel('Focal Plane X (mm)',fontsize=18)\n",
    "    ax.set_ylabel('Focal Plane Y (mm)',fontsize=18)\n",
    "    if savePath is not None:\n",
    "        plt.savefig(savePath)\n",
    "    if showFig:\n",
    "        plt.show()\n",
    "\n",
    "# Plot the focal plane        \n",
    "plotFancyFocalPlane(camera)"
   ]
  }
 ],
 "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
}