Skip to content

HSM: measure PSF/star shapes directly in sky coords (use_sky_coords) #812

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
cailmdaley wants to merge 5 commits into
base: develop
Choose a base branch
from
Open

HSM: measure PSF/star shapes directly in sky coords (use_sky_coords) #812

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
e8f5ff7
psf shapes: measure HSM adaptive moments in sky coords
cailmdaley Jul 12, 2026
71d502f
convert_psf_pix2world: strip now-redundant HSM shape rotation
cailmdaley Jul 12, 2026
071a604
test: HSM sky-coords moments match old Jacobian rotation
cailmdaley Jul 12, 2026
c08cdb0
convert_psf_pix2world: keep HSM rotation for external-mccd MOM_LIST b…
cailmdaley Jul 12, 2026
f153ea2
psf interp: use origin=1 positions for local_wcs_list
cailmdaley Jul 12, 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
12 changes: 9 additions & 3 deletions docs/source/pipeline_canfar.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -333,10 +333,16 @@ shapepipe_run -c $SP_CONFIG/config_Ms_$psf.ini
shapepipe_run -c $SP_CONFIG/config_Pl_$psf.ini
```

#### Convert star catalogue to wCS
#### Collate star catalogues per patch

Convert all input validation PSF files and create directories per patch `P?`.
Create files `validation_psf_conv-<patchnum>-<idx>.fits` (for the v1.4 setup only one file):
Collate all input validation PSF files and create directories per patch `P?`,
producing files `validation_psf_conv-<patchnum>-<idx>.fits` (for the v1.4 setup
only one file). This gathers positions (X/Y/RA/DEC) and the MCCD CCD id.

Note: HSM shapes are no longer rotated into world coordinates at this step. The
PSF/star ellipticities and sizes are now measured directly in sky coordinates
during PSF interpolation (galsim `FindAdaptiveMom(use_sky_coords=True)`), so
`convert_psf_pix2world.py` only collates and passes the shapes through.

@martinkilbinger martinkilbinger Jul 14, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Update to reflect removal of the patch concept in v2.0

```bash
cd /path/to/version
Expand Down
104 changes: 44 additions & 60 deletions scripts/python/convert_psf_pix2world.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,25 @@
#! /usr/bin/env python3

"""CONVERT PSF CATALOGUE.

Collate the per-exposure PSF validation catalogues into per-patch star
catalogues: gather positions (X/Y/RA/DEC), assign the MCCD focal-plane CCD id,
and merge into ``validation_psf_conv-<patch>-<idx>.fits``.

HSM ellipticities and sizes are no longer rotated here. PSFEx and the in-repo
MCCD interpolation now measure adaptive moments directly in world coordinates
(galsim ``FindAdaptiveMom(use_sky_coords=True)``), so the WCS-Jacobian shape
rotation this script used to perform is redundant and has been removed.

Caveat: the MCCD ``PSF_MOM_LIST``/``STAR_MOM_LIST`` columns are produced by the
external ``mccd`` fit-validation code (``mccd.auxiliary_fun.mccd_validation``),
which still measures HSM moments in the pixel frame. Those shapes are therefore
still rotated into world coordinates here, via the WCS-Jacobian rotation, until
``mccd`` itself adopts ``use_sky_coords``. This is the one branch that keeps the
rotation; the in-repo PSFEx and MCCD-interpolation paths measure adaptive
moments directly in world coordinates upstream and pass them straight through.
"""

import sys
import os
import re
Expand Down Expand Up @@ -555,7 +575,6 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
if self._params["psf"] == "psfex":
psf_file_hdus = fits.open(psf_file_path, memmap=False)
psf_file = psf_file_hdus[2].data
header_file = psf_file_hdus[1].data
psf_file_hdus.close()
mod = "RA"
else:
Expand All @@ -564,60 +583,10 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
except Exception:
continue

new_e1_psf = np.zeros_like(psf_file[mod])
new_e2_psf = np.zeros_like(psf_file[mod])
new_sig_psf = np.zeros_like(psf_file[mod])
new_e1_star = np.zeros_like(psf_file[mod])
new_e2_star = np.zeros_like(psf_file[mod])
new_sig_star = np.zeros_like(psf_file[mod])

if self._params["psf"] == "psfex":
header = fits.Header.fromstring(
"\n".join(header_file[0][0]), sep="\n"
)
wcs = galsim.AstropyWCS(header=header)

k = 0
for ind, obj in enumerate(psf_file):
try:
# jac = wcs.jacobian(world_pos=galsim.CelestialCoord(
# ra=obj["RA"]*galsim.degrees,
# dec=obj["DEC"]*galsim.degrees
# ))
jac = wcs.jacobian(
image_pos=galsim.PositionD(
obj["X"],
obj["Y"],
)
)
except Exception:
continue
g1_psf_tmp, g2_psf_tmp, sig_psf_tmp = transform_shape(
[
obj["E1_PSF_HSM"],
obj["E2_PSF_HSM"],
obj["SIGMA_PSF_HSM"],
],
jac,
)

new_e1_psf[ind] = g1_psf_tmp
new_e2_psf[ind] = g2_psf_tmp
new_sig_psf[ind] = sig_psf_tmp

g1_star_tmp, g2_star_tmp, sig_star_tmp = transform_shape(
[
obj["E1_STAR_HSM"],
obj["E2_STAR_HSM"],
obj["SIGMA_STAR_HSM"],
],
jac,
)
new_e1_star[ind] = g1_star_tmp
new_e2_star[ind] = g2_star_tmp
new_sig_star[ind] = sig_star_tmp
k += 1

# HSM ellipticities and sizes are measured directly in world
# coordinates upstream (FindAdaptiveMom use_sky_coords=True), so
# they are passed straight through; only positions are collated.
exp_cat = np.array(
list(
map(
Expand All @@ -628,13 +597,13 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
psf_file["Y"],
psf_file["RA"],
psf_file["DEC"],
new_e1_psf,
new_e2_psf,
new_sig_psf,
psf_file["E1_PSF_HSM"],
psf_file["E2_PSF_HSM"],
psf_file["SIGMA_PSF_HSM"],
psf_file["FLAG_PSF_HSM"],
new_e1_star,
new_e2_star,
new_sig_star,
psf_file["E1_STAR_HSM"],
psf_file["E2_STAR_HSM"],
psf_file["SIGMA_STAR_HSM"],
psf_file["FLAG_STAR_HSM"],
np.ones_like(psf_file["RA"], dtype=int)
* ccd_id,
Expand All @@ -661,8 +630,23 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
]
)

# Local-to-CCD position: subtract each CCD's focal-plane shift.
new_x = np.zeros_like(psf_file[mod])
new_y = np.zeros_like(psf_file[mod])

# The MCCD PSF_MOM_LIST/STAR_MOM_LIST columns come from the
# external mccd fit-validation code, which still measures HSM
# moments in the pixel frame; rotate them into world coordinates
# via the per-CCD WCS Jacobian. This rotation stays until mccd
# itself adopts use_sky_coords (see the module docstring). The
# in-repo PSFEx / MCCD-interpolation paths are already in world
# coordinates and are passed through unrotated.
new_e1_psf = np.zeros_like(psf_file[mod])
new_e2_psf = np.zeros_like(psf_file[mod])
new_sig_psf = np.zeros_like(psf_file[mod])
new_e1_star = np.zeros_like(psf_file[mod])
new_e2_star = np.zeros_like(psf_file[mod])
new_sig_star = np.zeros_like(psf_file[mod])
new_flag_psf = np.zeros_like(psf_file[mod])
new_flag_star = np.zeros_like(psf_file[mod])
for ccd_id in range(40):
Expand Down
67 changes: 57 additions & 10 deletions src/shapepipe/modules/mccd_package/mccd_interpolation_script.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
from sqlitedict import SqliteDict

from shapepipe.pipeline import file_io
from shapepipe.modules.psfex_interp_package.psfex_interp import local_wcs_list

try:
import galsim.hsm as hsm
Expand Down Expand Up @@ -204,18 +205,34 @@ def _get_galaxy_positions(self):
raise KeyError(pos_param_err)
galcat.close()

def _get_psfshapes(self):
def _get_psfshapes(self, wcs_list=None):
"""Get PSF Shapes.

Compute shapes of PSF at galaxy positions using HSM.

Parameters
----------
wcs_list : list of galsim.BaseWCS, optional
Local WCS at each object position. When provided, HSM adaptive
moments are measured in world coordinates (``use_sky_coords=True``),
so the returned shapes and sizes are already in sky coordinates.

"""
if import_fail:
raise ImportError("Galsim is required to get shapes information")

psf_moms = [
hsm.FindAdaptiveMom(Image(psf), strict=False) for psf in self.interp_PSFs
]
if wcs_list is not None:
psf_moms = [
hsm.FindAdaptiveMom(
Image(psf, wcs=wcs), strict=False, use_sky_coords=True
)
for psf, wcs in zip(self.interp_PSFs, wcs_list)
]
else:
psf_moms = [
hsm.FindAdaptiveMom(Image(psf), strict=False)
for psf in self.interp_PSFs
]

self.psf_shapes = np.array(
[
Expand Down Expand Up @@ -253,7 +270,7 @@ def _write_output(self):
data = {"VIGNET": self.interp_PSFs}
output.save_as_fits(data, sex_cat_path=self._galcat_path)

def _get_starshapes(self, star_vign):
def _get_starshapes(self, star_vign, wcs_list=None):
"""Get Star Shapes.

Compute shapes of stars at stars positions using HSM.
Expand All @@ -262,6 +279,10 @@ def _get_starshapes(self, star_vign):
----------
star_vign : numpy.ndarray
Array containing the star's vignets
wcs_list : list of galsim.BaseWCS, optional
Local WCS at each star position; see :meth:`_get_psfshapes`. When
provided, moments are measured in world coordinates
(``use_sky_coords=True``).

"""
if import_fail:
Expand All @@ -270,10 +291,23 @@ def _get_starshapes(self, star_vign):
masks = np.zeros_like(star_vign)
masks[np.where(star_vign == -1e30)] = 1

star_moms = [
hsm.FindAdaptiveMom(Image(star), badpix=Image(mask), strict=False)
for star, mask in zip(star_vign, masks)
]
if wcs_list is not None:
star_moms = [
hsm.FindAdaptiveMom(
Image(star, wcs=wcs),
badpix=Image(mask),
strict=False,
use_sky_coords=True,
)
for star, mask, wcs in zip(star_vign, masks, wcs_list)
]
else:
star_moms = [
hsm.FindAdaptiveMom(
Image(star), badpix=Image(mask), strict=False
)
for star, mask in zip(star_vign, masks)
]

self.star_shapes = np.array(
[
Expand Down Expand Up @@ -451,6 +485,15 @@ def _interpolate_me(self):
0,
)
).T
# 1-indexed (FITS) positions for local_wcs_list, which wraps
# galsim.AstropyWCS; gal_pos (origin 0) stays for MCCD interp.
gal_pos_wcs = np.array(
self._f_wcs_file[exp_name][ccd]["WCS"].all_world2pix(
self.gal_pos[:, 0][ind_obj],
self.gal_pos[:, 1][ind_obj],
1,
)
).T

self.interp_PSFs = interp_MCCD(mccd_model_path, gal_pos, ccd)

Expand Down Expand Up @@ -497,7 +540,11 @@ def _interpolate_me(self):
array_id = np.concatenate((array_id, np.copy(obj_id)))

if self._compute_shape:
self._get_psfshapes()
self._get_psfshapes(
local_wcs_list(
self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos_wcs
)
)
if array_shape is None:
array_shape = np.copy(self.psf_shapes)
else:
Expand Down
18 changes: 16 additions & 2 deletions src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
from astropy.io import fits
from cs_util import size as cs_size

from shapepipe.modules.psfex_interp_package.psfex_interp import local_wcs_list
from shapepipe.pipeline import file_io

NOT_ENOUGH_STARS = "Not enough stars to train the model."
Expand Down Expand Up @@ -403,6 +404,10 @@ def mccd_interpolation_pipeline(
x_pos = galcat[2].data[pos_params[0]]
y_pos = galcat[2].data[pos_params[1]]
interp_pos = np.array([x_pos, y_pos]).T
# CCD image WCS, stored as a card list in the FITS_LDAC LDAC_IMHEAD HDU
wcs_header = fits.Header.fromstring(
"\n".join(galcat[1].data[0][0]), sep="\n"
)
# Close catalog
galcat.close()

Expand All @@ -411,9 +416,18 @@ def mccd_interpolation_pipeline(

if interp_PSFs is not None:
if get_shapes:
# Local WCS at each object position lets HSM measure adaptive
# moments directly in world coordinates (use_sky_coords=True).
wcs_list = local_wcs_list(
galsim.AstropyWCS(header=wcs_header), interp_pos
)
PSF_moms = [
galsim.hsm.FindAdaptiveMom(galsim.Image(psf), strict=False)
for psf in interp_PSFs
galsim.hsm.FindAdaptiveMom(
galsim.Image(psf, wcs=wcs),
strict=False,
use_sky_coords=True,
)
for psf, wcs in zip(interp_PSFs, wcs_list)
]

PSF_shapes = np.array(
Expand Down
Loading