From d3b6a4827e38aa492b19a9a6ac8fee6f32637903 Mon Sep 17 00:00:00 2001 From: MelissaGraham Date: Thu, 3 Sep 2026 22:45:37 +0000 Subject: [PATCH] mlg suggested minor changes --- .../201_Alerts/201_1_Alert_packets.ipynb | 200 ++++++++---------- 1 file changed, 88 insertions(+), 112 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 fda0b407..b67478ff 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: r29.2.0\\\n", - "Last verified to run: 2026-06-17\\\n", + "LSST Science Pipelines version: r30.0.11\\\n", + "Last verified to run: 2026-09-03\\\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)" ] @@ -108,7 +108,7 @@ "source": [ "import fastavro\n", "import matplotlib.pyplot as plt\n", - "import numpy as np\n", + "# import numpy as np\n", "import pandas as pd\n", "import io\n", "from astropy.io import fits\n", @@ -890,13 +890,14 @@ " mew=2, alpha=1, color=filter_colors[filt])\n", " df_pick_band = df_lc.loc[df_lc['band'] == filt]\n", " if len(df_pick_band) > 0:\n", - " ax[0].plot(df_pick_band['midpointMjdTai']-61000, df_pick_band['psfFlux']+offsets[f], filter_symbols[filt],\n", - " ms=5, mew=1, alpha=0.7, mec=filter_colors[filt], color='None')\n", + " ax[0].plot(df_pick_band['midpointMjdTai']-61000, df_pick_band['psfFlux']+offsets[f],\n", + " filter_symbols[filt], ms=5, mew=1, alpha=0.7, mec=filter_colors[filt],\n", + " color='None')\n", " df_pick_band = df_flc.loc[df_flc['band'] == filt]\n", " if len(df_pick_band) > 0:\n", - " ax[1].plot(df_pick_band['midpointMjdTai']-61000, df_pick_band['psfFlux']+offsets[f], filter_symbols[filt],\n", - " ms=5, mew=1, alpha=0.7, mec=filter_colors[filt], color='None',\n", - " label=filt + \" +\" + str(offsets[f]))\n", + " ax[1].plot(df_pick_band['midpointMjdTai']-61000, df_pick_band['psfFlux']+offsets[f],\n", + " filter_symbols[filt], ms=5, mew=1, alpha=0.7, mec=filter_colors[filt],\n", + " color='None', label=filt + \" +\" + str(offsets[f]))\n", "ax[0].set_xlabel('MJD-61000')\n", "ax[1].set_xlabel('MJD-61000')\n", "ax[0].set_ylabel('difference flux [nJy]')\n", @@ -1096,62 +1097,71 @@ "> **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", + "#### 2.7.1. Rotate stamps to display north-up, east-left\n", + "In Figures 2 and 3 above, the image stamps are displayed at the camera's rotation angle.\n", + "This section demonstrates how to rotate the stamps to display with the convention of north-up, east-left.\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" + "Use the three image stamps for the SNIa." ] }, { "cell_type": "code", "execution_count": null, - "id": "cf9543b7-3a91-49cb-8276-6281526e76f5", + "id": "50a94caa-e013-4ffc-84c6-a647f6e97f8e", "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\"]" + "stamps = [snia_record['cutoutScience'],\n", + " snia_record['cutoutTemplate'],\n", + " snia_record['cutoutDifference']]\n", + "stamp_names = ['Science', 'Template', 'Difference']" + ] + }, + { + "cell_type": "markdown", + "id": "72061b85-78d3-4c43-9249-be5d4a9ab84e", + "metadata": {}, + "source": [ + "Use the `data` and `wcs` from the first stamp, `stamps[0]`, and 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", + "As above, correct for missing values in the WCS `ctype`.\n", + "\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": "7ee8adac-c846-4ef6-8e63-5d5f12914029", + "id": "4e2c8db9-03ed-4976-bf91-14ffa7377e0d", "metadata": {}, "outputs": [], "source": [ + "hdul = fits.HDUList.fromstring(stamps[0])\n", + "data = hdul[0].data\n", + "wcs = WCS(hdul[0].header)\n", + "wcs.wcs.ctype = [\"RA---TAN\", \"DEC--TAN\"]\n", + "\n", "wcs_out, shape_out = find_optimal_celestial_wcs([(data, wcs)])" ] }, { "cell_type": "markdown", - "id": "64f3121f-f324-4d6e-90ef-ce5384f12f3a", + "id": "8e2e2ab3-bd04-41f9-a4e5-eb1bc424391b", "metadata": {}, "source": [ - "Find the pixel coordinates of the target and the North and East vectors in this new, reprojected WCS." + "Get the pixel coordinates of the target and the North and East vectors in this new, reprojected WCS.\n", + "Then rename the " ] }, { "cell_type": "code", "execution_count": null, - "id": "89aab4a8-a90c-4589-a4fd-aa8cb3b9314c", + "id": "fa688b1e-85d0-491e-a7b0-53d3b9732a65", "metadata": {}, "outputs": [], "source": [ @@ -1165,53 +1175,53 @@ }, { "cell_type": "markdown", - "id": "58a834af-bac2-4b5e-b6b7-ef3e52b36826", + "id": "9de12e04-70f0-4d1b-98db-c9bf61c33991", "metadata": {}, "source": [ - "Use `reproject_exact` to get a new image aligned with the new WCS. " + "Recreate Figure 2, but with the images displayed in the new north-up, east-left projection\n", + "* Omit the colorbar because it is large; uncomment the line to display the colorbar again.\n", + "* Use of shorter variable names for the reprojected pixels is only to enable clearer, shorter plotting commands." ] }, { "cell_type": "code", "execution_count": null, - "id": "10dc8d1f-71a3-4de5-9d34-9696057f3658", + "id": "c9707382-b089-44fc-b90b-25340d6b6ef1", "metadata": {}, "outputs": [], "source": [ - "stamps = [snia_record['cutoutScience'],\n", - " snia_record['cutoutTemplate'],\n", - " snia_record['cutoutDifference']]\n", - "stamp_names = ['Science', 'Template', 'Difference']\n", + "pix = pixels_reproj\n", + "pix_N = temp_pixels_N_reproj\n", + "pix_E = temp_pixels_E_reproj\n", "\n", - "fig, ax = plt.subplots(1, 3, figsize=(9, 3), subplot_kw=dict(projection=wcs_out))\n", + "fig, ax = plt.subplots(1, 3, figsize=(12, 4), 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", + " data_reproj, footprint = reproject_exact(input_data=(data, wcs), output_projection=wcs_out,\n", + " shape_out=shape_out)\n", + " image = CCDData(data_reproj.astype(\"float32\"), unit=\"nJy\")\n", + " norm = ImageNormalize(image, interval=MinMaxInterval(), 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", + " # plt.colorbar()\n", + "\n", + " plt.plot(pix[0], pix[1], 'o', ms=20, mew=1, color='None', mec='yellow')\n", + " plt.plot([pix[0], pix_N[0]], [pix[1], pix_N[1]], color='cyan')\n", + " plt.text(pix_N[0]-1, pix_N[1]+1, 'N', color='cyan', fontsize=12)\n", + " plt.plot([pix[0], pix_E[0]], [pix[1], pix_E[1]], color='magenta')\n", + " plt.text(pix_E[0]-1, pix_E[1]+1, 'E', color='magenta', fontsize=12)\n", + "\n", " ax[s].grid(color='white', ls='solid')\n", + " ax[s].set_xlabel('RA', fontsize=10)\n", + " ax[s].tick_params(axis='x', labelsize=10)\n", + " ax[s].set_ylabel('Dec', fontsize=10)\n", + " ax[s].tick_params(axis='y', labelsize=10)\n", + " if s > 0:\n", + " ax[s].tick_params(axis='y', labelleft=False)\n", + " ax[s].set_title(name)\n", + "\n", "plt.suptitle('Image stamps from the SNIa alert')\n", - "plt.tight_layout()\n", "plt.show()" ] }, @@ -1228,31 +1238,15 @@ "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", + "#### 2.7.2. Update the FITS header\n", + "\n", + "Above, the WCS `ctype` is corrected in order to display the image stamp.\n", + "The following demonstrates how to propagate that change to the WCS into the FITS header `hdul`.\n", + "Every change made to the header of each image stamp is printed as output of the next cell.\n", + "\n", + "A FITS file with the updated header can be saved to file by uncommenting the \"`hdule.writeto()`\" line.\n", + "The filenames will be, e.g., \"170301167371288762_Difference.fits\".\n", + "\n", "> **Warning:** `output_verify='ignore'` is required to save the modified header to file." ] }, @@ -1274,19 +1268,19 @@ " data = hdul[0].data\n", " hdr = hdul[0].header\n", " wcs = WCS(hdr)\n", - " \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", + " 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", + " 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()" ] }, @@ -1295,30 +1289,12 @@ "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", - "# )" + "With the header data unit corrected, it is possible to then use the `reproject` functions on the `hdul` directly, with the example code below, which may be more convenient for certain workflows.\n", + "\n", + "```\n", + "wcs_out, shape_out = find_optimal_celestial_wcs([hdul], hdu_in=0)\n", + "data_reproj, footprint = reproject_exact(input_data=hdul, output_projection=wcs_out, shape_out=shape_out)\n", + "```" ] }, {