Issue 276 (PSF header update)

[afaisst]

PR for issue #276.

Changes:

• Added warning that PSF header needs to be updated for data versions <=6.5.5. Added link to PSF Erratum (webpage still needs to be added)
• Added learning goal: learn which version of the data you are looking at, and how to use this version information to correctly interpret the PSF extension of the SPHEREx spectral images.
• Added check for SPHEREx image version to notify user if the PSF header needs to be updated.
• Added function to update PSF header to fix incorrectly described zone centers. Tell use that PSF needs to be updated if version <=6.5.5. The function corrects the PSF header and updates version in primary header.
• Added note that number of PSF zones can change in later pipeline versions

Todo:

• move the function that updates the PSF header (pretty long, might want to put it somewhere else and then import) to another python file within the repo

[vandesai1]

The documentation page is going to be at:

https://irsa.ipac.caltech.edu/data/SPHEREx/docs/psfhdrerr.html

[afaisst]

The documentation page is going to be at:

https://irsa.ipac.caltech.edu/data/SPHEREx/docs/psfhdrerr.html

Added!

[gpdf]

For future reference, could you update the PR description above with the “<= 6.5.5” correction?

[afaisst]

Made following changes:

• removed warning in beginning
• mentioned that not intended for QR-1 data
• added link to function to update header

[tgoldina]

I checked that the updated hdul is correct.

For IRSA, I am using the slightly different version, which produces the same result:

    def build_corrected_psf_header(hdul) -> fits.Header:
        psf_header = hdul["PSF"].header.copy()
        # Count zones
        zone_id = 1
        while f'XCTR_{zone_id}' in psf_header:
            zone_id += 1
        n_zones = zone_id - 1

        if n_zones == 0:
            raise ValueError("No XCTR_* keywords found in PSF header")

        # Get unique centers to determine grid dimensions
        x_ctr_vals = [psf_header[f'XCTR_{i}'] for i in range(1, n_zones + 1)]
        y_ctr_vals = [psf_header[f'YCTR_{i}'] for i in range(1, n_zones + 1)]

        x_centers = np.array(sorted(set(x_ctr_vals)))
        y_centers = np.array(sorted(set(y_ctr_vals)))

        x_zones = len(x_centers)
        y_zones = len(y_centers)

        if x_zones * y_zones != n_zones:
            raise ValueError(f"Zone count mismatch: {x_zones}*{y_zones} != {n_zones}")

        # Extract widths using Y-fast pattern (how they were incorrectly written)
        x_widths = np.array([psf_header[f'XWID_{ii * y_zones + 1}'] for ii in range(x_zones)])
        y_widths = np.array([psf_header[f'YWID_{jj + 1}'] for jj in range(y_zones)])

        # Now fix the header to X-fast ordering
        zone_id = 1
        for iy in range(y_zones):
            for ix in range(x_zones):
                location = f"({ix}, {iy})"
                psf_header[f'XCTR_{zone_id}'] = (x_centers[ix], f"Center of x zone {location}")
                psf_header[f'YCTR_{zone_id}'] = (y_centers[iy], f"Center of y zone {location}")
                psf_header[f'XWID_{zone_id}'] = (x_widths[ix], f"Width of x zone {location}")
                psf_header[f'YWID_{zone_id}'] = (y_widths[iy], f"Width of y zone {location}")
                zone_id += 1
        return psf_header

[tgoldina]

Everything before + is the public version. For example, Version("6.5.4+psffix1").publicreturns6.5.4`. If comparisons are performed using the public version, no wrapper function is needed.

[troyraen]

Andreas is out for a couple days and I'm working on wrapping this up. After looking at your comment here and also the one below, I removed the wrapper function and replaced this version comparison with:

this_version = Version(image_hdul['PRIMARY'].header["VERSION"])
contains_psffix1 = this_version.local is not None and "psffix1" in this_version.local
if this_version <= Version("6.5.5") and not contains_psffix1:
    print("PSF header needs to be updated! -> Go to Section 5.1 :(")

I also updated the surrounding text. You can see all the changes I made based on your feedback here: 0f249d0. Please let me know if anything seems off.

[tgoldina]

This can be replaced with:

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

[troyraen]

Done

[christinanelson]

in section 7, I think there is a bug:

xctr = {}
yctr = {}

for key, val in psf_header.items():
    # Look for keys like XCTR* or YCTR*
    xm = re.match(r'(XCTR*)', key)
    if xm:
        xplane = int(key.split("_")[1])
        xctr[xplane] = val
    ym = re.match(r'(YCTR*)', key)
    if ym:
        yplane = int(key.split("_")[1])
        yctr[xplane] = val

the last line should say:
yctr[yplane] = val

[troyraen]

Fixed

[christinanelson]

The last line, should probably be?
elif this_version > parse_version("6.5.5"):

[troyraen]

I replaced this with the suggestion from above: #288 (comment)

[christinanelson]

This notebook, the EstimatePSFTask, and the PSFZoneData class are all consistent; they assume x-fast upon write, when reading interprets under the assumption that the data is stored x-fast and will allow reading of a header x-fast or y-fast. The class PSFZoneData class now has a function that allows a header fix, which is the same as in this notebook.

---

• Here is another confusion, the comments are referring to PSF zones, but the code snippet seems to be calculating pixel distance?

# Note that the x,y coordinates of the PSF zone centers are in 1,1 convention,
# therefore we have to subtract 1 pixels.
tab["distance"] = np.sqrt((tab["x"]-1 - xpix_orig)**2 + (tab["y"]-1 - ypix_orig)**2)

• Another bug is that the PSFZones are 0-based but should be 1-based. We also need to make this consistent. @tgoldina should we do the 1-based zone fix at the same time when fixing the headers?

[troyraen]

Thanks for your reviews @tgoldina and @christinanelson. Andreas is out for a couple days and I'm working on wrapping this up. I've implemented all of your feedback (see 0f249d0 and f551b24) except the two below. Let me know if there's anything else you'd like to change/add. I plan to merge this soon after the erratum page gets published at https://irsa.ipac.caltech.edu/data/SPHEREx/docs/psfhdrerr.html, which I expect to happen sometime today. We can always do a quick followup PR if change requests come in after it merges.

---

Another bug is that the PSFZones are 0-based but should be 1-based. We also need to make this consistent. @tgoldina should we do the 1-based zone fix at the same time when fixing the headers?

I was sure we had an open issue for this that Andreas was planning to complete soon, but it looks like I never opened the issue. If you want to tell me what to change I can add it to this PR. Otherwise I'll open the issue and ask Andreas to take care of it when he gets back.

---

Here is another confusion, the comments are referring to PSF zones, but the code snippet seems to be calculating pixel distance?

# Note that the x,y coordinates of the PSF zone centers are in 1,1 convention,
# therefore we have to subtract 1 pixels.
tab["distance"] = np.sqrt((tab["x"]-1 - xpix_orig)**2 + (tab["y"]-1 - ypix_orig)**2)

I wonder if you were looking at an older version of the code. The current version is the following. It makes sense to me but I'm relatively naive here so let me know if you still think something should be changed.

Once we have created this dictionary with zone pixel coordinates, we can simply search for the closest zone center to the coordinates of interest. For this we first add the distance between zone center coordinates and coordinates of interest to the table. (Note that the x,y coordinates of the PSF zone centers are in 1,1 convention, therefore we have to subtract 1 pixels.)

tab["distance"] = np.sqrt((tab["x"]-1 - xpix_orig)**2 + (tab["y"]-1 - ypix_orig)**2)

[troyraen]

Thanks @afaisst! I made the changes requested above and did some minor cleanup. Good to go.