Cut postage stamps in-house and reuse their exact centroid in ngmix (drop sf_tools)

docs/source/dependencies.md

@@ -28,7 +28,6 @@ Scientific stack:
 | [mpi4py](https://mpi4py.readthedocs.io/en/stable/) | {cite:p}`dalcin:05,dalcin:08,dalcin:11` |
 | [reproject](https://reproject.readthedocs.io/) | |
 | [h5py](https://www.h5py.org/) | |
-| [sf_tools](https://github.com/sfarrens/sf_tools) | |
 
 Data access & infrastructure (CANFAR / UNIONS):
 

pyproject.toml

@@ -36,7 +36,6 @@ dependencies = [
     "PyQt5",
     "pyqtgraph",
     "reproject>=0.19",
-    "sf_tools>=2.0.4",
     "skaha>=1.7",
     "sqlitedict>=2.0",
     "termcolor",

src/shapepipe/modules/ngmix_package/ngmix.py

@@ -181,11 +181,10 @@ def __init__(
         self.flags = []
         self.bkg_rms = []
         self.jacobs = []
-        # Per-epoch full WCS and the object's sky position, used only by the
-        # "wcs" centroid source (skipped for the default "hsm" path).
-        self.wcs = []
-        self.ra = []
-        self.dec = []
+        # Per-epoch sub-pixel [row, col] coadd-centroid offset propagated from
+        # the stamp extractor, used only by the "wcs" centroid source (skipped
+        # for the default "hsm" path).
+        self.offsets = []
         self.bkg_sub = bkg_sub
         self.megacam_flip = megacam_flip
 
@@ -856,10 +855,14 @@ def prepare_postage_stamps(vignet, obj_id, i_tile, tile_cat):
         stamp.flags.append(flag_vign)
         stamp.bkg_rms.append(bkg_rms_vign_scaled)
         stamp.jacobs.append(jacob)
-        # For the "wcs" centroid source (see make_ngmix_observation).
-        stamp.wcs.append(epoch_wcs)
-        stamp.ra.append(tile_cat.ra[i_tile])
-        stamp.dec.append(tile_cat.dec[i_tile])
+        # Coadd-centroid offset the stamp extractor used, propagated on the
+        # galaxy vignette. Consumed only by the "wcs" centroid source (see
+        # make_ngmix_observation), which raises if it is missing; the "hsm"
+        # path ignores it, so read it leniently rather than coupling hsm to a
+        # field it never uses.
+        stamp.offsets.append(
+            vignet.gal_vign_cat[str(obj_id)][expccd_name].get('OFFSET')
+        )
 
     return stamp
 
@@ -1055,7 +1058,7 @@ def prepare_ngmix_weights(gal, weight, flag, rng, bkg_rms=None):
 
 def make_ngmix_observation(
     gal, weight, flag, psf, wcs, rng,
-    bkg_rms=None, centroid_source="hsm", wcs_full=None, ra=None, dec=None,
+    bkg_rms=None, centroid_source="hsm", offset=None,
 ):
     """Build an ngmix Observation for a single galaxy epoch.
 
@@ -1067,11 +1070,15 @@ def make_ngmix_observation(
       measured from the stamp. Robust for **galaxies**: it follows the actual
       light and so the centroid prior (centered at the Jacobian origin) does
       not bias an object that is offset from the stamp center.
-    * ``"wcs"`` — place the origin at the object's catalog sky position,
-      projected through the WCS to a sub-pixel pixel offset from the stamp
+    * ``"wcs"`` — the **coadd centroid** (Axel, #767): place the origin at the
+      object's catalog sky position as a sub-pixel offset from the stamp
       center, with no shape measurement. Better for **stars**: their HSM
       moments are noisy, so trusting the astrometry is more stable than
-      re-measuring the centroid.
+      re-measuring the centroid. The offset is not recomputed here — it is the
+      value the stamp extractor already used to round the extraction pixel
+      (:func:`shapepipe.modules.vignetmaker_package.vignetmaker.get_stamps`),
+      propagated on the vignette. One projection, one rounding: the extraction
+      and the centroid prior cannot disagree near a rounding tie.
 
     Parameters
     ----------
@@ -1088,12 +1095,10 @@ def make_ngmix_observation(
         Per-pixel background RMS map.
     centroid_source : {"hsm", "wcs"}, optional
         How to place the galaxy Jacobian origin; the default is ``"hsm"``.
-    wcs_full : astropy.wcs.WCS, optional
-        Full exposure WCS for the object's CCD. Required for
-        ``centroid_source="wcs"`` (ignored for ``"hsm"``).
-    ra, dec : float, optional
-        Object sky position in degrees. Required for
-        ``centroid_source="wcs"`` (ignored for ``"hsm"``).
+    offset : array_like, optional
+        Sub-pixel ``[row, col]`` coadd-centroid offset propagated from the
+        stamp extractor. Required for ``centroid_source="wcs"`` (ignored for
+        ``"hsm"``).
 
     Returns
     -------
@@ -1137,16 +1142,20 @@ def make_ngmix_observation(
         except Exception:
             cen_row, cen_col = 0.0, 0.0
     elif centroid_source == "wcs":
-        # Place the origin at the catalog sky position projected through the
-        # WCS — no shape measurement. Stars have noisy HSM moments, so trust
-        # the astrometry instead of re-measuring the centroid.
-        g_wcs = galsim.fitswcs.AstropyWCS(wcs=wcs_full)
-        world_pos = galsim.CelestialCoord(
-            ra * galsim.degrees, dec * galsim.degrees
-        )
-        pos = g_wcs.toImage(world_pos)
-        cen_col = pos.x - np.round(pos.x).astype(int)
-        cen_row = pos.y - np.round(pos.y).astype(int)
+        # Coadd centroid: use the sub-pixel offset the stamp extractor already
+        # computed and rounded against (propagated on the vignette), rather
+        # than re-projecting the sky position through the WCS and re-rounding.
+        # This is the fix for #767 — one projection, one rounding, so a
+        # milli-pixel WCS disagreement can no longer flip a rounding tie and
+        # put the prior a whole pixel off.
+        if offset is None:
+            raise ValueError(
+                "centroid_source='wcs' requires the coadd-centroid offset "
+                + "propagated from the stamp extractor, but none was found on "
+                + "the vignette (were the stamps made before in-house "
+                + "extraction landed?)"
+            )
+        cen_row, cen_col = float(offset[0]), float(offset[1])
     else:
         raise ValueError(
             f"Unknown centroid_source '{centroid_source}'; expected"
@@ -1386,9 +1395,7 @@ def do_ngmix_metacal(stamp, prior, flux_guess, rng, centroid_source="hsm"):
             rng,
             bkg_rms=bkg_rms,
             centroid_source=centroid_source,
-            wcs_full=stamp.wcs[n_e] if n_e < len(stamp.wcs) else None,
-            ra=stamp.ra[n_e] if n_e < len(stamp.ra) else None,
-            dec=stamp.dec[n_e] if n_e < len(stamp.dec) else None,
+            offset=stamp.offsets[n_e] if n_e < len(stamp.offsets) else None,
         )
         gal_obs_list.append(gal_obs)
 

src/shapepipe/modules/python_example_runner.py

@@ -24,7 +24,6 @@
         "pandas",
         "pysap",
         "scipy",
-        "sf_tools",
         "sqlitedict",
     ],
     run_method="parallel",

src/shapepipe/modules/vignetmaker_package/vignetmaker.py

@@ -11,12 +11,80 @@
 
 import numpy as np
 from astropy.wcs import WCS
-from sf_tools.image.stamp import FetchStamps
 from sqlitedict import SqliteDict
 
 from shapepipe.pipeline import file_io
 
 
+def get_stamps(image, positions, rad):
+    """Get Stamps.
+
+    Extract postage stamps and record their sub-pixel centring.
+
+    The image is zero-padded by ``rad`` on every side and a
+    ``(2 * rad + 1, 2 * rad + 1)`` stamp is sliced around the integer pixel
+    nearest each position (``numpy.round``). The stamp values are
+    bit-identical to the ``sf_tools.image.stamp.FetchStamps`` this replaces
+    for in-bounds positions; unlike ``FetchStamps``, an out-of-bounds centre
+    raises rather than being silently wrapped modulo the image shape.
+
+    Parameters
+    ----------
+    image : numpy.ndarray
+        2D image array to extract from
+    positions : numpy.ndarray
+        ``(N, 2)`` array of 0-indexed ``[row, col]`` float positions
+    rad : int
+        Stamp radius; stamps are ``2 * rad + 1`` pixels on a side
+
+    Returns
+    -------
+    stamps : numpy.ndarray
+        ``(N, 2 * rad + 1, 2 * rad + 1)`` array of extracted stamps
+    int_pos : numpy.ndarray
+        ``(N, 2)`` integer array of the rounded ``[row, col]`` pixel each
+        stamp is centred on
+    offsets : numpy.ndarray
+        ``(N, 2)`` float array of the sub-pixel ``[row, col]`` offset
+        (``positions - int_pos``), in ``[-0.5, 0.5]``. This is the centroid
+        prior consumed downstream by the ngmix ``"wcs"`` (coadd) centroid
+        source, so the stamp and the prior share a single rounding and cannot
+        disagree.
+
+    Raises
+    ------
+    ValueError
+        If ``positions`` is not ``(N, 2)`` or any stamp centre rounds to a
+        pixel outside the image.
+
+    """
+    rad = int(rad)
+    positions = np.asarray(positions, dtype=float)
+    if positions.ndim != 2 or positions.shape[1] != 2:
+        raise ValueError("positions must have shape (N, 2)")
+
+    shape = np.array(image.shape)
+    int_pos = np.round(positions).astype(int)
+    offsets = positions - int_pos
+
+    out_of_bounds = np.any((int_pos < 0) | (int_pos >= shape), axis=1)
+    if np.any(out_of_bounds):
+        raise ValueError(
+            f"Stamp centre(s) fall outside the image (shape {tuple(shape)}): "
+            + f"{positions[out_of_bounds].tolist()}"
+        )
+
+    padded = np.pad(image, ((rad, rad), (rad, rad)), mode="constant")
+    stamps = np.array(
+        [
+            padded[row : row + 2 * rad + 1, col : col + 2 * rad + 1]
+            for row, col in int_pos
+        ]
+    )
+
+    return stamps, int_pos, offsets
+
+
 class VignetMaker(object):
     """Vignet Maker.
 
@@ -91,7 +159,10 @@ def process(self, image_path_list, rad, prefix):
                     + "'SPHE' (spherical world coordinates))"
                 )
 
-            vign = self._get_stamp(image_path, pos - 1, rad)
+            # Single-epoch (tile) path: the coadd-centroid offset is not
+            # consumed here (only the multi-epoch path feeds the ngmix
+            # centroid prior), so the extraction bookkeeping is discarded.
+            vign, _, _ = self._get_stamp(image_path, pos - 1, rad)
 
             save_vignet(
                 vign,
@@ -186,8 +257,13 @@ def _get_stamp(self, img_path, pos, rad):
 
         Returns
         -------
-        numpy.ndarray
+        stamps : numpy.ndarray
             Array containing the vignets
+        int_pos : numpy.ndarray
+            ``(N, 2)`` rounded ``[row, col]`` extraction pixels
+        offsets : numpy.ndarray
+            ``(N, 2)`` sub-pixel ``[row, col]`` offsets (see
+            :func:`get_stamps`)
 
         """
         try:
@@ -211,12 +287,7 @@ def _get_stamp(self, img_path, pos, rad):
             raise
         img_file.close()
 
-        fs = FetchStamps(img, int(rad))
-        fs.get_pixels(np.round(pos).astype(int))
-
-        vign = fs.scan()
-
-        return vign
+        return get_stamps(img, pos, int(rad))
 
     def _get_stamp_me(self, image_dirs, image_pattern):
         """Get Stamp Multi-Epoch.
@@ -266,6 +337,8 @@ def _get_stamp_me(self, image_dirs, image_pattern):
             array_vign = None
             array_id = None
             array_exp_name = None
+            array_offset = None
+            array_int_pos = None
 
             for ccd in ccd_list:
 
@@ -298,13 +371,29 @@ def _get_stamp_me(self, image_dirs, image_pattern):
                 ).T
                 pos[:, [0, 1]] = pos[:, [1, 0]]
 
-                tmp_vign = self._get_stamp(img_path, pos - 1, self._rad)
+                tmp_vign, tmp_int_pos, tmp_offset = self._get_stamp(
+                    img_path, pos - 1, self._rad
+                )
 
                 if array_vign is None:
                     array_vign = np.copy(tmp_vign)
                 else:
                     array_vign = np.concatenate((array_vign, np.copy(tmp_vign)))
 
+                if array_offset is None:
+                    array_offset = np.copy(tmp_offset)
+                else:
+                    array_offset = np.concatenate(
+                        (array_offset, np.copy(tmp_offset))
+                    )
+
+                if array_int_pos is None:
+                    array_int_pos = np.copy(tmp_int_pos)
+                else:
+                    array_int_pos = np.concatenate(
+                        (array_int_pos, np.copy(tmp_int_pos))
+                    )
+
                 if array_id is None:
                     array_id = np.copy(obj_id)
                 else:
@@ -319,7 +408,15 @@ def _get_stamp_me(self, image_dirs, image_pattern):
                     array_exp_name = np.concatenate((array_exp_name, exp_name_tmp))
 
             if array_id is not None:
-                final_list.append([array_id, array_vign, array_exp_name])
+                final_list.append(
+                    [
+                        array_id,
+                        array_vign,
+                        array_exp_name,
+                        array_offset,
+                        array_int_pos,
+                    ]
+                )
             else:
                 self._w_log.info(
                     f"All CCDs skipped for exposure {exp_name} "
@@ -344,6 +441,20 @@ def _get_stamp_me(self, image_dirs, image_pattern):
                     output_dict[id_tmp][index]["VIGNET"] = final_list[j][1][
                         where_res[0]
                     ]
+                    # Sub-pixel [row, col] offset of the coadd centroid from
+                    # the integer pixel this stamp was extracted around, plus
+                    # that integer pixel itself. Both ride with the vignette
+                    # (no sidecar file); OFFSET is consumed by the ngmix "wcs"
+                    # (coadd) centroid source so the extraction and the
+                    # centroid prior share one rounding (see #767), and
+                    # INT_POS + OFFSET reconstruct the absolute exposure
+                    # centroid for diagnostics.
+                    output_dict[id_tmp][index]["OFFSET"] = final_list[j][3][
+                        where_res[0]
+                    ]
+                    output_dict[id_tmp][index]["INT_POS"] = final_list[j][4][
+                        where_res[0]
+                    ]
                     counter += 1
 
             if counter == 0:

src/shapepipe/modules/vignetmaker_runner.py

@@ -18,7 +18,7 @@
     input_module="sextractor_runner",
     file_pattern=["galaxy_selection", "image"],
     file_ext=[".fits", ".fits"],
-    depends=["numpy", "astropy", "sf_tools", "sqlitedict"],
+    depends=["numpy", "astropy", "sqlitedict"],
 )
 def vignetmaker_runner(
     input_file_list,