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

Open
jrob93 wants to merge 1 commit into
base: main
Choose a base branch
from
Open

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
1 commit
Select commit
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
241 changes: 234 additions & 7 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: 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)"
]
Expand Down Expand Up @@ -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",
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,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": []
Expand Down
Loading