Skip to content

SP-3383: update an alert tutorial to demo how to rotate stamps to N-up E-left #205

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jrob93 merged 3 commits into from
Sep 17, 2026
Merged

SP-3383: update an alert tutorial to demo how to rotate stamps to N-up E-left #205

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3e9bec0
add fits reprojection
jrob93 Aug 28, 2026
952a86c
mlg suggested minor changes
MelissaGraham Sep 3, 2026
28e6ea1
minor text updates
jrob93 Sep 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
229 changes: 215 additions & 14 deletions Prompt/200_Data_products/201_Alerts/201_1_Alert_packets.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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: 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)"
]
Expand Down Expand Up @@ -96,7 +96,7 @@
"\n",
"Import the `fastavro` package for reading alerts from the alert retrieval service in the default Avro format.\n",
"Import the `RSPClient` and `get_service_url` in order to access the alert retrieval service.\n",
"Also import a range of other plotting and analyis packages."
"Also import a range of other plotting and analysis packages."
]
},
{
Expand All @@ -108,7 +108,6 @@
"source": [
"import fastavro\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import io\n",
"from astropy.io import fits\n",
Expand All @@ -119,6 +118,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",
Expand Down Expand Up @@ -888,13 +889,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",
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -1094,10 +1096,209 @@
"> **Figure 3:** Image stamps from the MBA alert packet."
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "4faede3f-8c24-4c25-9000-0142b7cab6a5",
"metadata": {},
"source": [
"#### 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",
"Use the three image stamps for the SNIa."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "50a94caa-e013-4ffc-84c6-a647f6e97f8e",
"metadata": {},
"outputs": [],
"source": [
"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": "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": "8e2e2ab3-bd04-41f9-a4e5-eb1bc424391b",
"metadata": {},
"source": [
"Get the pixel coordinates of the target and the North and East vectors in this new, reprojected WCS."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fa688b1e-85d0-491e-a7b0-53d3b9732a65",
"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": "9de12e04-70f0-4d1b-98db-c9bf61c33991",
"metadata": {},
"source": [
"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": "c9707382-b089-44fc-b90b-25340d6b6ef1",
"metadata": {},
"outputs": [],
"source": [
"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=(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",
" 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",
"\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.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.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."
]
},
{
"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, 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",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "695f926f-8df0-4f16-9f2f-c7b5afef7e85",
"id": "c8c988fb-ba48-4bc9-8a61-253f93da06b9",
"metadata": {},
"outputs": [],
"source": []
Expand Down
Loading