DP02_14_Injecting_Synthetic_Sources.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Injecting Synthetic Sources Into Single-Visit Images\n",
    "\n",
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\" alt=\"Rubin Observatory logo, a graphical representation of turning stars into data.\">\n",
    "<br>\n",
    "Contact author: Jeff Carlin <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: advanced <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9da1a210-d858-42fe-8591-570965b8be1a",
   "metadata": {},
   "source": [
    "**Description:** This tutorial demonstrates a method to inject artificial sources (stars and galaxies) into calexp images using the measured point-spread function of the given calexp image. Confirmation that the synthetic sources were correctly injected into the image is done by running a difference imaging task from the pipelines."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80a0baf5-51ad-40ec-8991-060a7b27c289",
   "metadata": {},
   "source": [
    "**Skills:** Use the `source_injection` tools to inject synthetic sources into images. Create a difference image from a `calexp` with injected sources."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393da88f-7978-4920-aa4a-a9830df6eed9",
   "metadata": {},
   "source": [
    "**LSST Data Products:** Butler calexp images and corresponding src catalogs, goodSeeingDiff_templateExp images, and injection_catalogs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c67fab9-136a-4adc-bb42-142b91ab69dd",
   "metadata": {},
   "source": [
    "**Packages:** astropy.io, astropy.coordinates, astropy.table, astropy.units, matplotlib, numpy, lsst.afw.display, lsst.afw.table, lsst.daf.base, lsst.daf.butler, lsst.ip.diffim, lsst.source.injection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f72b27f",
   "metadata": {},
   "source": [
    "**Credit:**\n",
    "Developed by Jeff Carlin in collaboration with Lee Kelvin and the Rubin Community Science Team. Much of the material is based on the documentation for the LSST [source_injection package](https://pipelines.lsst.io/v/daily/modules/lsst.source.injection/index.html)."
   ]
  },
  {
   "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>.\n",
    "Questions are welcome as new topics in the \n",
    "<a href=\"https://community.lsst.org/c/support/dp0\">Support - Data Preview 0 Category</a> \n",
    "of the Rubin Community Forum. \n",
    "Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc73be0",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notebook provides a brief introduction to the `source_injection` package from the LSST Science Pipelines. The `source_injection` tools can be used to inject sources into images produced during various stages of Pipelines processing, including visit-level images (e.g., `calexp`s), any dataset with a `datasetType` of `Exposure` (e.g., `postISRCCD` images), and coadd images (e.g., `deepCoadd`s). The main driver that both creates and injects synthetic sources into images is based on [Galsim](https://galsim-developers.github.io/GalSim/), so that the `source_injection` tools enable injection of many types of sources. These include a variety of parameterized galaxy models, stars, and postage stamp images.\n",
    "\n",
    "This notebook will teach some basics of using the `source_injection` package to insert artificial sources into images. Read the extensive [documentation](https://pipelines.lsst.io/v/daily/modules/lsst.source.injection/index.html) for more details and examples. Furthermore, the [Galsim](https://galsim-developers.github.io/GalSim/) documentation provides much more detail on the types of sources that can be injected."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc36f107",
   "metadata": {},
   "source": [
    "### 1.1 Import packages"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccb4fe61-835a-47ac-a691-e086ce766e79",
   "metadata": {},
   "source": [
    "<a href=\"https://numpy.org\">numpy</a> is a commonly used package for scientific computing with arrays in Python.\n",
    "\n",
    "<a href=\"https://matplotlib.org/\">matplotlib</a> is a library for creating static, animated, and\n",
    "interactive visualizations in Python \n",
    "(see also the <a href=\"https://matplotlib.org/stable/gallery/index.html\">matplotlib gallery</a>).\n",
    "\n",
    "<a href=\"https://www.astropy.org/\">astropy</a> is a package with many common astronomical utilities; here it is used for importing a FITS image, and working with coordinates and tables.\n",
    "\n",
    "From the `lsst` package, modules for accessing the butler, image display functions, and the source injection utilities are imported (see <a href=\"https://pipelines.lsst.io/\">pipelines.lsst.io</a> for documentation of the Science Pipelines codebase)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddc1458",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from astropy.io import fits\n",
    "import astropy.units as u\n",
    "from astropy.table import Table, vstack\n",
    "from astropy.coordinates import SkyCoord\n",
    "import getpass\n",
    "\n",
    "import lsst.afw.display as afwDisplay\n",
    "import lsst.afw.geom as afwGeom\n",
    "import lsst.afw.image as afwImage\n",
    "import lsst.afw.math as afwMath\n",
    "from lsst.daf.butler import Butler\n",
    "import lsst.geom as geom\n",
    "from lsst.source.injection import generate_injection_catalog\n",
    "from lsst.source.injection import VisitInjectConfig, VisitInjectTask\n",
    "from lsst.ip.diffim.subtractImages import AlardLuptonSubtractTask, AlardLuptonSubtractConfig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "662aa3ce-e9e5-4703-b3e3-15f554c7c325",
   "metadata": {},
   "source": [
    "### 1.2 Define functions and parameters\n",
    "\n",
    "Set `afwDisplay` to use `matplotlib`, and use colorblind-friendly colors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e75e7d7-4af9-4e2f-b25f-08b84d01823a",
   "metadata": {},
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib')\n",
    "plt.style.use('tableau-colorblind10')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1963fee-45c7-4689-852c-702f058af326",
   "metadata": {},
   "source": [
    "The following function is a modified version of [this code](https://github.com/lsst/atmospec/blob/1e7d6e8e5655cc13d71b21ba866001e6d49ee04e/python/lsst/atmospec/utils.py#L259-L301)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3885e7c5-7a1f-485e-bdcb-ed6be31aa44f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotate_exposure(exp, n_degrees):\n",
    "    \"\"\"Rotate an exposure by nDegrees clockwise.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    exp : `lsst.afw.image.exposure.Exposure`\n",
    "        The exposure to rotate\n",
    "    n_degrees : `float`\n",
    "        Number of degrees clockwise to rotate by\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    rotated_exp : `lsst.afw.image.exposure.Exposure`\n",
    "        A copy of the input exposure, rotated by nDegrees\n",
    "    \"\"\"\n",
    "    n_degrees = n_degrees % 360\n",
    "\n",
    "    wcs = exp.getWcs()\n",
    "\n",
    "    warper = afwMath.Warper('lanczos4')\n",
    "\n",
    "    affine_rot_transform = geom.AffineTransform.makeRotation(n_degrees*geom.degrees)\n",
    "    transform_p2top2 = afwGeom.makeTransform(affine_rot_transform)\n",
    "    rotated_wcs = afwGeom.makeModifiedWcs(transform_p2top2, wcs, False)\n",
    "\n",
    "    rotated_exp = warper.warpExposure(rotated_wcs, exp)\n",
    "    return rotated_exp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6e9358e-b3cb-4e87-84a6-c3b761ab4ab2",
   "metadata": {},
   "source": [
    "## 2. Get an image to inject sources into\n",
    "\n",
    "First, retrieve an image to inject sources into. To do so, instantiate a butler pointing to the DP0.2 dataset, identify `calexp` images overlapping a particular `tract`, and retrieve one of those images (plus some of its ancillary data) for later use.\n",
    "\n",
    "### 2.1 Instantiate a butler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4e00487-b9e6-4025-991d-533695c68dc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "butler_config = 'dp02'\n",
    "collections = '2.2i/runs/DP0.2'\n",
    "butler = Butler(butler_config, collections=collections)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a1115a8-4781-4f0b-9cb8-19ed4607d711",
   "metadata": {},
   "source": [
    "### 2.2 Identify images overlapping a particular tract and select one\n",
    "\n",
    "By default the next cell selects images from tract 3828, but you can change this to any tract number that exists in DP0.2 data. (See the map in Figure 15 of [\"The DC2 Simulated Sky Survey\"](https://ui.adsabs.harvard.edu/abs/2021ApJS..253...31L/abstract) overview paper; also visible in the [DP0.2 documentation](https://dp0-2.lsst.io/data-products-dp0-2/index.html#images)).\n",
    "\n",
    "Because a `calexp` image only spans one detector, select a single detector (in this case, number 19, but change this if you like) and find images from the chosen tract that contain that detector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15097b5f-5296-4530-86a3-acd5ffea0687",
   "metadata": {},
   "outputs": [],
   "source": [
    "tract = 3828\n",
    "where = f\"instrument='LSSTCam-imSim' AND skymap='DC2' AND \\\n",
    "          tract={tract} AND detector=19 AND band='g'\"\n",
    "calexp_g_DatasetRefs = sorted(list(set(butler.registry.queryDatasets('calexp', where=where))))\n",
    "print(f'Identified {len(calexp_g_DatasetRefs)} calexp DatasetRefs')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "217d700a-de66-488b-bdd6-ae90728993fb",
   "metadata": {},
   "source": [
    "Use index \"5\" to select an arbitrary `dataId`.\n",
    "Change this index to another integer to select a different image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12a63560-6f9c-4c21-954e-2163b8a6be43",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId_g = calexp_g_DatasetRefs[5].dataId\n",
    "print(f\"{dataId_g = }\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "485c5322-c0e0-42f4-b792-36c2e2dbefbe",
   "metadata": {},
   "source": [
    "#### 2.2.1 Retrieve the calexp image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "769b7130-0ec8-47c9-a8cb-765041e59bb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "calexp_g = butler.get('calexp', dataId=dataId_g)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dedd1af-f623-4e4c-b265-c3d663f016b0",
   "metadata": {},
   "source": [
    "#### 2.2.2 Retrieve additional information about the image\n",
    "\n",
    "To generate synthetic sources to be injected into the image, the coordinates of its bounding box are necessary.\n",
    "\n",
    "Retrieve the image's WCS and bounding box, and print its central coordinate to the screen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fe24c71-ca55-43df-92e0-b6da9964ce77",
   "metadata": {},
   "outputs": [],
   "source": [
    "wcs = calexp_g.getWcs()\n",
    "bbox = calexp_g.getBBox()\n",
    "print('bounding box: ', bbox)\n",
    "\n",
    "boxcen = bbox.getCenter()\n",
    "cen = wcs.pixelToSky(boxcen)\n",
    "sc_cen = SkyCoord(ra=cen[0].asDegrees()*u.deg, dec=cen[1].asDegrees()*u.deg)\n",
    "\n",
    "print(sc_cen)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dde6d3b1-0b90-4f2b-a272-e6f608cadcc1",
   "metadata": {},
   "source": [
    "#### 2.2.3 Figure out how large the image is on the sky\n",
    "\n",
    "Note that above the bounding box is roughly 4000 x 4000 pixels.\n",
    "\n",
    "Use the `bbox.getDimensions` and `wcs.getPixelScale` methods to estimate how large this image is in sky coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76060374-b3e1-475f-bf2f-6e8442c0b08b",
   "metadata": {},
   "outputs": [],
   "source": [
    "imsize = bbox.getDimensions()[0]*wcs.getPixelScale().asDegrees()\n",
    "print('Size of calexp in degrees: ', imsize)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96734edf-8d93-49f3-8ac0-1149ae835a4b",
   "metadata": {},
   "source": [
    "## 3. Make a catalog of synthetic sources\n",
    "\n",
    "Having retrieved a `calexp` image to inject into, now set up a simple synthetic source catalog.\n",
    "\n",
    "### 3.1 Make a catalog of galaxies and stars\n",
    "\n",
    "The `source_injection` package provides tools to create catalogs of synthetic sources. Here, use the `generate_injection_catalog` function to create the sources to inject.\n",
    "\n",
    "Note that the \"inject_size\" is selected to be 0.1 degrees, or slightly smaller than the size of the image as determined above (inject_size is a radius, so it equals 0.2/2 degrees)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28f7afc1-b476-4d42-89bf-117b0030817d",
   "metadata": {},
   "outputs": [],
   "source": [
    "inject_size = 0.2/2"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "91b9231d-a926-4c16-a491-313a0c3c5f75",
   "metadata": {},
   "source": [
    "#### 3.1.1 Make a catalog of galaxies to inject\n",
    "\n",
    "The simplest form of a galaxy in Galsim is parameterized by a [Sersic model](https://en.wikipedia.org/wiki/S%C3%A9rsic_profile): $I(r) = I_e~{\\rm exp}\\{-b_n [(\\frac{r}{r_e})^{1/n}-1]\\}$, which defines the shape of the galaxy's light profile as a function of radius (r) in terms of the \"Sersic index\" (n) and the \"half-light radius\" ($r_e$). (Note that $n = 1$ corresponds to an exponential profile.)\n",
    "\n",
    "The above equation results in a circular galaxy. To create elongated (elliptical) Sersic profiles, additionally specify the minor-to-major axis ratio (q) with a rotation angle (beta).\n",
    "\n",
    "`generate_injection_catalog` will automatically generate sources with all possible permutations of the parameters you provide. For example, the cell below specifies \"number=2\" to create 2 synthetic galaxies, but then specifies a single magnitude (mag), 3 values of Sersic index (n), 3 values of axis ratio (q), 2 values of beta, and a single value for half_light_radius. This should result in $2*1*3*3*2*1 = 36$ different combinations of those parameters.\n",
    "\n",
    "Finally, provide minimum and maximum RA and Dec coordinates in degrees. In this case, `generate_injection_catalog` will randomly select positions within those limits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bafc845-2fc1-42bf-a5e4-cce3330ecf52",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_injection_catalog_galaxies = generate_injection_catalog(\n",
    "    ra_lim=[sc_cen.ra.value-inject_size, sc_cen.ra.value+inject_size],\n",
    "    dec_lim=[sc_cen.dec.value-inject_size, sc_cen.dec.value+inject_size],\n",
    "    number=2,\n",
    "    seed='3210',\n",
    "    source_type=\"Sersic\",\n",
    "    mag=[15.0],\n",
    "    n=[1, 2, 4],\n",
    "    q=[0.9, 0.5, 0.1],\n",
    "    beta=[31.0, 144.0],\n",
    "    half_light_radius=[15.0],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2004011-6b1e-42f6-b910-1d33c739f2a2",
   "metadata": {},
   "source": [
    "#### 3.1.2 Make a catalog of stars to inject\n",
    "\n",
    "Using a similar method, create a catalog of stars to be injected. Specify 5 values of magnitude, with 7 instances, which will result in 35 stars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27c2e170-604c-41c3-8ca7-2ea45de324fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_injection_catalog_stars = generate_injection_catalog(\n",
    "    ra_lim=[sc_cen.ra.value-inject_size, sc_cen.ra.value+inject_size],\n",
    "    dec_lim=[sc_cen.dec.value-inject_size, sc_cen.dec.value+inject_size],\n",
    "    number=7,\n",
    "    seed='432',\n",
    "    source_type=\"Star\",\n",
    "    mag=[16.0, 17.0, 18.0, 19.0, 20.0],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e23a2012-fa6f-41ba-8df6-993e5e76d64b",
   "metadata": {},
   "source": [
    "#### 3.1.3 Combine the galaxies and stars into a single synthetic source catalog\n",
    "\n",
    "In particular, use the [vstack](https://docs.astropy.org/en/stable/table/operations.html#id4) method from Astropy's [Table](https://docs.astropy.org/en/stable/table/index.html) class to combine them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f31a1ad1-1a5f-409b-aaec-4d8baa61a683",
   "metadata": {},
   "outputs": [],
   "source": [
    "inject_cat = vstack([my_injection_catalog_galaxies, my_injection_catalog_stars])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f89ee26c-6d5c-4118-bd92-4f150dc886bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "inject_cat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e96a109c-fb7e-4158-a1c5-e3983edf63dd",
   "metadata": {},
   "source": [
    "### 3.2 Add a postage stamp of an image to the injection catalog\n",
    "\n",
    "In addition to parameterized sources of many types, a postage-stamp image (for example, an output image from a simulation) can also be injected. This can be any image; the only requirement is that it must be in FITS format and have a valid WCS (see a fun example [here](https://pipelines.lsst.io/v/d_2024_02_27/modules/lsst.source.injection/faqs.html#frequently-asked-questions-and-answers-about-synthetic-source-injection)). This notebook demonstrates injection of an [SDSS g-band image](https://ned.ipac.caltech.edu/byname?objname=NGC+4152&hconst=67.8&omegam=0.308&omegav=0.692&wmap=4&corr_z=1) of the galaxy PGC 038749 downloaded from the NASA/IPAC Extragalactic Database (NED).\n",
    "\n",
    "Load the image, then display it to see what it depicts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6e73816-f9fc-4a8f-97a5-f62a41d4db17",
   "metadata": {},
   "outputs": [],
   "source": [
    "stamp_img_hdu = fits.open('data/PGC_038749_I_g_bbl2011.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "681bacda-8908-4b1c-ba9f-94acba452a82",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.subplot()\n",
    "im = plt.imshow(stamp_img_hdu[0].data, cmap='gray', vmin=-20.0, vmax=100, origin='lower')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d12fcf9-6e9c-48be-afd8-27eee0524263",
   "metadata": {},
   "source": [
    "> Figure 1: The SDSS g-band image of PGC 038749 displayed in grayscale, with axes in pixel coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d09ab430-d2ac-4f76-98e5-5cb7bb77d226",
   "metadata": {},
   "source": [
    "#### 3.2.1 Create an entry in the injection catalog for the postage stamp image\n",
    "\n",
    "Add a postage stamp entry into the injection catalog, using `vstack` to add an astropy `Table` containing a single row. Specify its position, magnitude, and \"source_type = Stamp,\" in addition to the filename of the image to inject. (It is OK to leave unnecessary columns empty.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c914dbd4-80e5-48c5-a396-20ed199c0ee5",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_injection_catalog_stamp = Table(\n",
    "    {\n",
    "        'injection_id': [9999],\n",
    "        'ra': [56.622],\n",
    "        'dec': [-36.488],\n",
    "        'source_type': ['Stamp'],\n",
    "        'mag': [14.8],\n",
    "        'stamp': ['data/PGC_038749_I_g_bbl2011.fits'],\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6ecfcdd-b00e-4860-be58-d8b5acd2b25b",
   "metadata": {},
   "outputs": [],
   "source": [
    "inject_cat = vstack([inject_cat, my_injection_catalog_stamp])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd751614-9424-4ebf-9165-39cec1773fdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "inject_cat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca79ab5c-b8c1-4768-bbbe-cab766300df9",
   "metadata": {},
   "source": [
    "#### 3.2.2 Create a rotated version of the stamp image, and add an entry for it into the injection catalog\n",
    "\n",
    "To apply a rotation to the postage stamp image before injecting it into the `calexp`, use the `rotateExposure` function created above.\n",
    "\n",
    "The process is to read the image, apply a rotation to it, write it out as a new image, and then follow the same steps as the previous section to add it to the injection catalog.\n",
    "\n",
    "First, read the postage stamp image using the `afwImage` package. It's OK to ignore the warning about an expected extension type that will pop up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10d00f8c-ab25-4478-9f00-d23c454d2f11",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "stamp_img_orig = afwImage.ExposureF.readFits('data/PGC_038749_I_g_bbl2011.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51ae9220-d6e4-466b-9ef1-e1b1115a9f66",
   "metadata": {},
   "source": [
    "Apply a rotation to the image using the `rotateExposure` function defined above. This will apply a rotation to the WCS, and then \"warp\" the image (i.e., resample it) onto the new, rotated WCS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0fa1c3e-af7c-422d-9a94-20934fceb6e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "rotation_angle = 57.0\n",
    "stamp_img_rotated = rotate_exposure(stamp_img_orig, rotation_angle)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2369cc9-dcbd-4012-bc43-d85dd8fc9f70",
   "metadata": {},
   "source": [
    "Confirm that the image has been rotated by displaying them side-by-side:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71ef330e-853d-42c8-b15e-84b0d0343f32",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(6, 4), dpi=150)\n",
    "\n",
    "plt.sca(ax[0])\n",
    "display0 = afwDisplay.Display(frame=fig)\n",
    "display0.scale('linear', min=-20, max=150)\n",
    "display0.mtv(stamp_img_orig.image)\n",
    "plt.title('original stamp image')\n",
    "\n",
    "plt.sca(ax[1])\n",
    "display1 = afwDisplay.Display(frame=fig)\n",
    "display1.scale('linear', min=-20, max=150)\n",
    "display1.mtv(stamp_img_rotated.image)\n",
    "plt.title('rotated stamp image')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e3eae25-5043-4554-8fcf-83d8b7bea68f",
   "metadata": {},
   "source": [
    "> Figure 2: At left, the same image as shown in Figure 1.\n",
    "> At right, the image rotated by 57 degrees clockwise (west of north), as if the telescope had a rotation angle of 57 degrees east of north."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a18309c-5ab6-4d22-a40f-60c87f801f0b",
   "metadata": {},
   "source": [
    "When the image is rotated, the \"empty\" pixels are given NaN values. Replace those with zeros (if you leave them as NaNs, the image values will not be read correctly when reading in the resulting image)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a84caa50-f84f-4013-92a3-5d32e1a0b7d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "stamp_img_rotated.image.array[np.where(np.isnan(stamp_img_rotated.image.array))] = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3e476fa-b40e-4d3b-a9c8-f0d0a05f99c1",
   "metadata": {},
   "source": [
    "Write the rotated stamp to an output file in the home directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64f29cde-d9cc-42a5-9581-695225000097",
   "metadata": {},
   "outputs": [],
   "source": [
    "home_directory = '/home/' + getpass.getuser() + '/'\n",
    "stamp_img_rotated.writeFits(home_directory+'stamp_rotated.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5987196-89b5-4b10-91c4-445b86c4275f",
   "metadata": {},
   "source": [
    "Add an entry to the injection catalog for the rotated stamp."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f693331-af36-4d18-b280-b8b39eccc970",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_injection_catalog_rotated_stamp = Table(\n",
    "    {\n",
    "        'injection_id': [99999],\n",
    "        'ra': [56.642],\n",
    "        'dec': [-36.488],\n",
    "        'source_type': ['Stamp'],\n",
    "        'mag': [14.8],\n",
    "        'stamp': [home_directory+'stamp_rotated.fits'],\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6311ed07-9674-43e3-a2cf-d33f3630f09d",
   "metadata": {},
   "outputs": [],
   "source": [
    "inject_cat = vstack([inject_cat, my_injection_catalog_rotated_stamp])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02faa0d1-dd40-4522-9812-6ec90485157e",
   "metadata": {},
   "source": [
    "## 4. Inject sources into a calexp image\n",
    "\n",
    "### 4.1 Use source_injection tools to inject sources\n",
    "\n",
    "In Section 3, a catalog specifying the synthetic sources to inject was created.\n",
    "\n",
    "Now run the task from `source_injection` that injects sources into the `calexp` image that was retrieved earlier.\n",
    "\n",
    "First, extract the point spread function (PSF), photometric calibration object, and the WCS (World Coordinate System) object associated with the `calexp` image. These will be passed to the injection task so that sources can be injected using the properties of the image itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbcabdf1-aa2f-49ca-b686-b2e321832fb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "psf = calexp_g.getPsf()\n",
    "photo_calib = calexp_g.getPhotoCalib()\n",
    "wcs = calexp_g.getWcs()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "435f4b84-f7ce-4262-9b50-9d1afbefa94b",
   "metadata": {},
   "source": [
    "### 4.2 Instantiate the injection classes\n",
    "\n",
    "First, instantiate the `VisitInjectConfig` class. The `VisitInjectConfig` class is where configuration of the injection task occurs, allowing for modifications to be made to how the task operates.\n",
    "\n",
    "Then instantiate the `VisitInjectTask`, using `inject_config` as the configuration argument.\n",
    "\n",
    "NOTE: For injections into other dataset types, use the appropriate option from the following list:\n",
    "```\n",
    "from lsst.source.injection import CoaddInjectConfig, CoaddInjectTask\n",
    "from lsst.source.injection import ExposureInjectConfig, ExposureInjectTask\n",
    "from lsst.source.injection import VisitInjectConfig, VisitInjectTask\n",
    "```\n",
    "Both `ExposureInject` and `VisitInject` accept dimensions of `(\"instrument\", \"visit\", \"detector\")`. They differ in that `ExposureInject` takes a `postISRCCD` exposure (i.e., a raw image that has had \"instrument signature removal,\" or ISR, applied, and no further processing) as an input, while `VisitInject` expects to operate on a `calexp` (an image that has been astrometrically and photometrically calibrated to an external source)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c64c1e6-02ca-41c5-8d52-f416352d8b7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "inject_config = VisitInjectConfig()\n",
    "\n",
    "inject_task = VisitInjectTask(config=inject_config)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "115d7b80-d060-499a-8f63-ce733fc63b6b",
   "metadata": {},
   "source": [
    "### 4.3 Inject sources into an image\n",
    "\n",
    "Next, execute the `run` method of the `inject_task`.\n",
    "\n",
    "The `run` method needs the following inputs:\n",
    "\n",
    "* the input injection catalog\n",
    "* the input `calexp`\n",
    "* the PSF of the input exposure\n",
    "* the WCS information\n",
    "* the photometric calibration information\n",
    "\n",
    "The PSF, WCS, and photo_calib inputs were already retrieved directly from the `calexp` they are associated with.\n",
    "\n",
    "The inject task provides two outputs:\n",
    "\n",
    "* the output exposure with sources injected\n",
    "* the output source injection catalog\n",
    "\n",
    "The output source injection catalog is identical to the input, but with two additional columns (x and y) denoting the pixel coordinates of these sources. Note that this catalog is NOT the science catalog containing the full suite of LSST Science Pipelines outputs. To get that, this source injected image will need to be processed by additional Science Pipelines tasks.\n",
    "\n",
    "> Note: it is best to use a clone of the input `calexp`. This is because the `calexp` that is input for source injection will be edited in-place. Inputting a clone makes it possible to compare the output image to the original `calexp` later in this notebook.\n",
    "\n",
    "#### 4.3.1 Run the source injection task\n",
    "Run the source injection task and extract the \"output_exposure\" and \"output_catalog\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2da88cf6-bfda-4b3f-ad37-932ab896f7ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "injected_output = inject_task.run(\n",
    "    injection_catalogs=inject_cat,\n",
    "    input_exposure=calexp_g.clone(),\n",
    "    psf=psf,\n",
    "    photo_calib=photo_calib,\n",
    "    wcs=wcs,\n",
    ")\n",
    "injected_exposure = injected_output.output_exposure\n",
    "injected_catalog = injected_output.output_catalog"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f6029cf-fbee-4b54-a74c-2738bfd46797",
   "metadata": {},
   "source": [
    "#### 4.3.2 Compare the images before and after injection\n",
    "\n",
    "Plot the images side-by-side to confirm that the source injection task successfully added sources to the image.\n",
    "\n",
    "Note that the injected image will have different flux values than the original, which would cause the image scaling of the two images to be slightly different by default (to confirm this, uncomment the \"display0.scale('linear', 'zscale')\" lines below and comment out the lines below them that explicitly set the min/max values). For a direct comparison, explicitly set the minimum and maximum pixel values of the color scale.\n",
    "\n",
    "Additionally, to zoom in on the injected postage stamp to see how it looks in the image, uncomment the line below where it says \"To zoom on the PGC 038749 stamp:\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e96452d-b815-4355-8de8-22983d216804",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_injected_calexp = injected_exposure.clone()\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10, 6), dpi=150)\n",
    "\n",
    "plt.sca(ax[0])\n",
    "display0 = afwDisplay.Display(frame=fig)\n",
    "# display0.scale('linear', 'zscale')\n",
    "display0.scale('linear', min=-20, max=150)\n",
    "display0.mtv(calexp_g.image)\n",
    "plt.title('calexp image')\n",
    "\n",
    "plt.sca(ax[1])\n",
    "display1 = afwDisplay.Display(frame=fig)\n",
    "# display1.scale('linear', 'zscale')\n",
    "display1.scale('linear', min=-20, max=150)\n",
    "display1.mtv(plot_injected_calexp.image)\n",
    "# To zoom on the PGC 038749 stamp:\n",
    "# display1.mtv(plot_injected_calexp.image[3300:3800, 1500:2000])\n",
    "plt.title('injected_calexp image')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d10f921-4f59-47d3-8bb3-5b13cbdcd33f",
   "metadata": {},
   "source": [
    "> Figure 3: The `calexp` image displayed in grayscale (with scale bar at right and axes in pixels) before (left) and after (right) synthetic source injection."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2a56083-87b6-4b4a-a2be-5f9373d0408b",
   "metadata": {},
   "source": [
    "### 4.4 Create a difference image to see the injected sources\n",
    "\n",
    "One reason for injecting sources into a `calexp` image might be to test whether they would be recovered in a difference image. This section demonstrates the use of an image differencing task from the LSST Science Pipelines to create a difference image.\n",
    "\n",
    "#### 4.4.1 Load the template image\n",
    "Template images are created as coadds of the images with the best seeing, and are called `goodSeeingDiff_templateExp`. Load the appropriate template corresponding to the `calexp` used throughout this notebook (note the use of the same dataId used to load the `calexp`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72da2718-71f0-4653-a7f4-228267be99a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "difftemp = 'goodSeeingDiff_templateExp'\n",
    "templateExposure = butler.get(difftemp, dataId=dataId_g.required)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fc58a81-8230-49dc-926d-fae60a462f65",
   "metadata": {},
   "source": [
    "#### 4.4.2 Initialize and run the image subtraction task\n",
    "\n",
    "The first two lines below load the default configuration for the `AlardLuptonSubtract` task, and then initialize the task with that configuration. The task requires (1) a template exposure, (2) the science exposure, and (3) the catalog of sources from the science exposure.\n",
    "\n",
    "Load the source catalog (`src`) and run the task. Override the default configuration to point the \"sourceSelector\" to \"base_ClassificationExtendedness_value.\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ddb4207-6ee0-4687-b7f9-df9e1c1db690",
   "metadata": {},
   "source": [
    "As of w_2024_38, the source selector within `AlardLuptonSubtractTask` has been changed. In the new selector, the default value that is used does not exist in DP0.2 catalogs, so we need to override the configured value to use base_ClassificationExtendedness_value instead. However, that override will fail on all pipeline versions prior to w_38, so we wrap it in a try:except block so that this cell will run with any pipelines version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03039124-03de-49b3-b998-e59ae1a74a7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "config = AlardLuptonSubtractConfig()\n",
    "try:\n",
    "    config.sourceSelector.value.unresolved.name = 'base_ClassificationExtendedness_value'\n",
    "except:\n",
    "    pass\n",
    "\n",
    "alTask = AlardLuptonSubtractTask(config=config)\n",
    "\n",
    "scienceExposure = injected_exposure\n",
    "\n",
    "sources = butler.get('src', dataId=dataId_g)\n",
    "result = alTask.run(templateExposure, scienceExposure, sources)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e23e0b4d-4adf-4a62-8aa7-3118ad6ea248",
   "metadata": {},
   "source": [
    "#### 4.4.3 Compare the original and the difference image\n",
    "\n",
    "Plot the original `calexp`, the injected image, and the difference image side by side.\n",
    "\n",
    "Set `xmin`, `xmax`, `ymin`, and `ymax` to zoom in on an arbitrary section of the image.\n",
    "Change the values to zoom in on a different region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5958c29-554b-4497-a6df-e301f0856be5",
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, xmax = 500, 2000\n",
    "ymin, ymax = 1000, 2500\n",
    "\n",
    "fig, ax = plt.subplots(1, 3, figsize=(14, 11))\n",
    "plt.sca(ax[0])  # set the first axis as current\n",
    "display1 = afwDisplay.Display(frame=fig)\n",
    "display1.scale('linear', 'zscale')\n",
    "display1.mtv(calexp_g.image)\n",
    "plt.title('coadd image (template)')\n",
    "plt.xlim(xmin, xmax)\n",
    "plt.ylim(ymin, ymax)\n",
    "plt.sca(ax[1])  # set the second axis as current\n",
    "display2 = afwDisplay.Display(frame=fig)\n",
    "display2.scale('linear', 'zscale')\n",
    "display2.mtv(plot_injected_calexp.image)\n",
    "plt.title('calexp with synthetic sources')\n",
    "plt.xlim(xmin, xmax)\n",
    "plt.ylim(ymin, ymax)\n",
    "plt.sca(ax[2])  # set the third axis as current\n",
    "display3 = afwDisplay.Display(frame=fig)\n",
    "display3.scale('linear', 'zscale')\n",
    "display3.mtv(result.difference.image)\n",
    "plt.xlim(xmin, xmax)\n",
    "plt.ylim(ymin, ymax)\n",
    "plt.title('diffim')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0121d5e0-8d4e-49dc-9597-f035fa20c893",
   "metadata": {},
   "source": [
    "> Figure 4: At left, \"zoomed-in\" cutout of the `deepCoadd` image displayed in grayscale (with scale bar at right and axes in pixels). In the center, the cutout of the `calexp` with synthetic sources injected (as in the right panel of Figure 3). At right, the difference image which results from subtracting the `coadd` from the `calexp`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d5af646-bb96-48c1-9343-8979f93fa40d",
   "metadata": {},
   "source": [
    "# 5. Exercises for the learner\n",
    "\n",
    "Further explorations could include:\n",
    "\n",
    "1.  Injecting sources into coadd images.\n",
    "2.  Exploring other types of objects that can be injected into images (i.e., other Galsim parameters for more complicated sources).\n",
    "3.  Injecting variable objects into multiple calexp images and testing their recoverability.\n",
    "4.  Running detection and measurement tasks from the LSST Science Pipelines to test the recovery and measurement accuracy of injected sources (e.g., following DP0.2 tutorial [Notebook 05: Introduction to Source Detection](https://github.com/rubin-dp0/tutorial-notebooks/blob/main/DP02_05_Introduction_to_Source_Detection.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "098ae337-1e58-409c-88ba-114ca4ba3a1b",
   "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
}