obsTAP-raw-advanced.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bd1cf3c8-c8d1-4867-a73b-ce50d38163f8",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=170 style=\"padding: 10px\"> \n",
    "<b>Little Demo: Raw images from ObsTAP</b> <br>\n",
    "Contact author(s): Melissa Graham <br>\n",
    "Last verified to run: 2024-07-20 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_16 <br>\n",
    "Container Size: medium"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5e67efc-f77e-4e4e-abf0-6f58bb3aeb43",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "The simulted images for Rubin's Data Preview 0 (DP0) are available via both the\n",
    "[data butler](https://pipelines.lsst.io/modules/lsst.daf.butler/index.html)\n",
    "and the TAP (Table Access Protocol) service.\n",
    "\n",
    "Use of the data butler to query and retrieve images\n",
    "is covered in other tutorials, and is generally the recommended\n",
    "way to access DP0 images while working within the Notebook Aspect of the \n",
    "Rubin Science Platform, especially if any image reprocessing is being done.\n",
    "\n",
    "The Portal Aspect of the Rubin Science Platform offers a graphical user interface to the TAP service,\n",
    "and images can be displayed interactively in the Portal.\n",
    "No code experience is necessary to use the Portal.\n",
    "\n",
    "However, raw image access via the TAP service is also possible.\n",
    "This tutorial demonstrates how to retrieve and display raw images in a Notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "678c2593-c625-4129-8168-8ec111e81910",
   "metadata": {},
   "source": [
    "### 1.1. Import packages\n",
    "\n",
    "Import packages a few standard python packages,\n",
    "[lsst packages](https://pipelines.lsst.io/) for image display and data access,\n",
    "[astropy packages](https://www.astropy.org/) for working with images,\n",
    "and a module from the [pvyo.dal package](https://pyvo.readthedocs.io/en/latest/dal/index.html)\n",
    "for data access."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e29ed16-5f40-4118-8c60-74c1afae8e7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from urllib.request import urlretrieve\n",
    "\n",
    "import lsst.afw.display as afwDisplay\n",
    "import lsst.geom as geom\n",
    "from lsst.rsp import get_tap_service\n",
    "from lsst.afw.image import ExposureF\n",
    "\n",
    "import astropy.visualization as vis\n",
    "from astropy.io import fits\n",
    "from astropy.wcs import WCS\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy import units as u\n",
    "\n",
    "import warnings\n",
    "from astropy.wcs import FITSFixedWarning\n",
    "\n",
    "from pyvo.dal.adhoc import DatalinkResults"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cfee068-f967-4e45-83a2-2e085297ece8",
   "metadata": {},
   "source": [
    "### 1.2. Define functions and paramters\n",
    "\n",
    "Instantiate the RSP's TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d818a521-b241-45a8-8ca6-0addf8b38051",
   "metadata": {},
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c58604ec-879e-4b26-ba25-97cbeceaf71c",
   "metadata": {},
   "source": [
    "Define `auth_session` in order to authorize the TAP service to\n",
    "access LSST images via Datalink in Section 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad507d6f-537c-472c-864e-17fe2ea9acb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "auth_session = service._session"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ea0ac5d-63e2-4715-82a5-913163b0c693",
   "metadata": {},
   "source": [
    "Define a point near the center of the Data Preview 0 area with\n",
    "coordintes Right Ascension (RA) = 62 degrees\n",
    "and Declination (Dec) = -37 degrees.\n",
    "Use several different formats for the different use cases\n",
    "throughout this tutorial\n",
    "\n",
    " * `target_ra`, `target_dec` : floats, to be used in calculations\n",
    " * `target_str_ra`, `target_str_dec` : strings, to be inserted into query statements\n",
    " * `targetPoint` : a `SpherePoint` object from the LSST Science Pipelines' `geom` package\n",
    " * `targetCoord` : an astropy `SkyCoord` class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62d01de1-7103-4a24-bf0e-1159d1af4ae0",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_ra = 62.0\n",
    "target_dec = -37.0\n",
    "target_str_ra = '62.0'\n",
    "target_str_dec = '-37.0'\n",
    "targetPoint = geom.SpherePoint(target_ra*geom.degrees,\n",
    "                               target_dec*geom.degrees)\n",
    "targetCoord = SkyCoord(ra=target_ra*u.degree,\n",
    "                       dec=target_dec*u.degree,\n",
    "                       frame='icrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7aa9221a-586d-4218-8ba9-ec01e8585f99",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-07-04T02:06:34.143333Z",
     "iopub.status.busy": "2024-07-04T02:06:34.142887Z",
     "iopub.status.idle": "2024-07-04T02:06:34.279290Z",
     "shell.execute_reply": "2024-07-04T02:06:34.278550Z",
     "shell.execute_reply.started": "2024-07-04T02:06:34.143312Z"
    }
   },
   "source": [
    "Define a function `detsec_to_x1x2y1y2` to be used in Section 4.3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66d86776-20e6-4bed-ab11-ce91daf5c65a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def detsec_to_x1x2y1y2(detsec: str):\n",
    "    \"\"\"Convert the DETSEC from a FITS header into the\n",
    "    start and end pixels in x and y.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    detsec: `str`\n",
    "        String formatted as detector section, e.g.,:\n",
    "        \"[509:1,1:2000]\" or \"[1018:510,1:2000]\"\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    x1: `int`\n",
    "        The start x-pixel.\n",
    "    x2: `int`\n",
    "        The end x-pixel.\n",
    "    y1: `int`\n",
    "        The start y-pixel.\n",
    "    y2: `int`\n",
    "        The end y-pixel.\n",
    "    \"\"\"\n",
    "    temp1 = detsec.strip('[').strip(']')\n",
    "    temp2 = temp1.split(',')\n",
    "    temp3 = temp2[0].split(':')\n",
    "    temp4 = temp2[1].split(':')\n",
    "    x1 = int(temp3[0])\n",
    "    x2 = int(temp3[1])\n",
    "    y1 = int(temp4[0])\n",
    "    y2 = int(temp4[1])\n",
    "    return x1, x2, y1, y2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dde71a05-9803-419e-9693-08e423029bc2",
   "metadata": {},
   "source": [
    "## 2. Raw Images\n",
    "\n",
    "Follow a similar series of steps as in Section 4.1 to retrieve\n",
    "the pixel data for raw images.\n",
    "\n",
    "### 2.1. Query and retrieve\n",
    "\n",
    "One difference between here and Section 4.1 is to constrain the\n",
    "`obs_id` parameter to be the same as the `calexp`, in order\n",
    "to obtain the `raw` version of the `calexp` investigated above.\n",
    "\n",
    "Define the query, submit the job, and retrieve the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdb94829-9613-4f26-9c45-00e3f855e050",
   "metadata": {},
   "outputs": [],
   "source": [
    "use_obs_id = '214437-R20_S20'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bae656bf-7c76-4dd7-9995-956841aa3fd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "query = \"SELECT * FROM ivoa.ObsCore \"\\\n",
    "        \"WHERE CONTAINS(POINT('ICRS', \" + \\\n",
    "        target_str_ra + \", \" + target_str_dec + \"), s_region) = 1 \"\\\n",
    "        \"AND lsst_band = 'i' \"\\\n",
    "        \"AND calib_level = 1 \"\\\n",
    "        \"AND obs_id = '\" + use_obs_id + \"'\"\n",
    "print(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "351c228e-1a6d-4adf-b4f4-a165a04c8787",
   "metadata": {},
   "outputs": [],
   "source": [
    "job = service.submit_job(query)\n",
    "job.run()\n",
    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
    "print('Job phase is', job.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75ccb48-c07a-4da0-a46b-775c1ecd0db6",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = job.fetch_result()\n",
    "print(len(results))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f05691eb-2962-4241-8aa4-1761f4e981ea",
   "metadata": {},
   "source": [
    "Take the first (only!) row and print the same metadata as was printed for the `calexp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54c20dc3-335f-48cb-81a0-0ad1598367c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "results_0 = results[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6e0098a-5983-4f19-ade8-3ee4d002ae60",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_0['dataproduct_subtype'])\n",
    "print(results_0['lsst_ccdvisitid'])\n",
    "print(results_0['lsst_visit'], results_0['lsst_detector'])\n",
    "print('%7.4f %8.4f' % (results_0['s_ra'], results_0['s_dec']))\n",
    "print(results_0['lsst_band'])\n",
    "print(results_0['t_min'])\n",
    "print(results_0['s_region'])\n",
    "print(results_0['obs_id'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "351b9bb1-c6e4-4dc8-9311-0b63ba89090b",
   "metadata": {},
   "source": [
    "Above, notice that for the `raw` image, the `ccdvisitid` and `visit` are not defined\n",
    "(this is a temporary, DP0-era issue).\n",
    "Also notice that the `raw` image search results are sorted differently from the `calexp`, and \n",
    "this is a different image (different date, different observation id).\n",
    "\n",
    "Use `Datalink` to obtain the `image_url`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78a953f5-ef80-464e-b3d2-fed4bbc11059",
   "metadata": {},
   "outputs": [],
   "source": [
    "dl_results = DatalinkResults.from_result_url(results_0['access_url'],\n",
    "                                             session=auth_session)\n",
    "image_url = dl_results['access_url'][0]\n",
    "print(image_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22d1fb03-a9e7-4c91-ba63-9809db4688bb",
   "metadata": {},
   "source": [
    "Print header componenents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c2e4eb0-6173-440d-9d99-826a557382fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "hdulist = fits.open(image_url)\n",
    "for hdu in hdulist:\n",
    "    print(hdu.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "307172fe-501d-410d-b5ee-0b77460eeaa5",
   "metadata": {},
   "source": [
    "The above shows that `raw` images (which are one of the LSST Science Camera's 189 detectors) are divided into 16 segments (one for each of the detector's amplifiers).\n",
    "\n",
    "**Optional:** print the primary header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f49d329-2195-4993-8603-88b54e1fc64a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# hdulist[0].header"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96e246b3-3bd8-4aff-ba02-719f88e9aea5",
   "metadata": {},
   "source": [
    "### 2.2. Display one segment\n",
    "\n",
    "Display one segment with `matplotlib`.\n",
    "\n",
    "The following cell will return a few error messages about `astropy` changing\n",
    "some of the header keywords. It is ok to ignore these for now, and they\n",
    "will be fixed as the LSST-output FITS files continue to evolve.\n",
    "\n",
    "To suppress the output of these warnings, use:\n",
    "```\n",
    "import warnings\n",
    "from astropy.wcs import FITSFixedWarning\n",
    "warnings.filterwarnings('ignore', category=FITSFixedWarning)\n",
    "```\n",
    "\n",
    "Get image header and data for the first segment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "090b75f3-79c7-4c97-a131-9b01403cb6ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_hdr = hdulist[1].header\n",
    "img_wcs = WCS(img_hdr)\n",
    "img_seg = hdulist[1].name\n",
    "img_data = hdulist[1].data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "047caab3-5e9f-4517-9841-505609df12a9",
   "metadata": {},
   "source": [
    "**Option** to view the segment header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e6b6eda-1370-4777-af63-a69d49f6ade1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# img_hdr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50e91929-b7e4-42b5-acd8-e47d98a1cae4",
   "metadata": {},
   "outputs": [],
   "source": [
    "zscale = vis.ZScaleInterval()\n",
    "zlimits = zscale.get_limits(img_data)\n",
    "print(zlimits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33d732f3-6800-4e5b-87e9-5736a2877eba",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(3, 5))\n",
    "plt.subplot(projection=img_wcs)\n",
    "plt.imshow(img_data, cmap='grey',\n",
    "           vmin=zlimits[0], vmax=zlimits[1])\n",
    "plt.axis('on')\n",
    "plt.grid(color='white', ls='solid')\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.title(img_seg)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36e4dc88-a1b7-4bc8-83a3-7acb5c91f5e9",
   "metadata": {},
   "source": [
    "Figure 4: One of the segments of the `raw` image, plotted in the same manner as the `calexp` in Section 4.1.5.\n",
    "\n",
    "What's with the black pixels? They represent the image's bias (a pedestal level of counts added by the amplifier).\n",
    "In the header, `NAXIS` and `NAXIS2` represent the shape of the image array, and\n",
    "`DATASEC` provides the dimensions of the pixels with flux.\n",
    "Note that the `numpy.shape` is (rows, columns), whereas for the FITS file the (axis 1, axis 2) is (columns, rows)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fc1dde1-c36d-482e-9d32-6706d2d050c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(img_hdr['NAXIS1'], img_hdr['NAXIS2'])\n",
    "print(np.shape(img_data))\n",
    "print(img_hdr['DATASEC'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e713d808-1e1b-4423-b1cf-1999d9dbdc35",
   "metadata": {},
   "source": [
    "The overscan (bias) pixels are in columns 1, 2, 3, and 513 through 544, and in rows 2001 through 2048.\n",
    "\n",
    "Compare the mean flux of columns 2000 (flux) and 2001 (overscan).\n",
    "Index the image array with N-1 because python indexes from 0 whereas FITS indexes from 1.\n",
    "The bias level is ~1000 counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "556e8642-9dc9-48e7-9864-ae9d6203705b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('mean counts in column 2000: %7.2f,  column 2001: %7.2f' %\n",
    "      (np.mean(img_data[1999]), np.mean(img_data[2000])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5832a34-c201-47a5-b337-d949ef650122",
   "metadata": {},
   "source": [
    "### 2.3. Display mosaic\n",
    "\n",
    "Now that the warnings have been seen once, igore them going forward."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "945b0b6b-65a3-464f-8233-f5d0be413745",
   "metadata": {},
   "outputs": [],
   "source": [
    "warnings.filterwarnings('ignore', category=FITSFixedWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdc0c322-5ab4-4241-a41c-0626f3348741",
   "metadata": {},
   "outputs": [],
   "source": [
    "h_indices = np.arange(16, dtype='int') + 1\n",
    "print(h_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f28a1be-f5db-445b-9df5-ca99bb6917d7",
   "metadata": {},
   "source": [
    "Print the boundary information for each segment.\n",
    "\n",
    " * `CRVAL` is the reference coordinate for the mosaic.\n",
    " * `CRPIX` are the pixel values relative to the reference.\n",
    " * `DATASEC` is the chip pixels that have data.\n",
    " * `DETSEC` defines where the segement is in the mosaic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe1658a9-e201-4a03-b520-275cf6de8c2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "for h in h_indices:\n",
    "    hdr = hdulist[h].header\n",
    "    wcs = WCS(hdr)\n",
    "    print('%2i %10s %6s   %5i %5i  %8.4f %8.4f %15s %15s' %\n",
    "          (h, hdr['EXTNAME'],\n",
    "           wcs.footprint_contains(targetCoord),\n",
    "           np.floor(hdr['CRPIX1']), np.floor(hdr['CRPIX2']),\n",
    "           hdr['CRVAL1'], hdr['CRVAL2'],\n",
    "           hdr['DATASEC'], hdr['DETSEC']))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9848c4a1-f82a-4904-917e-5c53714af5f7",
   "metadata": {},
   "source": [
    "Show the segment boundaries in a plot.\n",
    "\n",
    "Use [qualitative colormap](https://matplotlib.org/stable/users/explain/colors/colormaps.html) `tab20c`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7641137e-ab6c-486c-a63b-6c8d1ef0845f",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.colormaps[\"tab20c\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f041f04-4b71-4141-8a66-3122a7bb2bb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 6))\n",
    "for h in h_indices:\n",
    "    hdr = hdulist[h].header\n",
    "    ext = hdr['EXTNAME']\n",
    "    seg = ext[7:9]\n",
    "    x1, x2, y1, y2 = detsec_to_x1x2y1y2(hdr['DETSEC'])\n",
    "    xvals = [x1, x2, x2, x1, x1]\n",
    "    yvals = [y1, y1, y2, y2, y1]\n",
    "    plt.plot(xvals, yvals, lw=1, color=cmap(h))\n",
    "    plt.plot(x1, y1, 'o', ms=h+3, color='None', mec=cmap(h))\n",
    "    plt.plot(x2, y2, 's', ms=h+3, color='None', mec=cmap(h))\n",
    "    plt.text(x1-300, y1, seg, color=cmap(h))\n",
    "plt.xlabel('DETSEC X')\n",
    "plt.ylabel('DETSEC Y')\n",
    "plt.title('Segment Map')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f9c804d-7ae3-4829-bcb3-e76be9c63eeb",
   "metadata": {},
   "source": [
    "Figure 5: Boundaries of the 16 segments of the raw image, each in a different color.\n",
    "Squares and circle of ascending size mark the `DETSEC` corners."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a6c5975-3beb-4f3e-a248-20bdd84d8c7d",
   "metadata": {},
   "source": [
    "Create array `mosaic` to hold all segements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8581dd4-f303-4725-9b8d-cc7ba733351b",
   "metadata": {},
   "outputs": [],
   "source": [
    "mosaic = np.zeros((4000, 4072), dtype='float')\n",
    "print(np.shape(mosaic))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d099f627-bb9e-4278-ac68-a9a4fa3b6c32",
   "metadata": {},
   "source": [
    "Segment rows 1:2000 indexed at 0.\n",
    "\n",
    "Segment columns 4:512 indexed at 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df4f547f-80fa-4edf-904e-28fc57734b73",
   "metadata": {},
   "outputs": [],
   "source": [
    "seg_row = np.arange(2000, dtype='int')\n",
    "seg_col = np.arange(509, dtype='int') + 3"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf7cb9b3-8225-4ce9-aee2-aac1b15d253b",
   "metadata": {},
   "source": [
    "Paste each segment into the `mosaic`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5d09099-1a54-43b2-96cf-e17e22441050",
   "metadata": {},
   "outputs": [],
   "source": [
    "for h in h_indices:\n",
    "    hdr = hdulist[h].header\n",
    "    x1, x2, y1, y2 = detsec_to_x1x2y1y2(hdr['DETSEC'])\n",
    "    if y2 < y1:\n",
    "        mos_row = np.arange(start=y1-1, stop=y2-2, step=-1)\n",
    "    else:\n",
    "        mos_row = np.arange(start=y1-1, stop=y2, step=1)\n",
    "    if x2 < x1:\n",
    "        mos_col = np.arange(start=x1-1, stop=x2-2, step=-1)\n",
    "    else:\n",
    "        mos_col = np.arange(start=x1-1, stop=x2, step=1)\n",
    "\n",
    "    img_data = np.asarray(hdulist[h].data, dtype='float')\n",
    "    for y, n in zip(seg_row, mos_row):\n",
    "        for x, m in zip(seg_col, mos_col):\n",
    "            mosaic[n, m] = img_data[y, x]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b28da782-b462-4f85-a8e0-33d62b129d73",
   "metadata": {},
   "source": [
    "Get header, WCS, and define scale limits for image display."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7284b53-5f10-4664-b883-e64bd93107b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "img_hdr = hdulist[0].header\n",
    "img_wcs = WCS(img_hdr)\n",
    "zlimits = zscale.get_limits(mosaic)\n",
    "print(zlimits)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0aa1f2bb-020b-422f-82c9-aeb1466a0e17",
   "metadata": {},
   "source": [
    "**Note:** the header WCS is for a FITS-based origin and row-col,\n",
    "so when the pixels come out of the `img_wcs` here they're row-col (y,x) not col-row (x,y),\n",
    "so there's a flip that has to happen when plotting col-row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f91ea069-ea9d-4be5-8d8e-fa6274bcb37c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(6, 6))\n",
    "plt.subplot(projection=img_wcs)\n",
    "plt.imshow(mosaic, cmap='grey',\n",
    "           vmin=2250, vmax=2683)\n",
    "\n",
    "cra = 61.9699 * u.degree\n",
    "cdec = -36.8508 * u.degree\n",
    "offset = (90. / 3600.) * u.degree\n",
    "cPix = img_wcs.world_to_pixel(SkyCoord(cra, cdec, frame='icrs'))\n",
    "oraPix = img_wcs.world_to_pixel(SkyCoord(cra + offset, cdec, frame='icrs'))\n",
    "odecPix = img_wcs.world_to_pixel(SkyCoord(cra, cdec + offset, frame='icrs'))\n",
    "plt.plot([cPix[1], oraPix[1]], [cPix[0], oraPix[0]], \n",
    "         ls='solid', lw=1, color='cyan')\n",
    "plt.plot([cPix[1], odecPix[1]], [cPix[0], odecPix[0]], \n",
    "         ls='solid', lw=1, color='yellow')\n",
    "plt.text(oraPix[1], oraPix[0], 'E', color='cyan')\n",
    "plt.text(odecPix[1], odecPix[0], 'N', color='yellow')\n",
    "del cra, cdec, offset, cPix, oraPix, odecPix\n",
    "\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "plt.xlabel('DETSEC X')\n",
    "plt.ylabel('DETSEC Y')\n",
    "plt.axis('on')\n",
    "plt.grid(color='white', ls='solid')\n",
    "plt.title('Mosaic with 1,1 at lower left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eba3e65-a45e-4ac0-9b2e-8487687f1f0f",
   "metadata": {},
   "source": [
    "Figure 6: A mosaic of the raw image, with a compass showing north and east in the middle."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a480f12d-ead9-4c37-945f-02f9c78738c8",
   "metadata": {},
   "source": [
    "Save this file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b653d9d0-496d-41ec-9d46-55f5809747df",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = os.path.join(os.getenv(\"HOME\"), 'dp02_02c_image_raw.fits')\n",
    "urlretrieve(image_url, filename)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}