Address comments by CW

Author: plazas

DP1/100_How_to_Use_RSP_Tools/103_Image_access_and_display/103_7_Big_deep_coadd_cutout.ipynb

@@ -68,8 +68,6 @@
    "source": [
     "### 1.1. Import packages\n",
     "\n",
-    "Geometry classes from `lsst.geom` handle pixel and sky coordinates: `Box2I`, `Point2I`, and `Extent2I` define integer-based image regions, `Point2D` represents floating-point coordinates, and `SpherePoint` with `degrees` specifies positions on the celestial sphere.\n",
-    "`lsst.sphgeom.Box` defines regions on the sky.\n",
     "`GetTemplateTask` (`lsst.ip.diffim`) retrieves template images for image differencing, while `lsst.pipe.tasks.makeDiscreteSkyMap` generates discrete sky maps for defining sky coverage."
    ]
   },
@@ -152,6 +150,7 @@
    "metadata": {},
    "source": [
     "Use the galaxy cluster Abell 360 (A360), a rich galaxy cluster at redshift 0.22 that is included in the DP1 data release in the `SV_38_7` field.\n",
+    "The angular extent of A360 is larger than a single coadd patch, so creating a large cutout is necessary to capture the full cluster within a single image.\n",
     "\n",
     "Define the coordinates of the brightest cluster galaxy (BCG), in degrees."
    ]
@@ -164,7 +163,7 @@
    "outputs": [],
    "source": [
     "ra_bcg = 37.865\n",
-    "dec_bcg = 6.98"
+    "dec_bcg = 6.982"
    ]
   },
   {
@@ -219,7 +218,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "bands = \"('r')\""
+    "my_band = ['r']"
    ]
   },
   {
@@ -232,7 +231,7 @@
     "dataset_refs = set(\n",
     "    butler.query_datasets(\n",
     "        \"deep_coadd\",\n",
-    "        where=f\"patch.region OVERLAPS(region) and band in {bands}\",\n",
+    "        where=f\"patch.region OVERLAPS(region) and band in ('{my_band[0]}')\",\n",
     "        bind={\"region\": box_region}\n",
     "    )\n",
     ")"
@@ -256,14 +255,38 @@
     "dataset_refs"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "cb4c5d01-ffee-4029-b1f7-b16eb7b7a56d",
+   "metadata": {},
+   "source": [
+    "## 3. Create the cutout"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6cf9f404-8f17-4d21-85d9-6c2bcfd9c014",
+   "metadata": {},
+   "source": [
+    "Create an image of an arbitrary sky region that is not aligned with the native Legacy Survey of Space and Time (LSST) sky tessellation.\n",
+    "This is similar to the process in alert production, where past LSST observations are combined to generate static-sky templates for detecting transients in new exposures.\n",
+    "The same template-generation code (`GetTemplateTask`) can be used to transform input data aligned on the LSST tessellation into the required geometry.\n",
+    "\n",
+    "To use the template code, a single astrometric solution or World Coordinate System (WCS) that encompasses all input data is needed, instead of having a separate tangent-plane WCS for each tract.\n",
+    "A new and custom discrete skymap is therefore created for the input data, producing one tract with a single WCS projection.\n",
+    "A bounding box is then defined in these new coordinates, centered on the brightest cluster galaxy (BCG) and spanning the specified width.\n",
+    "\n",
+    "With this setup, images are generated per band by combining the inputs that fall within the bounding box. The `GetTemplateTask` creates coadds for each band."
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "b6e864a7-47cf-4b58-8d87-d22874d71b80",
    "metadata": {},
    "source": [
-    "Prepare the inputs for building a custom skymap.\n",
+    "Prepare the inputs for building the custom skymap.\n",
     "\n",
-    "For each dataset reference, the WCS and bounding box of its patch are stored, so the collection of positions can later define a skymap that includes only the data in use.\n",
+    "For each dataset reference, the WCS and bounding box of its patch are stored, so the collection of positions can later define the custom skymap that includes only the data in use.\n",
     "At the same time, the dataset references are organized into a nested dictionary: first grouped by band, then within each band grouped by tract.\n",
     "Each tract key contains all the dataset refs for that tract band combination."
    ]
@@ -288,30 +311,6 @@
     "    tract_container.append(ds_ref)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "id": "cb4c5d01-ffee-4029-b1f7-b16eb7b7a56d",
-   "metadata": {},
-   "source": [
-    "## 3. Create the cutout"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "6cf9f404-8f17-4d21-85d9-6c2bcfd9c014",
-   "metadata": {},
-   "source": [
-    "Create an image of an arbitrary sky region that is not aligned with the native Legacy Survey of Space and Time (LSST) sky tessellation.\n",
-    "This is similar to the process in alert production, where past LSST observations are combined to generate static-sky templates for detecting transients in new exposures.\n",
-    "The same template-generation code (`GetTemplateTask`) can be used to transform input data aligned on the LSST tessellation into the required geometry.\n",
-    "\n",
-    "To use the template code, a single astrometric solution or World Coordinate System (WCS) that encompasses all input data is needed, instead of having a separate tangent-plane WCS for each tract.\n",
-    "A discrete skymap is therefore created for the input data, producing one tract with a single WCS projection.\n",
-    "A bounding box is then defined in these new coordinates, centered on the brightest cluster galaxy (BCG) and spanning the specified width.\n",
-    "\n",
-    "With this setup, images are generated per band by combining the inputs that fall within the bounding box. The `GetTemplateTask` creates coadds for each band."
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "0fea3ad8-0fe7-45da-808a-a49373058101",
@@ -335,15 +334,16 @@
    "id": "b7ce29de-5122-4939-9a46-9550596ef8a3",
    "metadata": {},
    "source": [
-    "Override the default configuration values for LSSTComCam, in particular, the pixel scale and the WCS projection."
+    "The default configuration values of this custom skymap for the pixel scale and WCS projection type are not appropriate for LSSTComCam.\n",
+    "The following cells override these parameters with appropriate values."
    ]
   },
   {
    "cell_type": "markdown",
    "id": "549fa673-002d-4f6e-8ab9-7027a00ccd13",
    "metadata": {},
    "source": [
-    "Use the BCG location to obtain the pixel scale.\n",
+    "Use the BCG location to obtain the pixel scale there.\n",
     "First, get the WCS of the patch and tract where the BCG is located.\n",
     "From the 300 DP1 tutorials, the patch and tract numbers are 10463 and 62 respectively."
    ]
@@ -384,7 +384,7 @@
    "id": "58613b7a-9f1f-4156-b8ba-cb08c915996d",
    "metadata": {},
    "source": [
-    "Get the WCs from the coadd."
+    "Get the WCS from the coadd."
    ]
   },
   {
@@ -413,7 +413,7 @@
    "id": "f1bb7b34-66e6-4bbf-9381-a66a2bcc021c",
    "metadata": {},
    "source": [
-    "Use the plate sclate at the BCG location."
+    "Use the plate scale at the BCG location."
    ]
   },
   {
@@ -449,7 +449,7 @@
    "id": "b63a608c-3646-4913-ada2-5c961f82b0f7",
    "metadata": {},
    "source": [
-    "Set the WCS projection.\n",
+    "Set the WCS projection type.\n",
     "The WCS of `deep_coadd` images are `TAN` projections ([see the DP1 documentation](https://dp1.lsst.io/products/images/deep_coadd.html#wcs))."
    ]
   },
@@ -537,48 +537,46 @@
    "id": "f6bd5018-ed88-4df7-abfc-deeb3b17091a",
    "metadata": {},
    "source": [
-    "Loop over bands (just the `r` band in this case) to create images that are combinations of the inputs contained within the defined box. Note that the run method modifies the bounding box in place, so a copy is made."
+    "Create a new coadd image in the band defined above (`r`) that is a combination of the inputs contained within the defined box.\n",
+    "Note that the run method modifies the bounding box in place, so a copy is made."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "683bdacb-1fda-4588-96c7-dca359d0ea36",
+   "id": "18354f90-d2fe-4659-b021-ff212ba7380a",
    "metadata": {},
    "outputs": [],
    "source": [
     "getTemplateTask = GetTemplateTask()\n",
-    "coadds = {}\n",
-    "for band, mapping in sorted_data.items():\n",
-    "    small_box_copy = Box2I(small_box)\n",
-    "    new_image_tmp = getTemplateTask.run(\n",
-    "        coaddExposures={\n",
-    "            tract: [butler.get(ref) for ref in refs] for tract, refs in mapping.items()\n",
-    "        },\n",
-    "        bbox=small_box_copy,\n",
-    "        wcs=new_wcs,\n",
-    "        dataIds={tract: [ref.dataId for ref in refs] for tract, refs in mapping.items()},\n",
-    "        physical_filter=band,\n",
-    "    )\n",
-    "    coadds[band] = new_image_tmp"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "6870187a-7aba-4863-8570-00f63ad7c5fe",
-   "metadata": {},
-   "source": [
-    "## 4. Visualize the custom cutout"
+    "mapping = sorted_data[my_band[0]]\n",
+    "\n",
+    "small_box_copy = Box2I(small_box)\n",
+    "new_image_r = getTemplateTask.run(\n",
+    "    coaddExposures={\n",
+    "        tract: [butler.get(ref) for ref in refs] for tract, refs in mapping.items()\n",
+    "    },\n",
+    "    bbox=small_box_copy,\n",
+    "    wcs=new_wcs,\n",
+    "    dataIds={tract: [ref.dataId for ref in refs] for tract, refs in mapping.items()},\n",
+    "    physical_filter=my_band[0],\n",
+    ")"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "672f1c76-aaba-46e6-882a-162e60fca0de",
+   "id": "79b9046c-f83d-4d33-ac7e-5d0f165468d8",
    "metadata": {},
    "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6870187a-7aba-4863-8570-00f63ad7c5fe",
+   "metadata": {},
    "source": [
-    "new_image_r = coadds['r']"
+    "## 4. Visualize the custom cutout"
    ]
   },
   {
@@ -621,6 +619,16 @@
     "> Figure 1: Cutout of the Abell 360 cluster, displayed using `matplotlib` as backend for `awf.display`. Note that the coordinates are in pixel units."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eb0982f3-7990-43f4-be0b-edd43cc35a28",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "new_image_r.template.image.array.shape"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "aa61bcdb-a195-4b98-a568-71bce6a2d5fe",
@@ -717,7 +725,7 @@
     "## 5. Compare PSFs\n",
     "\n",
     "Retrieve the Point Spread Function (PSF) objects for the new big coadd and for the coadd associated with the cluster BCG. \n",
-    "Compare the PSF size at the location of the BCG and at one corner of the big coadd region."
+    "Compare the PSF size at the location of the BCG and at near one corner of the big coadd region."
    ]
   },
   {
@@ -735,10 +743,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "margin = 0.01\n",
     "locations = [\n",
     "    (\"center\", bcg_point),\n",
-    "    (\"corner\", SpherePoint(ra_bcg - box_size_deg,\n",
-    "                           dec_bcg - box_size_deg, degrees)),\n",
+    "    (\"corner\", SpherePoint(ra_bcg - box_size_deg + margin,\n",
+    "                           dec_bcg - box_size_deg + margin, degrees)),\n",
     "]"
    ]
   },
@@ -747,7 +756,7 @@
    "id": "c9137e34-b303-4647-a0ef-8990012619cc",
    "metadata": {},
    "source": [
-    "Loop over locations and compare the PSF sizes."
+    "Loop over locations and compare the PSF Full Width at Half Maximum (FWHM) size for each one."
    ]
   },
   {
@@ -757,6 +766,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "SIGMA_TO_FWHM = 2.0*np.sqrt(2.0*np.log(2.0))\n",
     "for label, loc in locations:\n",
     "    center_pix_coadd = wcs_bcg.skyToPixel(loc)\n",
     "    point_coadd = Point2D(center_pix_coadd)\n",
@@ -769,8 +779,8 @@
     "    psf_sigma_new = newimage_psf.computeShape(point_new).getDeterminantRadius()\n",
     "\n",
     "    percent_diff = np.abs(psf_sigma_coadd - psf_sigma_new) / psf_sigma_new * 100\n",
-    "    print(f\"{label} coadd_bcg sigma = {psf_sigma_coadd: .3f} pixels\")\n",
-    "    print(f\"{label} newimage sigma = {psf_sigma_new: .3f} pixels\")\n",
+    "    print(f\"{label} coadd_bcg FWHM = {psf_sigma_coadd * SIGMA_TO_FWHM: .3f} pixels\")\n",
+    "    print(f\"{label} newimage FWHM = {psf_sigma_new * SIGMA_TO_FWHM: .3f} pixels\")\n",
     "    print(f\"{label} percent difference = {percent_diff: .2f}%\\n\")"
    ]
   }