Skip to content

SP-3404: minor updates to prompt/201-1 as part of review #209

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
MelissaGraham wants to merge 1 commit into
base: tickets/SP-3383
Choose a base branch
from
Open

SP-3404: minor updates to prompt/201-1 as part of review #209

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
200 changes: 88 additions & 112 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: 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)"
]
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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": [
Expand All @@ -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()"
]
},
Expand All @@ -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."
]
},
Expand All @@ -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()"
]
},
Expand All @@ -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",
"```"
]
},
{
Expand Down
Loading