diff --git a/tutorials/spherex/spherex_psf.md b/tutorials/spherex/spherex_psf.md index 0a78ccc9..92443c95 100644 --- a/tutorials/spherex/spherex_psf.md +++ b/tutorials/spherex/spherex_psf.md @@ -4,7 +4,7 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.17.3 + jupytext_version: 1.18.1 kernelspec: display_name: Python 3 (ipykernel) language: python @@ -24,11 +24,18 @@ authors: ## 1. Learning Goals * Determine how pixels in a SPHEREx cutout map to the pixels in the parent SPHEREx spectral image. -* Understand the structure of the PSF extension in a SPHEREx cutout (which is the same as the PSF extension in the parent spectral image) +* Understand the structure of the PSF extension in a SPHEREx cutout (which is the same as the PSF extension in the parent spectral image). +* Learn how to tell which version of the SPHEREx spectral image you are looking at, and how to interpret this information to obtain the correct PSF extension for the SPHEREx spectral images. * Learn which plane in a SPHEREx cutout PSF extension cube most accurately describes the coordinates you are interested in. +++ +```{note} +This notebook is not intended for use of QR-1 data. +``` + ++++ + ## 2. SPHEREx Overview SPHEREx is a NASA Astrophysics Medium Explorer mission that launched in March 2025. @@ -51,7 +58,7 @@ The following packages must be installed to run this notebook. ```{code-cell} ipython3 # Uncomment the next line to install dependencies if needed. -# !pip install astropy numpy pyvo +# !pip install astropy matplotlib numpy pyvo ``` ```{code-cell} ipython3 @@ -59,6 +66,8 @@ import http.client import re import time import urllib.error +import copy +from packaging.version import Version import astropy.units as u import matplotlib.pyplot as plt @@ -120,10 +129,10 @@ print("Time to do TAP query: {:2.2f} seconds.".format(time.time() - t1)) print("Number of images found: {}".format(len(results))) ``` -:::{note} +```{note} SPHEREx data are also available via SIA which can provide a simpler interface for many queries, as demonstrated in {ref}`spherex-intro`. An advantage of the method shown above is that it provides access to data immediately after ingestion (which occurs weekly) and is not subject to the same ~1 day delay as SIA. -::: +``` For this example, we focus on the first one of the retrieved SPHEREx spectral images. @@ -148,6 +157,7 @@ for attempt in range(max_retries): try: # Read the data. with fits.open(spectral_image_url) as hdul: + image_hdul = copy.deepcopy(hdul) cutout_header = hdul['IMAGE'].header psf_header = hdul['PSF'].header cutout = hdul['IMAGE'].data @@ -159,7 +169,7 @@ for attempt in range(max_retries): time.sleep(10 * (attempt + 1)) ``` -Examine the header. +Let's examine the HDU list info. ```{code-cell} ipython3 hdul.info() @@ -173,8 +183,14 @@ We have already loaded their data as well as their header. psfcube.shape ``` -The shape of the `psfcube` is (121,101,101). -This corresponds to a grid of 11x11 PSFs across the image, each of them of the size 101x101 pixels. +The shape of the `psfcube` is output above. + +```{note} +In the QR-2 data, the shape is (121,101,101), which corresponds to a grid of 11x11 PSF zones across the image. +The number of PSF zones may change in later versions of data products. +``` + +Each PSF has a size of 101x101 pixels. ```{note} Remember that the PSFs are oversampled by a factor of 10. @@ -183,18 +199,18 @@ This means that the actual size of the PSFs is about 10x10 SPHEREx pixels, which +++ -Let's look at a small part of the PSF header to understand its format: +Let's look at a small part of the PSF header to understand its format. ```{code-cell} ipython3 -psf_header[22:40] +psf_header[0:30] ``` We confirm that the oversampling factor (`OVERSAMP`) is 10. -The PSFs are distributed in an even grid with 11x11 zones. -Each of the 121 PSFs is responsible for one of these zones. +The PSFs are distributed in an even grid with NxM zones (in QR-2 data products it is N=M=11). +Each of the NxM PSFs is responsible for one of these zones. The PSF header therefore includes the center position of these zones as well as the width of the zones. -These center coordinate are specified with `XCTR_i` and `YCTR_i`, respectively, where i = 1...121. -The widths are specified with `XWID_i` and `YWID_i`, respectively, where again i = 1...121. +These center coordinate are specified with `XCTR_i` and `YCTR_i`, respectively, where i = 1...(NxM). +The widths are specified with `XWID_i` and `YWID_i`, respectively, where again i = 1...(NxM). The zones have approximately equal widths and are arranged in an even grid. The size of the zones is sufficient to capture well the changes of the PSF size and structure with wavelength and spatial coordinates. @@ -202,6 +218,303 @@ The goal of this tutorial now is to find the PSF corresponding to our input coor +++ +```{warning} +In the SPHEREx spectral image versions prior or equal to 6.5.5, there was a mismatch between the spatial layout of the PSF zones and the indexing of the PSF zones in the image header. This has now been fixed in versions 6.5.6 and beyond. + +For more information about these changes, see the following webpage: [PSF Erratum](https://irsa.ipac.caltech.edu/data/SPHEREx/docs/psfhdrerr.html) + +**Users using the old versions will need to implement an extra step to update the image header. A function to update the header is given [in Section 5.1 below](#update-psf-header-function).** +``` + +Let's first check here if a header update is necessary. We can do that by printing the `VERSION` keyword in the header. + +For comparing versions, we can use the Python-internal `Version()` function from the `packaging.version` package. Images that have already been reprocessed can have version names such as `6.5.4+psffix1` (which are superior to `6.5.4`, for example), and we can use `Version().local` to check for those. + +```{code-cell} ipython3 +this_version = Version(image_hdul['PRIMARY'].header["VERSION"]) +contains_psffix1 = this_version.local is not None and "psffix1" in this_version.local +print(f"Current version is {this_version}") + +if this_version <= Version("6.5.5") and not contains_psffix1: + print("PSF header needs to be updated! -> Go to Section 5.1 :(") +else: + print("PSF header is already up-to-date! -> Proceed to Section 6 :)") +``` + +If the version of the SPHEREx spectral image is less or equal than `6.5.5` and hasn't already been reprocessed, we will have to update the header. This is explained in Section 5.1. If the version is later than `6.5.5` or includes `"psffix1"`, the header is already updated and the PSF issue is fixed. In this case, proceed to Section 6 directly. + ++++ + +(update-psf-header-function)= +### 5.1 Updating PSF Header (for SPHEREx Spectral Image versions $\leq$ 6.5.5) + ++++ + +The function that can be used to update the header is shown below. The function +* first checks if a header update is necessary +* changes the PSF zone indexing and +* changes the version of the header such that it is consistent with the new released images + +Note that this function can work as standalone function to process many images. + +```{code-cell} ipython3 +def update_psf_header(old_hdul): + """ + Fix an old PSF FITS file header by rewriting only the per-plane header metadata + so that plane k corresponds to x-fast ordering: + k0 = iy * bins_x + ix + + The cube data are left untouched. + + Parameters + ---------- + old_hdul : fits.HDUList + Old SPHEREx Spectral Image HDUL + + Return + ---------- + new_hdul : fits.HDUList + New SPHEREx Spectral Image HDUL with updated PSF zone data in header and updated version number + """ + + VERSION_FIXED = Version("6.5.6") + PSF_FIX_TAG = "psffix1" + + def psf_fix_applied(hdul) -> bool: + """ + Return True if the PSF fix has been applied. + + Rules: + - If the VERSION header is missing in the primary HDU, the fix is not applied. + - If VERSION >= VERSION_FIXED, the fix is included in the software release. + - Otherwise the local version tag (+...) must contain PSF_FIX_TAG. + """ + header = hdul[0].header + + if "VERSION" not in header: + return False + + v = Version(header["VERSION"]) + + if v >= VERSION_FIXED: + return True + + return v.local is not None and PSF_FIX_TAG in v.local + + if psf_fix_applied(old_hdul): + return old_hdul + + ## Define some auxiliary functions ------- + def parse_ixiy_from_comment(comment): + _zone_pat = re.compile(r"\((\d+)\s*,\s*(\d+)\)") + m = _zone_pat.search(str(comment)) + if not m: + raise ValueError(f"Could not parse zone indices from comment: {comment!r}") + return int(m.group(1)), int(m.group(2)) + + def infer_grid_shape_from_header_comments(hdr, nzone): + max_ix = -1 + max_iy = -1 + + for k1 in range(1, nzone + 1): + key = f"XCTR_{k1}" + if key not in hdr: + raise KeyError(f"Missing required key: {key}") + ix, iy = parse_ixiy_from_comment(hdr.comments[key]) + max_ix = max(max_ix, ix) + max_iy = max(max_iy, iy) + + bins_x = max_ix + 1 + bins_y = max_iy + 1 + + if bins_x * bins_y != nzone: + raise ValueError( + f"Inconsistent grid inferred from comments: " + f"bins_x={bins_x}, bins_y={bins_y}, nzone={nzone}" + ) + + return bins_x, bins_y + + def collect_axis_values_by_zone(hdr, nzone): + """ + Read the old header and collect unique x/y centers and widths by zone index + labels found in the comments. + + This uses the old header only to recover the per-axis values for each ix, iy. + It does NOT use the old plane ordering as truth. + """ + x_center_by_ix = {} + y_center_by_iy = {} + x_width_by_ix = {} + y_width_by_iy = {} + + for k1 in range(1, nzone + 1): + ix, iy = parse_ixiy_from_comment(hdr.comments[f"XCTR_{k1}"]) + + xck = f"XCTR_{k1}" + yck = f"YCTR_{k1}" + xwk = f"XWID_{k1}" + ywk = f"YWID_{k1}" + + if xck in hdr: + val = hdr[xck] + if ix in x_center_by_ix and not np.isclose(x_center_by_ix[ix], val): + raise ValueError( + f"Inconsistent XCTR for ix={ix}: " + f"{x_center_by_ix[ix]} vs {val}" + ) + x_center_by_ix[ix] = val + + if yck in hdr: + val = hdr[yck] + if iy in y_center_by_iy and not np.isclose(y_center_by_iy[iy], val): + raise ValueError( + f"Inconsistent YCTR for iy={iy}: " + f"{y_center_by_iy[iy]} vs {val}" + ) + y_center_by_iy[iy] = val + + if xwk in hdr: + val = hdr[xwk] + if ix in x_width_by_ix and not np.isclose(x_width_by_ix[ix], val): + raise ValueError( + f"Inconsistent XWID for ix={ix}: " + f"{x_width_by_ix[ix]} vs {val}" + ) + x_width_by_ix[ix] = val + + if ywk in hdr: + val = hdr[ywk] + if iy in y_width_by_iy and not np.isclose(y_width_by_iy[iy], val): + raise ValueError( + f"Inconsistent YWID for iy={iy}: " + f"{y_width_by_iy[iy]} vs {val}" + ) + y_width_by_iy[iy] = val + + return x_center_by_ix, y_center_by_iy, x_width_by_ix, y_width_by_iy + ## End defining some auxiliary functions -------- + + ## Get Header + extname = "PSF" + hdu = old_hdul[extname] + cube = np.asarray(hdu.data) + hdr_in = hdu.header.copy() + + if cube.ndim != 3: + raise ValueError(f"Expected 3D PSF cube, got shape {cube.shape}") + + nzone = cube.shape[0] + bins_x, bins_y = infer_grid_shape_from_header_comments(hdr_in, nzone) + + print(f"Detected bins_x={bins_x}, bins_y={bins_y}, nzone={nzone}") + + x_center_by_ix, y_center_by_iy, x_width_by_ix, y_width_by_iy = collect_axis_values_by_zone( + hdr_in, nzone + ) + + # Validate that all needed axis values were recovered + missing = [] + for ix in range(bins_x): + if ix not in x_center_by_ix: + missing.append(f"x_center[{ix}]") + if ix not in x_width_by_ix: + missing.append(f"x_width[{ix}]") + for iy in range(bins_y): + if iy not in y_center_by_iy: + missing.append(f"y_center[{iy}]") + if iy not in y_width_by_iy: + missing.append(f"y_width[{iy}]") + + if missing: + raise ValueError(f"Missing axis metadata recovered from old header: {missing}") + + hdr_out = hdr_in.copy() + + # Rewrite only the per-plane metadata so plane k matches x-fast ordering. + # plane k0 should correspond to: + # ix = k0 % bins_x + # iy = k0 // bins_x + for k0 in range(nzone): + ix = k0 % bins_x + iy = k0 // bins_x + k1 = k0 + 1 + + hdr_out[f"XCTR_{k1}"] = ( + x_center_by_ix[ix], + f"Center of x zone ({ix}, {iy})" + ) + hdr_out[f"YCTR_{k1}"] = ( + y_center_by_iy[iy], + f"Center of y zone ({ix}, {iy})" + ) + hdr_out[f"XWID_{k1}"] = ( + x_width_by_ix[ix], + f"Width of x zone ({ix}, {iy})" + ) + hdr_out[f"YWID_{k1}"] = ( + y_width_by_iy[iy], + f"Width of y zone ({ix}, {iy})" + ) + + # Optional but useful provenance note + hdr_out["HISTORY"] = "Rewrote PSF per-plane zone metadata to x-fast ordering." + hdr_out["HISTORY"] = "Cube plane data left unchanged." + + + + new_hdu = fits.ImageHDU(data=cube, header=hdr_out, name=hdu.name) + + ext_index = old_hdul.index_of(extname) + new_hdul = fits.HDUList() + for i, old in enumerate(old_hdul): + if i == ext_index: + new_hdul.append(new_hdu) + else: + new_hdul.append(old.copy()) + + ## TO DO: UPDATE VERSION + new_hdul['PRIMARY'].header["VERSION"] = new_hdul['PRIMARY'].header["VERSION"] + "+psffix1" # SET NEW VERSION HERE + + return(new_hdul) +``` + +We now run this function to create a new HDU list that we will use later. + +```{code-cell} ipython3 +new_image_hdul = update_psf_header(old_hdul=image_hdul) +``` + +Let's check if the version keywords was updated: + +```{code-cell} ipython3 +print(f"Old version: {image_hdul['PRIMARY'].header['VERSION']}") +print(f"Updated version: {new_image_hdul['PRIMARY'].header['VERSION']}") +``` + +Let's compare the new and old PSF headers to see the difference. + +```{code-cell} ipython3 +image_hdul['PSF'].header[22:40] +``` + +```{code-cell} ipython3 +new_image_hdul['PSF'].header[22:40] +``` + +Now we have to update the variables we have set above. + +```{code-cell} ipython3 +cutout_header = new_image_hdul['IMAGE'].header +psf_header = new_image_hdul['PSF'].header +cutout = new_image_hdul['IMAGE'].data +psfcube = new_image_hdul['PSF'].data +``` + +With this fix, we are now ready to proceed! + ++++ + ## 6. Determine the Pixel Location on the Parent SPHEREx Image To identify the zone which covers the coordinates of interest, we first need to translate these coordinates to the pixel coordinates on the parent large SPHEREx image from which the cutout was created. @@ -245,7 +558,7 @@ for key, val in psf_header.items(): ym = re.match(r'(YCTR*)', key) if ym: yplane = int(key.split("_")[1]) - yctr[xplane] = val + yctr[yplane] = val ``` Check that we got all of them! @@ -317,7 +630,7 @@ To use this PSF for forward modeling or fitting, you must: ## About this notebook -**Updated:** 5 March 2026 +**Updated:** 13 March 2026 **Contact:** Contact [IRSA Helpdesk](https://irsa.ipac.caltech.edu/docs/help_desk.html) with questions or problems.