From e8f5ff744557fdc43c70514ba483dcf1e3609b34 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sun, 12 Jul 2026 23:34:13 +0200 Subject: [PATCH 1/6] 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/6] 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/6] 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/6] 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/6] 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: From baea03a3608058f46c0d20ae4945d02c72153f3e Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 20 Jul 2026 02:44:03 +0200 Subject: [PATCH 6/6] refactor(scripts): rename convert_psf_pix2world -> collate_star_cat, add patch-less v2.0 Re-scope the script to its actual remaining purpose: collating per-exposure PSF validation catalogues into star catalogues. The WCS-Jacobian shape rotation is gone (HSM moments are measured in sky coordinates upstream), so "convert" no longer describes it. The external-mccd PSF_MOM_LIST/STAR_MOM_LIST rotation branch stays until mccd adopts use_sky_coords. Add a -V/--version option mirroring get_ccds_with_psf. v1.x is unchanged (loop over -P patches, /P/output, validation_psf_conv--). v2.0 is patch-less: single /output root, patch token dropped from the output filename (validation_psf_conv-). The downstream combine_runs.bash glob (validation_psf_conv-*) matches both. Path/filename logic extracted into pure collate_paths/output_filename helpers with unit tests; docs collate section updated with both invocations. The script command name follows the file (Dockerfile symlinks by basename). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WqNFhq9uGpDpprYeuGeb8k --- docs/source/pipeline_canfar.md | 32 +++-- ...t_psf_pix2world.py => collate_star_cat.py} | 117 +++++++++++++++--- tests/module/test_collate_star_cat.py | 70 +++++++++++ tests/module/test_hsm_sky_coords.py | 4 +- 4 files changed, 197 insertions(+), 26 deletions(-) rename scripts/python/{convert_psf_pix2world.py => collate_star_cat.py} (87%) create mode 100644 tests/module/test_collate_star_cat.py diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index a3ec03e7b..96add5758 100644 --- a/docs/source/pipeline_canfar.md +++ b/docs/source/pipeline_canfar.md @@ -333,27 +333,43 @@ shapepipe_run -c $SP_CONFIG/config_Ms_$psf.ini shapepipe_run -c $SP_CONFIG/config_Pl_$psf.ini ``` -#### Collate star catalogues per patch +#### Collate star catalogues -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. +Collate all input validation PSF files into star catalogues, gathering 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. +`collate_star_cat.py` only collates and passes the shapes through. + +> **v2.0 is patch-less.** Runs up to `v1.6` are organised in sky patches +> `P1`..`P`, each patch a run directory. `v2.0` removes the patch concept: a +> single run root, outputs directly under it. `v2.0` is the default; select an +> older layout with `-V` (e.g. `-V v1.6`). For `v2.0` the patch loop and the +> `-P` option no longer apply, and the patch token drops from the output +> filename (`validation_psf_conv-.fits` instead of +> `validation_psf_conv--.fits`). ```bash cd /path/to/version -mkdir stat_car +mkdir star_cat cd star_cat ``` -For each patch run +For `v2.0` (the default), run once against the patch-less run root, producing +files `validation_psf_conv-.fits`: + +```bash +collate_star_cat.py -i .. -v +``` + +For `v1.x`, pass the version explicitly and run once per patch, creating a +directory per patch `P?` and producing files +`validation_psf_conv--.fits` (for the v1.4 setup only one file): ```bash -convert_psf_pix2world.py -i .. -P $patchnum -v +collate_star_cat.py -i .. -V v1.6 -P $patchnum -v ``` Combine previously created files as links within one ShapePipe run directory (for the v1.4 setup only one link). diff --git a/scripts/python/convert_psf_pix2world.py b/scripts/python/collate_star_cat.py similarity index 87% rename from scripts/python/convert_psf_pix2world.py rename to scripts/python/collate_star_cat.py index e1adc90df..f7ab25b41 100755 --- a/scripts/python/convert_psf_pix2world.py +++ b/scripts/python/collate_star_cat.py @@ -1,10 +1,17 @@ #! /usr/bin/env python3 -"""CONVERT PSF CATALOGUE. +"""COLLATE STAR CATALOGUES. -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``. +Collate the per-exposure PSF validation catalogues into star catalogues: gather +positions (X/Y/RA/DEC), assign the MCCD focal-plane CCD id, and merge into +``validation_psf_conv`` FITS files. + +Catalogue layout depends on the major version (``-V``). Runs up to ``v1.6`` are +organised in sky patches (``P1``..``P``): input runs live under +``/P/output`` and outputs are named +``validation_psf_conv--.fits``. ``v2.0`` removes the patch concept: +input runs live under a single ``/output`` root and outputs drop the patch +token (``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 @@ -36,6 +43,58 @@ from cs_util import logging +def collate_paths(input_base_dir, output_base_dir, patch): + """Collate Paths. + + Return the ``(input run dir, output dir)`` for a patch. ``patch`` is None + for the patch-less v2.0 layout, which drops the ``P`` token; v1.x + passes the patch number. + + Parameters + ---------- + input_base_dir : str + input base directory + output_base_dir : str + output base directory + patch : str or None + patch number, or None for the patch-less v2.0 layout + + Returns + ------- + tuple + input run directory and output directory + + """ + if patch is None: + return f"{input_base_dir}/output/", output_base_dir + return f"{input_base_dir}/P{patch}/output/", f"{output_base_dir}/P{patch}" + + +def output_filename(file_pattern, patch, idx): + """Output Filename. + + Build the collated catalogue filename. ``patch`` is None for the patch-less + v2.0 layout, which drops the patch token. + + Parameters + ---------- + file_pattern : str + input file pattern (e.g. ``validation_psf``) + patch : str or None + patch number, or None for the patch-less v2.0 layout + idx : int + exposure run index + + Returns + ------- + str + output catalogue file name + + """ + patch_token = "" if patch is None else f"{patch}-" + return f"{file_pattern}_conv-{patch_token}{idx}.fits" + + def transform_shape(mom_list, jac): """Transform Shape. @@ -398,6 +457,7 @@ def params_default(self): self._params = { "input_base_dir": ".", "output_base_dir": ".", + "version_cat": "v2.0", "mode": "merge", "patches": "", "psf": "psfex", @@ -406,6 +466,7 @@ def params_default(self): self._short_options = { "input_base_dir": "-i", + "version_cat": "-V", "mode": "-m", "psf": "-p", "patches": "-P", @@ -415,15 +476,19 @@ def params_default(self): self._help_strings = { "input_base_dir": ( - "input base dir, runs are expected in" - + " /P/output;" - " default is {}" + "input base dir; for v1.x runs are expected in" + + " /P/output, for v2.0 (patch-less) in" + + " /output; default is {}" + ), + "version_cat": ( + "catalogue major version, allowed are v1.3, v1.4, v1.5, v1.6," + + " v2.0; v2.0 is patch-less; default is {}" ), "mode": ( "run mode, allowed are 'merge', 'test'; default is" + " '{}'" ), "psf": "PSF model, allowed are 'psfex' and 'mccd'; default is {}", - "patches": "(list of) input patches", + "patches": "(list of) input patches; ignored for v2.0", } # Output column names with types @@ -474,7 +539,22 @@ def run(self): Main processing function. """ - if self._params["mode"] == "test": + # Guard against a mistyped version silently falling through to the + # v1.x patch loop (e.g. ``-V v2`` or ``-V 2.0``). + allowed_versions = ("v1.3", "v1.4", "v1.5", "v1.6", "v2.0") + if self._params["version_cat"] not in allowed_versions: + raise ValueError( + f"Invalid version {self._params['version_cat']}; allowed are" + + f" {', '.join(allowed_versions)}" + ) + + # v2.0 removes the patch concept: a single patch-less run root. For + # v1.x, iterate over the requested sky patches as before. ``patch`` is + # None in the patch-less case, which drops the patch token from the + # input path and the output filename. + if self._params["version_cat"] == "v2.0": + patch_nums = [None] + elif self._params["mode"] == "test": patch_nums = ["3", "4"] else: patch_nums = cs_args.my_string_split(self._params["patches"]) @@ -484,13 +564,16 @@ def run(self): # Loop over patches for patch in patch_nums: - output_dir = f"{self._params['output_base_dir']}/P{patch}" - if not os.path.isdir(output_dir): - os.makedirs(output_dir, exist_ok=False) + patch_dir, output_dir = collate_paths( + self._params["input_base_dir"], + self._params["output_base_dir"], + patch, + ) + print("Running patch-less (v2.0)" if patch is None else f"Running patch: {patch}") - print("Running patch:", patch) + if not os.path.isdir(output_dir): + os.makedirs(output_dir, exist_ok=True) - patch_dir = f"{self._params['input_base_dir']}/P{patch}/output/" subdirs = f"{patch_dir}/{self._params['sub_dir_pattern']}*" exp_run_dirs = glob.glob(subdirs) n_exp_runs = len(exp_run_dirs) @@ -536,8 +619,10 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir): """ output_path = ( - f"{output_dir}/{self._params['file_pattern_psfint']}_conv" - + f"-{patch}-{idx}.fits" + f"{output_dir}/" + + output_filename( + self._params["file_pattern_psfint"], patch, idx + ) ) if os.path.exists(output_path): print(f"Skipping transform_exposures, file {output_path} exists") diff --git a/tests/module/test_collate_star_cat.py b/tests/module/test_collate_star_cat.py new file mode 100644 index 000000000..20d52ce23 --- /dev/null +++ b/tests/module/test_collate_star_cat.py @@ -0,0 +1,70 @@ +"""UNIT TESTS FOR STAR-CATALOGUE COLLATION PATHS. + +Pin the patch vs patch-less (v2.0) path and filename convention of +``scripts/python/collate_star_cat.py``. Runs up to v1.6 carry a ``P`` +token in both the input run directory and the output filename; v2.0 is +patch-less (``patch is None``) and drops that token, reading from a single +``/output`` root and writing ``validation_psf_conv-.fits`` — the +name still matched by the downstream ``validation_psf_conv-*`` glob. +""" + +import importlib.util +from pathlib import Path + +import pytest + +# The collation script lives under scripts/python (not an importable package), +# so load it by path. +_SCRIPT = ( + Path(__file__).resolve().parents[2] + / "scripts" + / "python" + / "collate_star_cat.py" +) +_spec = importlib.util.spec_from_file_location("collate_star_cat", _SCRIPT) +collate_star_cat = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(collate_star_cat) + + +@pytest.mark.parametrize( + "patch, exp_input, exp_output", + [ + ("3", "in/P3/output/", "out/P3"), + (None, "in/output/", "out"), + ], +) +def test_collate_paths(patch, exp_input, exp_output): + """v1.x carries the P token; v2.0 (patch None) drops it.""" + assert collate_star_cat.collate_paths("in", "out", patch) == ( + exp_input, + exp_output, + ) + + +@pytest.mark.parametrize( + "patch, expected", + [ + ("3", "validation_psf_conv-3-0.fits"), + (None, "validation_psf_conv-0.fits"), + ], +) +def test_output_filename(patch, expected): + """The patch token is present for v1.x and absent for v2.0.""" + assert collate_star_cat.output_filename("validation_psf", patch, 0) == expected + + +def test_output_filename_matches_downstream_glob(): + """Both layouts stay under the downstream ``validation_psf_conv-*`` glob.""" + for patch in ("1", None): + assert collate_star_cat.output_filename( + "validation_psf", patch, 5 + ).startswith("validation_psf_conv-") + + +@pytest.mark.parametrize("bad", ["v2", "2.0", "v1.7", ""]) +def test_invalid_version_raises(bad): + """A mistyped -V is rejected rather than falling through to v1.x.""" + obj = collate_star_cat.Convert() + obj._params["version_cat"] = bad + with pytest.raises(ValueError): + obj.run() diff --git a/tests/module/test_hsm_sky_coords.py b/tests/module/test_hsm_sky_coords.py index afa5e5167..3d19c0af6 100644 --- a/tests/module/test_hsm_sky_coords.py +++ b/tests/module/test_hsm_sky_coords.py @@ -5,7 +5,7 @@ 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``, +math that ``scripts/python/collate_star_cat.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: @@ -29,7 +29,7 @@ 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`` + Copied from the removed ``collate_star_cat.transform_shape``. ``mom_list`` is ``[g1, g2, sigma]`` measured in the pixel frame; ``jac`` is the local WCS. """ scale, shear, theta, flip = jac.getDecomposition()