From e8f5ff744557fdc43c70514ba483dcf1e3609b34 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 12 Jul 2026 23:34:13 +0200 Subject: [PATCH 1/5] psf shapes: measure HSM adaptive moments in sky coords Attach the local (Jacobian) WCS at each object position to the PSF/star stamps and pass use_sky_coords=True to galsim FindAdaptiveMom, so PSFEx and MCCD HSM ellipticities and sizes come out already in world coordinates. This makes the downstream WCS-Jacobian shape rotation (convert_psf_pix2world) redundant. Threaded at every in-repo HSM PSF/star call site: - psfex_interp: multi-epoch (per-CCD WCS log), single-exposure validation and classic (WCS from the galcat FITS_LDAC image header via _galcat_wcs) - mccd_interpolation_script: multi-epoch (per-CCD WCS log) - shapepipe_auxiliary_mccd: interpolation (WCS from galcat LDAC header) local_wcs_list() builds the per-object local WCS: use_sky_coords evaluates the stamp WCS at true_center, so it must be the full WCS linearised at the object's image position, not at the stamp centre. SIGMA_*_HSM (hence HSM_T_* = 2 sigma^2) is now in sky units (arcsec) rather than pixels; this is the intended, physically correct size. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w --- .../mccd_package/mccd_interpolation_script.py | 58 ++++++-- .../mccd_package/shapepipe_auxiliary_mccd.py | 18 ++- .../psfex_interp_package/psfex_interp.py | 126 ++++++++++++++++-- 3 files changed, 176 insertions(+), 26 deletions(-) diff --git a/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py b/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py index bfd9f3417..468cf55ee 100644 --- a/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py +++ b/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py @@ -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 @@ -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( [ @@ -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. @@ -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: @@ -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( [ @@ -497,7 +531,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 + ) + ) if array_shape is None: array_shape = np.copy(self.psf_shapes) else: diff --git a/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py b/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py index 943433e8f..f0aea4f65 100644 --- a/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py +++ b/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py @@ -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." @@ -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() @@ -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( diff --git a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py index c14aeb136..96b2c2007 100644 --- a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py +++ b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py @@ -17,6 +17,7 @@ from shapepipe.pipeline import file_io try: + import galsim import galsim.hsm as hsm from galsim import Image except ImportError: @@ -30,6 +31,37 @@ FILE_NOT_FOUND = "File_not_found" +def local_wcs_list(wcs, positions): + """Local WCS List. + + Build a list of local (Jacobian) WCS from an astropy/galsim WCS, one per + object position. Attaching the local WCS to an object-centred stamp lets + ``FindAdaptiveMom(use_sky_coords=True)`` measure moments directly in world + coordinates: at the stamp's ``true_center`` galsim evaluates this local WCS, + which must therefore be the linearisation of the full WCS *at the object's + position in the original image*, not at the stamp centre. + + Parameters + ---------- + wcs : galsim.BaseWCS or astropy.wcs.WCS + Full (celestial) WCS for the CCD. + positions : numpy.ndarray + Object positions ``(N, 2)`` in image (pixel) coordinates. + + Returns + ------- + list of galsim.JacobianWCS + One local WCS per object position. + + """ + if not isinstance(wcs, galsim.BaseWCS): + wcs = galsim.AstropyWCS(wcs=wcs) + return [ + wcs.local(image_pos=galsim.PositionD(float(x), float(y))) + for x, y in positions + ] + + class PSFExInterpolator(object): """The PSFEx Interpolator Class. @@ -145,7 +177,9 @@ def process(self): self._w_log.info(f"Psf model file {self._dotpsf_path} not found.") else: if self._compute_shape: - self._get_psfshapes() + self._get_psfshapes( + local_wcs_list(self._galcat_wcs(), self.gal_pos) + ) self._write_output() @@ -307,18 +341,36 @@ def _interpolate(self): self.gal_pos, ) - 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 (see :func:`local_wcs_list`). When + provided, HSM adaptive moments are measured directly in world + coordinates (``use_sky_coords=True``), so the ``g1``/``g2``/``sigma`` + returned are already in sky coordinates and need no post-hoc + rotation. When ``None``, moments are measured in the pixel frame. + """ 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( [ @@ -407,21 +459,30 @@ def process_validation(self, psfex_cat_path): star_dict["SNR"] = np.copy(star_cat.get_data()["SNR_WIN"]) star_cat.close() - self._get_psfshapes() - self._get_starshapes(star_vign) + wcs_list = local_wcs_list( + self._galcat_wcs(), + np.column_stack((star_dict["X"], star_dict["Y"])), + ) + + self._get_psfshapes(wcs_list) + self._get_starshapes(star_vign, wcs_list) psfex_cat_dict = self._get_psfexcatdict(psfex_cat_path) self._write_output_validation(star_dict, psfex_cat_dict) - 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. Parameters ---------- - numpy.ndarray + 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: @@ -430,10 +491,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( [ @@ -447,6 +521,26 @@ def _get_starshapes(self, star_vign): ] ) + def _galcat_wcs(self): + """Galaxy Catalogue WCS. + + Read the CCD image WCS from the SExtractor FITS_LDAC galaxy catalogue. + The original image header is stored as a card list in the ``LDAC_IMHEAD`` + HDU (HDU 1); this is the same header the (now retired) pixel-to-world + conversion reconstructed to rotate shapes. + + Returns + ------- + galsim.AstropyWCS + WCS of the CCD the stars/PSFs live on. + + """ + with fits.open(self._galcat_path, memmap=False) as hdu_list: + header = fits.Header.fromstring( + "\n".join(hdu_list[1].data[0][0]), sep="\n" + ) + return galsim.AstropyWCS(header=header) + def _get_psfexcatdict(self, psfex_cat_path): """Get PSFEx Catalogue Dictionary. @@ -700,7 +794,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 + ) + ) if array_shape is None: array_shape = np.copy(self.psf_shapes) else: From 71d502f7b1fae443f90b4bf43252afcce612c010 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 12 Jul 2026 23:37:02 +0200 Subject: [PATCH 2/5] convert_psf_pix2world: strip now-redundant HSM shape rotation HSM ellipticities and sizes are measured directly in world coordinates upstream, so the WCS-Jacobian shape rotation (transform_shape, getDecomposition) is redundant. Remove it from both the PSFEx and MCCD branches; the script keeps its remaining job of collating positions (X/Y/ RA/DEC), assigning the MCCD focal-plane CCD id, and merging per-exposure catalogues into per-patch star catalogues. Drops the galsim dependency. Caveat documented in the module docstring: the MCCD PSF_MOM_LIST/STAR_MOM_LIST columns come from the external mccd fit-validation code, still measured in the pixel frame; they are now passed through unrotated. Update the CANFAR doc to describe the step as collation, not pix->world shape conversion. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w --- docs/source/pipeline_canfar.md | 12 +- scripts/python/convert_psf_pix2world.py | 206 +++++------------------- 2 files changed, 45 insertions(+), 173 deletions(-) diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index 7e8bb29c5..a3ec03e7b 100644 --- a/docs/source/pipeline_canfar.md +++ b/docs/source/pipeline_canfar.md @@ -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--.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--.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. ```bash cd /path/to/version diff --git a/scripts/python/convert_psf_pix2world.py b/scripts/python/convert_psf_pix2world.py index dea99656b..72faa34c8 100755 --- a/scripts/python/convert_psf_pix2world.py +++ b/scripts/python/convert_psf_pix2world.py @@ -1,5 +1,23 @@ #! /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--.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, which still measures in the pixel frame; +those shapes are passed through unrotated here. If the MCCD validation branch is +revived, its moments must be measured in sky coordinates upstream (as the PSFEx +and MCCD interpolation paths now are). +""" + import sys import os import re @@ -10,46 +28,11 @@ import numpy as np from astropy.io import fits -import galsim from cs_util import args as cs_args from cs_util import logging -def transform_shape(mom_list, jac): - """Transform Shape. - - Transform shape (ellipticity and size) using a Jacobian. - - Parameters - ---------- - mom_list : list - input moment measurements; each list element contains - first and second ellipticity component and size - jac : galsim.JacobianWCS - Jacobian transformation matrix information - - Returns - ------- - list - transformed shape parameters, which are - first and second ellipticity component and size - - """ - scale, shear, theta, flip = jac.getDecomposition() - - sig_tmp = mom_list[2] * scale - shape = galsim.Shear(g1=mom_list[0], g2=mom_list[1]) - if flip: - # The following output is not observed - print("FLIP!") - shape = galsim.Shear(g1=-shape.g1, g2=shape.g2) - shape = galsim.Shear(g=shape.g, beta=shape.beta + theta) - shape = shear + shape - - return shape.g1, shape.g2, sig_tmp - - class Loc2Glob(object): r"""Change from local to global coordinates. @@ -555,7 +538,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: @@ -564,60 +546,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( @@ -628,13 +560,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, @@ -661,10 +593,9 @@ 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]) - new_flag_psf = np.zeros_like(psf_file[mod]) - new_flag_star = np.zeros_like(psf_file[mod]) for ccd_id in range(40): m_ccd_id = new_ccd_id == ccd_id if sum(m_ccd_id) == 0: @@ -681,73 +612,8 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): - y_shift ) - header_file_path = ( - self._params["sub_dir_setools"] - + self._params["file_pattern_psfint"] - + f"{exp_name}-{ccd_id}.fits" - ) - try: - header_file = fits.getdata(header_file_path, 1) - except Exception: - continue - header = fits.Header.fromstring( - "\n".join(header_file[0][0]), sep="\n" - ) - wcs = galsim.AstropyWCS(header=header) - - g1_psf_tmp_l = [] - g2_psf_tmp_l = [] - sig_psf_tmp_l = [] - g1_star_tmp_l = [] - g2_star_tmp_l = [] - sig_star_tmp_l = [] - flag_psf_tmp_l = [] - flag_star_tmp_l = [] - - for obj in psf_file[m_ccd_id]: - try: - jac = wcs.jacobian( - world_pos=galsim.CelestialCoord( - ra=obj["RA_LIST"] * galsim.degrees, - dec=obj["DEC_LIST"] * galsim.degrees, - ) - ) - except Exception: - flag_star_tmp_l.append(16) - flag_psf_tmp_l.append(16) - g1_psf_tmp_l.append(0) - g2_psf_tmp_l.append(0) - sig_psf_tmp_l.append(0) - g1_star_tmp_l.append(0) - g2_star_tmp_l.append(0) - sig_star_tmp_l.append(0) - continue - g1_psf_tmp, g2_psf_tmp, sig_psf_tmp = transform_shape( - obj["PSF_MOM_LIST"], jac - ) - - g1_psf_tmp_l.append(g1_psf_tmp) - g2_psf_tmp_l.append(g2_psf_tmp) - sig_psf_tmp_l.append(sig_psf_tmp) - flag_psf_tmp_l.append(obj["PSF_MOM_LIST"][3]) - - g1_star_tmp, g2_star_tmp, sig_star_tmp = ( - transform_shape(obj["STAR_MOM_LIST"], jac) - ) - g1_star_tmp_l.append(g1_star_tmp) - g2_star_tmp_l.append(g2_star_tmp) - sig_star_tmp_l.append(sig_star_tmp) - flag_star_tmp_l.append(obj["STAR_MOM_LIST"][3]) - - new_e1_psf[m_ccd_id] = g1_psf_tmp_l - new_e2_psf[m_ccd_id] = g2_psf_tmp_l - new_sig_psf[m_ccd_id] = sig_psf_tmp_l - new_flag_psf[m_ccd_id] = flag_psf_tmp_l - new_e1_star[m_ccd_id] = g1_star_tmp_l - new_e2_star[m_ccd_id] = g2_star_tmp_l - new_sig_star[m_ccd_id] = sig_star_tmp_l - new_flag_star[m_ccd_id] = flag_star_tmp_l - + # MCCD moment lists carry [e1, e2, sigma, flag]; ellipticities + # and sizes are passed through unrotated (see note above). exp_cat = np.array( list( map( @@ -758,13 +624,13 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): new_y, psf_file["RA_LIST"], psf_file["DEC_LIST"], - new_e1_psf, - new_e2_psf, - new_sig_psf, + psf_file["PSF_MOM_LIST"][:, 0], + psf_file["PSF_MOM_LIST"][:, 1], + psf_file["PSF_MOM_LIST"][:, 2], psf_file["PSF_MOM_LIST"][:, 3], - new_e1_star, - new_e2_star, - new_sig_star, + psf_file["STAR_MOM_LIST"][:, 0], + psf_file["STAR_MOM_LIST"][:, 1], + psf_file["STAR_MOM_LIST"][:, 2], psf_file["STAR_MOM_LIST"][:, 3], new_ccd_id, psf_file["GLOB_POSITION_IMG_LIST"][:, 0], From 071a604427b5c8ab43d7e3b6415edab5e79da7d6 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 12 Jul 2026 23:38:13 +0200 Subject: [PATCH 3/5] test: HSM sky-coords moments match old Jacobian rotation Equivalence test for use_sky_coords=True: draw an elliptical Gaussian on a stamp with a nontrivial local WCS (rotation + shear + 0.187 arcsec/pix scale, positive determinant), and assert that measuring in the pixel frame then applying the removed convert_psf_pix2world.transform_shape rotation (copied in verbatim as reference) agrees with measuring directly in sky coordinates, to numerical precision, for e1/e2 and sigma. A second test pins the deliberate unit change: the sky-frame sigma is the pixel sigma times the WCS scale. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w --- tests/module/test_hsm_sky_coords.py | 128 ++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 tests/module/test_hsm_sky_coords.py diff --git a/tests/module/test_hsm_sky_coords.py b/tests/module/test_hsm_sky_coords.py new file mode 100644 index 000000000..afa5e5167 --- /dev/null +++ b/tests/module/test_hsm_sky_coords.py @@ -0,0 +1,128 @@ +"""UNIT TESTS FOR HSM SKY-COORDINATE MOMENTS. + +Pin the equivalence that motivates measuring PSF/star HSM adaptive moments +directly in world coordinates (``FindAdaptiveMom(use_sky_coords=True)``) instead +of measuring in the pixel frame and rotating afterwards with a WCS Jacobian. + +The reference ``_transform_shape`` below is a verbatim copy of the shape-rotation +math that ``scripts/python/convert_psf_pix2world.py`` used to apply (``transform_shape``, +removed in this change). The test draws an elliptical Gaussian on a stamp with a +nontrivial local WCS (rotation + shear + scale, matching a real CD matrix), then +checks that: + + pixel-frame moments → _transform_shape(local WCS Jacobian) + +agrees with the sky-coordinate moments to numerical precision. + +The local WCS here has positive determinant (no parity flip): production CFHT +runs decomposed with ``flip=False`` (the old ``transform_shape`` flip branch, +which is not a correct galsim reflection, never fired), so the faithful +equivalence to pin is the no-flip one. +""" + +import galsim +import numpy as np +import numpy.testing as npt +import pytest + + +def _transform_shape(mom_list, jac): + """Reference: the old pixel-frame → world shape rotation (verbatim). + + Copied from the removed ``convert_psf_pix2world.transform_shape``. ``mom_list`` + is ``[g1, g2, sigma]`` measured in the pixel frame; ``jac`` is the local WCS. + """ + scale, shear, theta, flip = jac.getDecomposition() + + sig_tmp = mom_list[2] * scale + shape = galsim.Shear(g1=mom_list[0], g2=mom_list[1]) + if flip: + shape = galsim.Shear(g1=-shape.g1, g2=shape.g2) + shape = galsim.Shear(g=shape.g, beta=shape.beta + theta) + shape = shear + shape + + return shape.g1, shape.g2, sig_tmp + + +# A CD-matrix-like local WCS: 0.187 arcsec/pix scale, ~28 deg rotation and a +# small shear, positive determinant (no parity flip). +_SCALE = 0.187 +_THETA = 28.0 * np.pi / 180.0 +_c, _s = np.cos(_THETA), np.sin(_THETA) +_SHEAR = galsim.Shear(g1=0.06, g2=-0.04) +_A = _SHEAR.getMatrix() * _SCALE @ np.array([[_c, -_s], [_s, _c]]) +LOCAL_WCS = galsim.JacobianWCS(_A[0, 0], _A[0, 1], _A[1, 0], _A[1, 1]) + + +@pytest.mark.parametrize( + "e1, e2, sigma_pix", + [ + (0.3, 0.1, 3.0), + (-0.2, 0.25, 2.5), + (0.15, -0.35, 4.0), + (0.0, 0.0, 3.5), + ], +) +def test_sky_coords_moments_match_jacobian_rotation(e1, e2, sigma_pix): + """use_sky_coords=True reproduces pixel-frame moments + Jacobian rotation. + + The local WCS has positive determinant, so ``getDecomposition`` returns + ``flip=False`` and the reference path is a pure scale+rotation+shear — the + transformation ``use_sky_coords`` performs internally. + """ + # Sanity: the reference decomposition sees no parity flip. + assert LOCAL_WCS.getDecomposition()[3] is False + + # An elliptical Gaussian PSF, drawn on a unit-pixel-scale stamp; sigma is in + # pixels because the pixel scale is 1. + stamp = ( + galsim.Gaussian(sigma=sigma_pix) + .shear(e1=e1, e2=e2) + .drawImage(nx=51, ny=51, scale=1.0, method="no_pixel") + ) + + # Pixel-frame measurement, then the old Jacobian rotation. + moms_pix = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array), strict=False + ) + g1_ref, g2_ref, sigma_ref = _transform_shape( + [ + moms_pix.observed_shape.g1, + moms_pix.observed_shape.g2, + moms_pix.moments_sigma, + ], + LOCAL_WCS, + ) + + # Sky-coordinate measurement: same stamp, local WCS attached. + moms_sky = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array, wcs=LOCAL_WCS), + strict=False, + use_sky_coords=True, + ) + + npt.assert_allclose(moms_sky.observed_shape.g1, g1_ref, atol=1e-9, rtol=1e-6) + npt.assert_allclose(moms_sky.observed_shape.g2, g2_ref, atol=1e-9, rtol=1e-6) + npt.assert_allclose(moms_sky.moments_sigma, sigma_ref, rtol=1e-6) + + +def test_sky_coords_sigma_is_scaled_to_world_units(): + """The sky-frame size is the pixel size times the WCS scale (arcsec/pix). + + Pins the deliberate unit change: SIGMA_*_HSM leaves the measurement in world + (arcsec) units rather than pixels once use_sky_coords is on. + """ + stamp = ( + galsim.Gaussian(sigma=3.0) + .drawImage(nx=51, ny=51, scale=1.0, method="no_pixel") + ) + sigma_pix = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array), strict=False + ).moments_sigma + sigma_sky = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array, wcs=LOCAL_WCS), + strict=False, + use_sky_coords=True, + ).moments_sigma + + npt.assert_allclose(sigma_sky, sigma_pix * _SCALE, rtol=1e-6) From c08cdb04330317db8b8873c10c711ef1fb0cca5e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 13 Jul 2026 00:07:09 +0200 Subject: [PATCH 4/5] convert_psf_pix2world: keep HSM rotation for external-mccd MOM_LIST branch The MCCD PSF_MOM_LIST/STAR_MOM_LIST columns come from the external mccd fit-validation code (mccd.auxiliary_fun.mccd_validation), which still measures HSM adaptive moments in the pixel frame. Stripping the WCS-Jacobian rotation there (71d502f7) would have handed MergeStarCatMCCD pixel-frame ellipticities. Restore transform_shape + the per-CCD getDecomposition rotation for the MCCD MOM_LIST branch only (and the galsim import it needs); the in-repo PSFEx and MCCD-interpolation paths stay unrotated, as they now measure in world coordinates upstream. Update the module docstring caveat: this branch keeps the rotation until mccd itself adopts use_sky_coords. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w --- scripts/python/convert_psf_pix2world.py | 142 ++++++++++++++++++++++-- 1 file changed, 130 insertions(+), 12 deletions(-) diff --git a/scripts/python/convert_psf_pix2world.py b/scripts/python/convert_psf_pix2world.py index 72faa34c8..e1adc90df 100755 --- a/scripts/python/convert_psf_pix2world.py +++ b/scripts/python/convert_psf_pix2world.py @@ -12,10 +12,12 @@ 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, which still measures in the pixel frame; -those shapes are passed through unrotated here. If the MCCD validation branch is -revived, its moments must be measured in sky coordinates upstream (as the PSFEx -and MCCD interpolation paths now are). +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 @@ -28,11 +30,46 @@ import numpy as np from astropy.io import fits +import galsim from cs_util import args as cs_args from cs_util import logging +def transform_shape(mom_list, jac): + """Transform Shape. + + Transform shape (ellipticity and size) using a Jacobian. + + Parameters + ---------- + mom_list : list + input moment measurements; each list element contains + first and second ellipticity component and size + jac : galsim.JacobianWCS + Jacobian transformation matrix information + + Returns + ------- + list + transformed shape parameters, which are + first and second ellipticity component and size + + """ + scale, shear, theta, flip = jac.getDecomposition() + + sig_tmp = mom_list[2] * scale + shape = galsim.Shear(g1=mom_list[0], g2=mom_list[1]) + if flip: + # The following output is not observed + print("FLIP!") + shape = galsim.Shear(g1=-shape.g1, g2=shape.g2) + shape = galsim.Shear(g=shape.g, beta=shape.beta + theta) + shape = shear + shape + + return shape.g1, shape.g2, sig_tmp + + class Loc2Glob(object): r"""Change from local to global coordinates. @@ -596,6 +633,22 @@ 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): m_ccd_id = new_ccd_id == ccd_id if sum(m_ccd_id) == 0: @@ -612,8 +665,73 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): - y_shift ) - # MCCD moment lists carry [e1, e2, sigma, flag]; ellipticities - # and sizes are passed through unrotated (see note above). + header_file_path = ( + self._params["sub_dir_setools"] + + self._params["file_pattern_psfint"] + + f"{exp_name}-{ccd_id}.fits" + ) + try: + header_file = fits.getdata(header_file_path, 1) + except Exception: + continue + header = fits.Header.fromstring( + "\n".join(header_file[0][0]), sep="\n" + ) + wcs = galsim.AstropyWCS(header=header) + + g1_psf_tmp_l = [] + g2_psf_tmp_l = [] + sig_psf_tmp_l = [] + g1_star_tmp_l = [] + g2_star_tmp_l = [] + sig_star_tmp_l = [] + flag_psf_tmp_l = [] + flag_star_tmp_l = [] + + for obj in psf_file[m_ccd_id]: + try: + jac = wcs.jacobian( + world_pos=galsim.CelestialCoord( + ra=obj["RA_LIST"] * galsim.degrees, + dec=obj["DEC_LIST"] * galsim.degrees, + ) + ) + except Exception: + flag_star_tmp_l.append(16) + flag_psf_tmp_l.append(16) + g1_psf_tmp_l.append(0) + g2_psf_tmp_l.append(0) + sig_psf_tmp_l.append(0) + g1_star_tmp_l.append(0) + g2_star_tmp_l.append(0) + sig_star_tmp_l.append(0) + continue + g1_psf_tmp, g2_psf_tmp, sig_psf_tmp = transform_shape( + obj["PSF_MOM_LIST"], jac + ) + + g1_psf_tmp_l.append(g1_psf_tmp) + g2_psf_tmp_l.append(g2_psf_tmp) + sig_psf_tmp_l.append(sig_psf_tmp) + flag_psf_tmp_l.append(obj["PSF_MOM_LIST"][3]) + + g1_star_tmp, g2_star_tmp, sig_star_tmp = ( + transform_shape(obj["STAR_MOM_LIST"], jac) + ) + g1_star_tmp_l.append(g1_star_tmp) + g2_star_tmp_l.append(g2_star_tmp) + sig_star_tmp_l.append(sig_star_tmp) + flag_star_tmp_l.append(obj["STAR_MOM_LIST"][3]) + + new_e1_psf[m_ccd_id] = g1_psf_tmp_l + new_e2_psf[m_ccd_id] = g2_psf_tmp_l + new_sig_psf[m_ccd_id] = sig_psf_tmp_l + new_flag_psf[m_ccd_id] = flag_psf_tmp_l + new_e1_star[m_ccd_id] = g1_star_tmp_l + new_e2_star[m_ccd_id] = g2_star_tmp_l + new_sig_star[m_ccd_id] = sig_star_tmp_l + new_flag_star[m_ccd_id] = flag_star_tmp_l + exp_cat = np.array( list( map( @@ -624,13 +742,13 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): new_y, psf_file["RA_LIST"], psf_file["DEC_LIST"], - psf_file["PSF_MOM_LIST"][:, 0], - psf_file["PSF_MOM_LIST"][:, 1], - psf_file["PSF_MOM_LIST"][:, 2], + new_e1_psf, + new_e2_psf, + new_sig_psf, psf_file["PSF_MOM_LIST"][:, 3], - psf_file["STAR_MOM_LIST"][:, 0], - psf_file["STAR_MOM_LIST"][:, 1], - psf_file["STAR_MOM_LIST"][:, 2], + new_e1_star, + new_e2_star, + new_sig_star, psf_file["STAR_MOM_LIST"][:, 3], new_ccd_id, psf_file["GLOB_POSITION_IMG_LIST"][:, 0], From f153ea2896fdce65013754704aba4d31e418ee0d Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 13 Jul 2026 00:07:18 +0200 Subject: [PATCH 5/5] psf interp: use origin=1 positions for local_wcs_list MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit local_wcs_list wraps galsim.AstropyWCS, which follows the 1-indexed FITS pixel convention, but the multi-epoch paths fed it positions from all_world2pix(..., 0) (origin 0) — evaluating the local Jacobian one pixel off. Compute a separate origin=1 position array (gal_pos_wcs) for local_wcs_list; the origin-0 gal_pos stays for the PSF-model interpolation (interpsfex / interp_MCCD), which is its only other use. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w --- .../modules/mccd_package/mccd_interpolation_script.py | 11 ++++++++++- .../modules/psfex_interp_package/psfex_interp.py | 11 ++++++++++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py b/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py index 468cf55ee..a01ca095d 100644 --- a/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py +++ b/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py @@ -485,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) @@ -533,7 +542,7 @@ def _interpolate_me(self): if self._compute_shape: self._get_psfshapes( local_wcs_list( - self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos + self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos_wcs ) ) if array_shape is None: diff --git a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py index 96b2c2007..b6da6623f 100644 --- a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py +++ b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py @@ -750,6 +750,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 PSFEx 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 = self.interpsfex( dot_psf_path, @@ -796,7 +805,7 @@ def _interpolate_me(self): if self._compute_shape: self._get_psfshapes( local_wcs_list( - self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos + self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos_wcs ) ) if array_shape is None: