refactor(pixel-scale): derive from image WCS instead of hardcoding it

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

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"

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

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
 

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
 

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

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
 

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)

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'

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

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

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,34 +98,32 @@ 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
     -------
     tuple
         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,15 +162,13 @@ def __init__(
         psf_cat_path,
         weight_cat_path,
         output_path,
-        pixel_scale,
         output_mode,
     ):
 
         self._sex_cat_path = sex_cat_path
         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(

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,
     )