From ef96bb74b42f755402efc944022bcf9bae615ae4 Mon Sep 17 00:00:00 2001 From: jrob93 Date: Fri, 28 Aug 2026 16:04:58 +0000 Subject: [PATCH] add fits reprojection --- .../201_Alerts/201_1_Alert_packets.ipynb | 241 +++++++++++++++++- 1 file changed, 234 insertions(+), 7 deletions(-) diff --git a/Prompt/200_Data_products/201_Alerts/201_1_Alert_packets.ipynb b/Prompt/200_Data_products/201_Alerts/201_1_Alert_packets.ipynb index 3ee7d420..fda0b407 100644 --- a/Prompt/200_Data_products/201_Alerts/201_1_Alert_packets.ipynb +++ b/Prompt/200_Data_products/201_Alerts/201_1_Alert_packets.ipynb @@ -22,8 +22,8 @@ "For the Rubin Science Platform at data.lsst.cloud.\\\n", "Data Release: [Prompt Products](https://prompt-products.lsst.io/)\\\n", "Container Size: Large\\\n", - "LSST Science Pipelines version: v30.0.11\\\n", - "Last verified to run: 2026-08-18\\\n", + "LSST Science Pipelines version: r29.2.0\\\n", + "Last verified to run: 2026-06-17\\\n", "Repository: [github.com/lsst/tutorial-notebooks](https://github.com/lsst/tutorial-notebooks)\\\n", "DOI: [10.11578/rubin/dc.20250909.20](https://doi.org/10.11578/rubin/dc.20250909.20)" ] @@ -119,6 +119,8 @@ " AsinhStretch, LinearStretch,\n", " ImageNormalize)\n", "from astropy.wcs.utils import skycoord_to_pixel\n", + "from reproject.mosaicking import find_optimal_celestial_wcs\n", + "from reproject import reproject_exact\n", "\n", "from lsst.rsp import RSPClient, get_service_url\n", "from lsst.utils.plotting import (get_multiband_plot_colors,\n", @@ -995,8 +997,8 @@ "outputs": [], "source": [ "temp_dec = snia_dec + 4.0/3600.0\n", - "temp_coord = SkyCoord(snia_ra, temp_dec, unit=\"deg\")\n", - "temp_pixels = skycoord_to_pixel(temp_coord, wcs, origin=0)" + "temp_coord_N = SkyCoord(snia_ra, temp_dec, unit=\"deg\")\n", + "temp_pixels_N = skycoord_to_pixel(temp_coord_N, wcs, origin=0)" ] }, { @@ -1035,8 +1037,8 @@ " plt.colorbar()\n", " ax[s].set_title(name)\n", " plt.plot(pixels[0], pixels[1], 'o', ms=20, mew=1, color='None', mec='yellow')\n", - " plt.plot([pixels[0], temp_pixels[0]], [pixels[1], temp_pixels[1]], color='cyan')\n", - " plt.text(temp_pixels[0], temp_pixels[1], 'N', color='cyan', fontsize=12)\n", + " plt.plot([pixels[0], temp_pixels_N[0]], [pixels[1], temp_pixels_N[1]], color='cyan')\n", + " plt.text(temp_pixels_N[0], temp_pixels_N[1], 'N', color='cyan', fontsize=12)\n", " ax[s].set_title(name)\n", "plt.suptitle('Image stamps from the SNIa alert')\n", "plt.tight_layout()\n", @@ -1094,10 +1096,235 @@ "> **Figure 3:** Image stamps from the MBA alert packet." ] }, + { + "cell_type": "markdown", + "id": "695f926f-8df0-4f16-9f2f-c7b5afef7e85", + "metadata": {}, + "source": [ + "## 2.8. WCS and reprojection of fits images" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "4faede3f-8c24-4c25-9000-0142b7cab6a5", + "metadata": {}, + "source": [ + "As noted above, when trying to create a WCS object from the alert fits image headers there are missing values which can be corrected for by defining `ctype` (a temporary solution).\n", + "\n", + "With a WCS one can plot the images in `matplotlib` using the `projection` keyword (astropy `WCSAxes`) to indicate the Right Ascension and Declination axes. In this case the images would remain aligned with the x, y pixel coordinates but these axes are not necessarily aligned with the RA and Dec axes.\n", + "Use the `find_optimal_celestial_wcs` from the `reproject` package to get a new WCS object and image shape that have RA and Dec aligned with the conventional North (up) and East (left) projection for the SNIa alert fits images.\n", + "> **Warning:** If the WCS `ctype` correction is not applied `find_optimal_celestial_wcs` will fail with TypeError: WCS does not have celestial components" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "695f926f-8df0-4f16-9f2f-c7b5afef7e85", + "id": "cf9543b7-3a91-49cb-8276-6281526e76f5", + "metadata": {}, + "outputs": [], + "source": [ + "hdul = fits.HDUList.fromstring(snia_record['cutoutScience'])\n", + "data = hdul[0].data\n", + "wcs = WCS(hdul[0].header)\n", + "wcs.wcs.ctype = [\"RA---TAN\", \"DEC--TAN\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ee8adac-c846-4ef6-8e63-5d5f12914029", + "metadata": {}, + "outputs": [], + "source": [ + "wcs_out, shape_out = find_optimal_celestial_wcs([(data, wcs)])" + ] + }, + { + "cell_type": "markdown", + "id": "64f3121f-f324-4d6e-90ef-ce5384f12f3a", + "metadata": {}, + "source": [ + "Find the pixel coordinates of the target and the North and East vectors in this new, reprojected WCS." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "89aab4a8-a90c-4589-a4fd-aa8cb3b9314c", + "metadata": {}, + "outputs": [], + "source": [ + "pixels_reproj = skycoord_to_pixel(coord, wcs_out, origin=0)\n", + "temp_pixels_N_reproj = skycoord_to_pixel(temp_coord_N, wcs_out, origin=0)\n", + "\n", + "temp_ra = snia_ra + 4.0/3600.0\n", + "temp_coord_E = SkyCoord(temp_ra, snia_dec, unit=\"deg\")\n", + "temp_pixels_E_reproj = skycoord_to_pixel(temp_coord_E, wcs_out, origin=0)" + ] + }, + { + "cell_type": "markdown", + "id": "58a834af-bac2-4b5e-b6b7-ef3e52b36826", + "metadata": {}, + "source": [ + "Use `reproject_exact` to get a new image aligned with the new WCS. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10dc8d1f-71a3-4de5-9d34-9696057f3658", + "metadata": {}, + "outputs": [], + "source": [ + "stamps = [snia_record['cutoutScience'],\n", + " snia_record['cutoutTemplate'],\n", + " snia_record['cutoutDifference']]\n", + "stamp_names = ['Science', 'Template', 'Difference']\n", + "\n", + "fig, ax = plt.subplots(1, 3, figsize=(9, 3), subplot_kw=dict(projection=wcs_out))\n", + "for s, (stamp, name) in enumerate(zip(stamps, stamp_names)):\n", + " hdul = fits.HDUList.fromstring(stamp)\n", + " data = hdul[0].data\n", + "\n", + " data_reproj, footprint = reproject_exact(input_data = (data, wcs),\n", + " output_projection = wcs_out, shape_out=shape_out,\n", + " )\n", + " \n", + " image = CCDData(data_reproj.astype(\"float32\"), unit='nJy')\n", + "# norm = ImageNormalize(image, interval=ZScaleInterval(),\n", + "# stretch=AsinhStretch())\n", + "# norm = ImageNormalize(image, interval=ZScaleInterval(),\n", + "# stretch=LinearStretch())\n", + " norm = ImageNormalize(image, interval=MinMaxInterval(),\n", + " stretch=LinearStretch())\n", + " plt.sca(ax[s])\n", + " plt.imshow(image, origin='lower', norm=norm, cmap='gray')\n", + " plt.colorbar()\n", + " ax[s].set_title(name)\n", + " plt.plot(pixels_reproj[0], pixels_reproj[1], 'o', ms=20, mew=1, color='None', mec='yellow')\n", + " plt.plot([pixels_reproj[0], temp_pixels_N_reproj[0]], [pixels_reproj[1], temp_pixels_N_reproj[1]], color='cyan')\n", + " plt.text(temp_pixels_N_reproj[0], temp_pixels_N_reproj[1], 'N', color='cyan', fontsize=12)\n", + " plt.plot([pixels_reproj[0], temp_pixels_E_reproj[0]], [pixels_reproj[1], temp_pixels_E_reproj[1]], color='magenta')\n", + " plt.text(temp_pixels_E_reproj[0], temp_pixels_E_reproj[1], 'E', color='magenta', fontsize=12)\n", + " ax[s].set_title(name)\n", + " ax[s].grid(color='white', ls='solid')\n", + "plt.suptitle('Image stamps from the SNIa alert')\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "b3aca59d-3cae-4bcf-8c33-879ac5eba450", + "metadata": {}, + "source": [ + "> **Figure 4:** Image stamps from the SNIa alert packet, reprojected such that North is up and East is left as indicated by the gridlines. The literature coordinates of the SNIa are marked with a yellow circle, and the North - East directions are indicated with cyan and magenta lines respectively." + ] + }, + { + "cell_type": "markdown", + "id": "ecf8eca7-289d-47be-a3d8-002823c84d5d", + "metadata": {}, + "source": [ + "### 2.8.1 Update the fits header with WCS terms\n", + "Inspect what has been added and/or changed by setting the WCS `ctype`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf8cf278-b97a-405c-90b4-2cde2242d26b", + "metadata": {}, + "outputs": [], + "source": [ + "hdul = fits.HDUList.fromstring(snia_record['cutoutScience'])\n", + "data = hdul[0].data\n", + "wcs = WCS(hdul[0].header)\n", + "print(\"original wcs:\\n{}\\n\".format(wcs))\n", + "wcs.wcs.ctype = [\"RA---TAN\", \"DEC--TAN\"]\n", + "print(\"corrected wcs:\\n{}\".format(wcs))" + ] + }, + { + "cell_type": "markdown", + "id": "111fbecb-b3ca-4dc2-8336-a3791a67b28d", + "metadata": {}, + "source": [ + "This can be used to update the fits header information and create a new header data unit, which could subsequently be saved to file.\n", + "> **Warning:** `output_verify='ignore'` is required to save the modified header to file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "789a0a15-a862-43d1-afba-1504ca2d80b1", + "metadata": {}, + "outputs": [], + "source": [ + "stamps = [snia_record['cutoutScience'],\n", + " snia_record['cutoutTemplate'],\n", + " snia_record['cutoutDifference']]\n", + "stamp_names = ['Science', 'Template', 'Difference']\n", + "\n", + "for s, (stamp, name) in enumerate(zip(stamps, stamp_names)):\n", + " print(\"{}:\".format(name))\n", + " hdul = fits.HDUList.fromstring(stamp)\n", + " data = hdul[0].data\n", + " hdr = hdul[0].header\n", + " wcs = WCS(hdr)\n", + " \n", + " wcs.wcs.ctype = [\"RA---TAN\", \"DEC--TAN\"]\n", + " hdr_wcs = wcs.to_header()\n", + " for x in list(hdr_wcs.keys()):\n", + " if x not in hdr:\n", + " print(\"add {} = {}\".format(x, hdr_wcs[x]))\n", + " hdr[x] = (hdr_wcs[x],hdr_wcs.comments[x])\n", + " else:\n", + " if hdr_wcs[x] != hdr[x]:\n", + " print(\"update {} = {} -> {}\".format(x, hdr[x], hdr_wcs[x]))\n", + " hdr[x] = (hdr_wcs[x],hdr_wcs.comments[x])\n", + " hdul[0] = fits.PrimaryHDU(data,header=hdr)\n", + " # hdul.writeto(\"{}_{}.fits\".format(snia_alert_id,name),output_verify='ignore')\n", + " print()" + ] + }, + { + "cell_type": "markdown", + "id": "cb70a544-9f22-4ad3-9ae2-46e292a1235a", + "metadata": {}, + "source": [ + "With the header data unit corrected, one can use the `reproject` functions on that object directly, rather than passing a tuple of data and (corrected) WCS, which may be more convenient for certain workflows." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "851fc3d2-09e8-4b66-8d15-93b38f3e828f", + "metadata": {}, + "outputs": [], + "source": [ + "# wcs_out, shape_out = find_optimal_celestial_wcs([hdul],\n", + "# hdu_in=0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "558eead5-e571-4ba4-b0ef-6b5d88054002", + "metadata": {}, + "outputs": [], + "source": [ + "# data_reproj, footprint = reproject_exact(input_data = hdul,\n", + "# output_projection = wcs_out, shape_out=shape_out,\n", + "# )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "70a1655b-55b8-4a4e-90bb-389294b1a52e", "metadata": {}, "outputs": [], "source": []