From c041cb3f83de41eddea1289e9141f9b5031ccd69 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Mon, 20 Jul 2026 22:52:50 +0200 Subject: [PATCH] refactor(pixel-scale): derive from image WCS, drop redundant config knobs The pixel scale lived in four places that could silently drift from the actual image astrometry: the ngmix and spread_model config values, a dead MESSIER_PIXEL_SCALE mask knob, and the images' own WCS. SExtractor and the mask code already read it from the WCS; only ngmix and spread_model carried a hardcoded constant. Unify on the image WCS as the single source of truth. - ngmix: derive the pixel scale from the image WCS when PIXEL_SCALE is absent or non-positive (mirrors SExtractor's `PIXEL_SCALE 0`). It stays an optional override. Only the centroid-prior width and noise window use it; the fit Jacobian is already built per object from the full WCS. - spread_model: drop the parameter entirely. The statistic is computed in pixel space, so the scale cancels exactly (verified byte-for-byte); the models now render in pixel units. - mask: delete MESSIER_PIXEL_SCALE from the configs -- nothing reads it (mask.py derives the Messier radius from the WCS pixel_scale_matrix). - configs: remove every hardcoded PIXEL_SCALE / MESSIER_PIXEL_SCALE from the example configs. Nothing left to drift. Supersedes the 0.185 hardcode (#858). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CQxJdyKjKSbHRudXLHfXZb --- .../cfis/config_tile_Ng_batch_psfex_uc.ini | 3 - example/cfis/config_tile_Ng_template.ini | 3 - .../cfis/config_tile_Ng_template_batch.ini | 3 - example/cfis/config_tile_onthefly.mask | 1 - example/cfis/config_tile_save.mask | 1 - .../config_tile_Ng_batch_psfex_sx.ini | 3 - .../config_tile_onthefly.mask_simu | 1 - .../modules/ngmix_package/__init__.py | 7 +- src/shapepipe/modules/ngmix_package/ngmix.py | 66 ++++++++++++++++++- src/shapepipe/modules/ngmix_runner.py | 9 ++- .../modules/spread_model_package/__init__.py | 2 - .../spread_model_package/spread_model.py | 29 ++++---- src/shapepipe/modules/spread_model_runner.py | 5 +- 13 files changed, 91 insertions(+), 42 deletions(-) diff --git a/example/cfis/config_tile_Ng_batch_psfex_uc.ini b/example/cfis/config_tile_Ng_batch_psfex_uc.ini index b223e05a1..f2696bd8f 100644 --- a/example/cfis/config_tile_Ng_batch_psfex_uc.ini +++ b/example/cfis/config_tile_Ng_batch_psfex_uc.ini @@ -73,9 +73,6 @@ BKG_RMS_VIGNET_PATH = $SP_RUN/output/run_sp_tile_PiViVi/vignetmaker_runner_run_2 # Magnitude zero-point MAG_ZP = 30.0 -# Pixel scale in arcsec -PIXEL_SCALE = 0.186 - SAVE_BATCH = 1000 ID_OBJ_MIN = -1 diff --git a/example/cfis/config_tile_Ng_template.ini b/example/cfis/config_tile_Ng_template.ini index 27f2a5aa2..bd2af3ea1 100644 --- a/example/cfis/config_tile_Ng_template.ini +++ b/example/cfis/config_tile_Ng_template.ini @@ -73,9 +73,6 @@ BKG_RMS_VIGNET_PATH = $SP_RUN/output/run_sp_tile_PiViVi/vignetmaker_runner_run_2 # Magnitude zero-point MAG_ZP = 30.0 -# Pixel scale in arcsec -PIXEL_SCALE = 0.186 - # CENTROID_SOURCE (optional): how to place the galaxy Jacobian origin for the # centroid prior. "wcs" (default) uses the catalog sky position projected # through the WCS, trusting the astrometry — the recommended choice. "hsm" diff --git a/example/cfis/config_tile_Ng_template_batch.ini b/example/cfis/config_tile_Ng_template_batch.ini index f15849674..655397156 100644 --- a/example/cfis/config_tile_Ng_template_batch.ini +++ b/example/cfis/config_tile_Ng_template_batch.ini @@ -77,8 +77,5 @@ SAVE_BATCH = 1000 # Magnitude zero-point MAG_ZP = 30.0 -# Pixel scale in arcsec -PIXEL_SCALE = 0.186 - ID_OBJ_MIN = X ID_OBJ_MAX = X diff --git a/example/cfis/config_tile_onthefly.mask b/example/cfis/config_tile_onthefly.mask index 18cf5db2d..69ad20769 100644 --- a/example/cfis/config_tile_onthefly.mask +++ b/example/cfis/config_tile_onthefly.mask @@ -51,7 +51,6 @@ SPIKE_REG_FILE = spike.reg MESSIER_MAKE = True MESSIER_CAT_PATH = $SP_CONFIG/mask_default/Messier_catalog_updated.fits -MESSIER_PIXEL_SCALE = 0.187 MESSIER_SIZE_PLUS = 0. MESSIER_FLAG_VALUE = 16 diff --git a/example/cfis/config_tile_save.mask b/example/cfis/config_tile_save.mask index 1de738301..82c4b66af 100644 --- a/example/cfis/config_tile_save.mask +++ b/example/cfis/config_tile_save.mask @@ -51,7 +51,6 @@ SPIKE_REG_FILE = spike.reg MESSIER_MAKE = True MESSIER_CAT_PATH = $SP_CONFIG/mask_default/Messier_catalog_updated.fits -MESSIER_PIXEL_SCALE = 0.187 MESSIER_SIZE_PLUS = 0. MESSIER_FLAG_VALUE = 16 diff --git a/example/cfis_image_sims/config_tile_Ng_batch_psfex_sx.ini b/example/cfis_image_sims/config_tile_Ng_batch_psfex_sx.ini index 11f11cea6..8ca298565 100644 --- a/example/cfis_image_sims/config_tile_Ng_batch_psfex_sx.ini +++ b/example/cfis_image_sims/config_tile_Ng_batch_psfex_sx.ini @@ -70,9 +70,6 @@ MAG_ZP = 30.0 # No background subtraction for image sims (background not simulated) BKG_SUB = False -# Pixel scale in arcsec -PIXEL_SCALE = 0.186 - SAVE_BATCH = 1000 # Position-seeded per-object fixnoise RNG (#796/#803): the SKiLLS shear diff --git a/example/cfis_image_sims/config_tile_onthefly.mask_simu b/example/cfis_image_sims/config_tile_onthefly.mask_simu index f42d05d39..2e404ed96 100644 --- a/example/cfis_image_sims/config_tile_onthefly.mask_simu +++ b/example/cfis_image_sims/config_tile_onthefly.mask_simu @@ -51,7 +51,6 @@ SPIKE_REG_FILE = spike.reg MESSIER_MAKE = False MESSIER_CAT_PATH = $SP_CONFIG/mask_default/Messier_catalog_updated.fits -MESSIER_PIXEL_SCALE = 0.187 MESSIER_SIZE_PLUS = 0. MESSIER_FLAG_VALUE = 16 diff --git a/src/shapepipe/modules/ngmix_package/__init__.py b/src/shapepipe/modules/ngmix_package/__init__.py index 230606bba..89727caff 100644 --- a/src/shapepipe/modules/ngmix_package/__init__.py +++ b/src/shapepipe/modules/ngmix_package/__init__.py @@ -33,8 +33,11 @@ MAG_ZP : float Photometric zero point -PIXEL_SCALE : float - Pixel scale in arcseconds +PIXEL_SCALE : float, optional + Pixel scale in arcsec. Optional override; when omitted (or non-positive) + it is read from the image WCS so it cannot drift from the pixels. Only + sets the centroid-prior width and noise window -- the fit Jacobian is + built per object from the full WCS. SAVE_BATCH : int, optional Save the output catalogue in batches of this size; default is ``-1`` (no batch saving) diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 1403d4d0d..e07e0dd7b 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -14,6 +14,7 @@ import galsim import numpy as np from astropy.io import fits +from astropy.wcs.utils import proj_plane_pixel_scales from cs_util import size as cs_size from modopt.math.stats import sigma_mad from ngmix.observation import Observation, ObsList @@ -352,6 +353,44 @@ def close(self): if self.bkg_rms_vign_cat is not None: self.bkg_rms_vign_cat.close() +def pixel_scale_from_wcs(f_wcs_file): + """Representative pixel scale (arcsec) from the tile's image WCS. + + The ngmix fit builds each object's Jacobian from the full per-epoch WCS, + so this scalar only sets the centroid-prior width and the noise-window + scale (see :func:`get_prior`, :func:`get_noise`). A single value read from + the astrometry is therefore sufficient -- and, unlike a hard-coded config + constant, it cannot silently drift from the pixels it describes (this + mirrors SExtractor's ``PIXEL_SCALE 0`` convention). + + Parameters + ---------- + f_wcs_file : dict-like + Mapping ``exposure -> {ccd -> {"WCS": astropy.wcs.WCS, ...}}`` (the + merged single-exposure headers opened by :class:`Vignet`). + + Returns + ------- + float + Pixel scale in arcsec, averaged over the two axes of the first WCS. + + Raises + ------ + ValueError + If no WCS can be found in ``f_wcs_file``. + """ + for exp_name in f_wcs_file: + for ccd_dict in f_wcs_file[exp_name].values(): + # proj_plane_pixel_scales returns deg/pixel per world axis; the + # two axes are equal to sub-per-mille for survey astrometry. + scales = proj_plane_pixel_scales(ccd_dict["WCS"]) + return float(np.mean(scales) * 3600.0) + raise ValueError( + "cannot derive PIXEL_SCALE from the WCS: the merged single-exposure " + "headers file contains no exposures" + ) + + class Ngmix(object): """Ngmix. @@ -367,8 +406,10 @@ class Ngmix(object): File numbering scheme zero_point : float Photometric zero point - pixel_scale : float - Pixel scale in arcsec + pixel_scale : float or None + Pixel scale in arcsec. Optional override: when ``None`` or + non-positive, the pixel scale is derived from the image WCS (see + :func:`pixel_scale_from_wcs`). f_wcs_path : str Path to merged single-exposure single-HDU headers w_log : logging.Logger @@ -492,7 +533,6 @@ def __init__( self._file_number_string = file_number_string self._zero_point = zero_point - self._pixel_scale = pixel_scale self._f_wcs_path = f_wcs_path @@ -513,6 +553,26 @@ def __init__( seed = int(''.join(re.findall(r'\d+', self._file_number_string))) self._rng = np.random.RandomState(seed) self._w_log.info(f'Random generator initialisation seed = {seed}') + + # Pixel scale: an explicit positive PIXEL_SCALE overrides; otherwise + # derive it from the image WCS so it can never drift from the pixels + # (mirrors SExtractor's ``PIXEL_SCALE 0`` convention). Only the + # centroid-prior width and noise window use it -- the fit Jacobian is + # built per object from the full WCS. + if pixel_scale is None or pixel_scale <= 0: + self._pixel_scale = pixel_scale_from_wcs( + self._vignet_cat.f_wcs_file + ) + self._w_log.info( + f'PIXEL_SCALE derived from image WCS = ' + f'{self._pixel_scale:.6f} arcsec' + ) + else: + self._pixel_scale = pixel_scale + self._w_log.info( + f'PIXEL_SCALE from config = {self._pixel_scale:.6f} arcsec' + ) + if self._seed_from_position: self._w_log.info( 'SEED_FROM_POSITION on: per-object RNG seeded from sky position' diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index d97f6739e..16403a543 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -48,8 +48,13 @@ def ngmix_runner( # Photometric zero point zero_point = config.getfloat(module_config_sec, "MAG_ZP") - # Pixel scale - pixel_scale = config.getfloat(module_config_sec, "PIXEL_SCALE") + # Pixel scale -- optional override. When absent (or non-positive) it is + # derived from the image WCS inside Ngmix, so it cannot drift from the + # pixels. Only the centroid-prior width and noise window use it. + if config.has_option(module_config_sec, "PIXEL_SCALE"): + pixel_scale = config.getfloat(module_config_sec, "PIXEL_SCALE") + else: + pixel_scale = None # Background subtraction: disable for image sims, where there is no # background image to subtract. The background vignet occupies one input diff --git a/src/shapepipe/modules/spread_model_package/__init__.py b/src/shapepipe/modules/spread_model_package/__init__.py index 336283e78..247f60308 100644 --- a/src/shapepipe/modules/spread_model_package/__init__.py +++ b/src/shapepipe/modules/spread_model_package/__init__.py @@ -25,8 +25,6 @@ PREFIX : str, optional Ouput file prefix -PIXEL_SCALE : float - Pixel scale in arcseconds OUTPUT_MODE : str Output mode, options are ``add`` to add the outputs to the input SExtractor catalogue or ``new`` to generte a new output catalogue diff --git a/src/shapepipe/modules/spread_model_package/spread_model.py b/src/shapepipe/modules/spread_model_package/spread_model.py index ed72aa58f..57d88dc43 100644 --- a/src/shapepipe/modules/spread_model_package/spread_model.py +++ b/src/shapepipe/modules/spread_model_package/spread_model.py @@ -77,13 +77,19 @@ def get_sm(obj_vign, psf_vign, model_vign, weight_vign): return sm, sm_err -def get_model(sigma, flux, img_shape, pixel_scale=0.186): +def get_model(sigma, flux, img_shape): """Get Model. This method computes - an exponential galaxy model with scale radius = 1/16 FWHM - a Gaussian model for the PSF + The spread-model statistic compares the object, PSF, and galaxy-model + vignets entirely in pixel space, so the models are rendered in pixel + units. A physical pixel scale would cancel out exactly (every length is + ``pixels x pixel_scale``, drawn at ``scale=pixel_scale``), so none is + needed here. + Parameters ---------- sigma : float @@ -92,8 +98,6 @@ def get_model(sigma, flux, img_shape, pixel_scale=0.186): Flux of the galaxy for the model img_shape : list Size of the output vignet ``[xsize, ysize]`` - pixel_scale : float, optional - Pixel scale to use for the model (in arcsec); default is ``0.186`` Returns ------- @@ -101,25 +105,25 @@ def get_model(sigma, flux, img_shape, pixel_scale=0.186): Vignet of the galaxy model and of the PSF model """ - # Get scale radius - scale_radius = 1 / 16 * cs_size.sigma_to_fwhm(sigma) * pixel_scale + # Get scale radius (pixels) + scale_radius = 1 / 16 * cs_size.sigma_to_fwhm(sigma) # Get galaxy model gal_obj = galsim.Exponential(scale_radius=scale_radius, flux=flux) - # Get PSF - psf_obj = galsim.Gaussian(sigma=sigma * pixel_scale) + # Get PSF (sigma in pixels) + psf_obj = galsim.Gaussian(sigma=sigma) # Convolve both gal_obj = galsim.Convolve(gal_obj, psf_obj) - # Draw galaxy and PSF on vignets + # Draw galaxy and PSF on vignets, in pixel units (scale = 1 pixel) gal_vign = gal_obj.drawImage( - nx=img_shape[0], ny=img_shape[1], scale=pixel_scale + nx=img_shape[0], ny=img_shape[1], scale=1.0 ).array psf_vign = psf_obj.drawImage( - nx=img_shape[0], ny=img_shape[1], scale=pixel_scale + nx=img_shape[0], ny=img_shape[1], scale=1.0 ).array return gal_vign, psf_vign @@ -138,8 +142,6 @@ class SpreadModel(object): Path to weight catalogue output_path : str Output file path of pasted catalog - pixel_scale : float - Pixel scale in arcsec output_mode : str Options are ``new`` or ``add`` @@ -160,7 +162,6 @@ def __init__( psf_cat_path, weight_cat_path, output_path, - pixel_scale, output_mode, ): @@ -168,7 +169,6 @@ def __init__( self._psf_cat_path = psf_cat_path self._weight_cat_path = weight_cat_path self._output_path = output_path - self._pixel_scale = pixel_scale self._output_mode = output_mode def process(self): @@ -228,7 +228,6 @@ def process(self): obj_sigma_tmp, obj_flux_tmp, obj_vign_tmp.shape, - self._pixel_scale, ) obj_sm, obj_sm_err = get_sm( diff --git a/src/shapepipe/modules/spread_model_runner.py b/src/shapepipe/modules/spread_model_runner.py index 16c9a2d3b..93b70b753 100644 --- a/src/shapepipe/modules/spread_model_runner.py +++ b/src/shapepipe/modules/spread_model_runner.py @@ -44,8 +44,8 @@ def spread_model_runner( else: prefix = "" - # Get pixel scale and output mode - pixel_scale = config.getfloat(module_config_sec, "PIXEL_SCALE") + # Get output mode. The spread model is computed entirely in pixel space, + # so no pixel scale is needed (see spread_model.get_model). output_mode = config.get(module_config_sec, "OUTPUT_MODE") # Set output file path @@ -58,7 +58,6 @@ def spread_model_runner( psf_cat_path, weight_cat_path, output_path, - pixel_scale, output_mode, )