PSF fourth-order moments in the sky frame (spin-2 + rho4)

src/shapepipe/modules/make_cat_package/make_cat.py

@@ -619,6 +619,12 @@ def _save_psf_data(self, galaxy_psf_path):
 
         Save the PSF catalogue into the final one.
 
+        @sc [label:schema] psf-epoch-slot-columns
+        The per-epoch families are ``psf_shape_cols`` plus ``EXP_ID``/``CCD``,
+        all slot-aligned by the producer's epoch enumeration; the names must
+        be a subset of what ``PSFExInterpolator._interpolate_me`` writes into
+        ``SHAPES`` (``test_hsm_column_seams``).
+
         Parameters
         ----------
         galaxy_psf_path : str
@@ -628,49 +634,30 @@ def _save_psf_data(self, galaxy_psf_path):
         galaxy_psf_cat = SqliteDict(galaxy_psf_path)
 
         max_epoch = np.max(self._final_cat_file.get_data()["N_EPOCH"]) + 1
+        n_obj = len(self._obj_id)
 
-        self._output_dict = {
-            f"HSM_G1_PSF_{idx + 1}": np.ones(len(self._obj_id)) * -10.0
-            for idx in range(max_epoch)
-        }
-        self._output_dict = {
-            **self._output_dict,
-            **{
-                f"HSM_G2_PSF_{idx + 1}": np.ones(len(self._obj_id)) * -10.0
-                for idx in range(max_epoch)
-            },
-        }
-        self._output_dict = {
-            **self._output_dict,
-            **{
-                f"HSM_T_PSF_{idx + 1}": np.zeros(len(self._obj_id))
-                for idx in range(max_epoch)
-            },
-        }
-        self._output_dict = {
-            **self._output_dict,
-            **{
-                f"HSM_FLAG_PSF_{idx + 1}": np.ones(
-                    len(self._obj_id), dtype="int16"
-                )
-                for idx in range(max_epoch)
-            },
-        }
+        # Per-epoch PSF shape columns copied from the producer's SHAPES dict:
+        # (column, empty-slot fill, dtype). Fills are out of physical range
+        # so an unmeasured slot cannot pass for a measurement. HSM_T_PSF
+        # already holds T (sigma_to_T applied at the producer's
+        # _interpolate_me).
+        psf_shape_cols = [
+            ("HSM_G1_PSF", -10.0, float),
+            ("HSM_G2_PSF", -10.0, float),
+            ("HSM_T_PSF", 0.0, float),
+            ("HSM_FLAG_PSF", 1, "int16"),
+            ("HSM_M4_1_PSF", -10.0, float),
+            ("HSM_M4_2_PSF", -10.0, float),
+            ("HSM_RHO4_PSF", -1.0, float),
+        ]
         # Per-epoch exposure ID and CCD number, slot-aligned with HSM_*_PSF_n;
         # -1 marks an empty slot (the CCD_N sentinel convention).
+        epoch_id_cols = [("EXP_ID", -1, "int32"), ("CCD", -1, "int32")]
+
         self._output_dict = {
-            **self._output_dict,
-            **{
-                f"EXP_ID_{idx + 1}": np.ones(len(self._obj_id), dtype="int32") * -1
-                for idx in range(max_epoch)
-            },
-        }
-        self._output_dict = {
-            **self._output_dict,
-            **{
-                f"CCD_{idx + 1}": np.ones(len(self._obj_id), dtype="int32") * -1
-                for idx in range(max_epoch)
-            },
+            f"{name}_{idx + 1}": np.full(n_obj, fill, dtype=dtype)
+            for name, fill, dtype in psf_shape_cols + epoch_id_cols
+            for idx in range(max_epoch)
         }
 
         for idx, id_tmp in enumerate(self._obj_id):
@@ -680,7 +667,7 @@ def _save_psf_data(self, galaxy_psf_path):
 
             for epoch, key in enumerate(galaxy_psf_cat[str(id_tmp)].keys()):
 
-                gpc_data = galaxy_psf_cat[str(id_tmp)][key]
+                shapes = galaxy_psf_cat[str(id_tmp)][key]["SHAPES"]
 
                 # `key` is "<exp>-<ccd>"; reading it in the enumeration that
                 # assigns `epoch` aligns EXP_ID_n/CCD_n with HSM_*_PSF_n by
@@ -690,28 +677,15 @@ def _save_psf_data(self, galaxy_psf_path):
                 self._add2dict(f"EXP_ID_{epoch + 1}", int(exp_name), idx)
                 self._add2dict(f"CCD_{epoch + 1}", int(ccd_n), idx)
 
-                if gpc_data["SHAPES"]["HSM_FLAG_PSF"] != 0:
+                if shapes["HSM_FLAG_PSF"] != 0:
                     continue
 
-                self._add2dict(
-                    f"HSM_G1_PSF_{epoch + 1}",
-                    gpc_data["SHAPES"]["HSM_G1_PSF"], idx
-                )
-                self._add2dict(
-                    f"HSM_G2_PSF_{epoch + 1}",
-                    gpc_data["SHAPES"]["HSM_G2_PSF"], idx
-                )
-
-                # HSM_T_PSF already holds T (sigma_to_T applied at the
-                # producer's _interpolate_me); read straight through.
-                self._add2dict(
-                    f"HSM_T_PSF_{epoch + 1}",
-                    gpc_data["SHAPES"]["HSM_T_PSF"], idx
-                )
-
-                self._add2dict(
-                    f"HSM_FLAG_PSF_{epoch + 1}",
-                    gpc_data["SHAPES"]["HSM_FLAG_PSF"], idx
-                )
+                for name, fill, _ in psf_shape_cols:
+                    # A SHAPES dict without the fourth-moment keys (MCCD, or
+                    # a producer predating them) leaves those slots at their
+                    # out-of-range fill while FLAG still reports the fit.
+                    self._add2dict(
+                        f"{name}_{epoch + 1}", shapes.get(name, fill), idx
+                    )
 
         galaxy_psf_cat.close()

src/shapepipe/modules/merge_starcat_package/merge_starcat.py

@@ -556,10 +556,17 @@ def process(self):
 
         Process merging.
 
+        @sc [label:schema] psfex-starcat-columns-strict
+        Every ``HSM_*`` column is read by name with no fallback, so the set
+        read here equals the set ``PSFExInterpolator._write_output_validation``
+        writes (``test_hsm_column_seams``).
+
         """
         x, y, ra, dec = [], [], [], []
         g1_psf, g2_psf, size_psf = [], [], []
+        m4_1_psf, m4_2_psf, rho4_psf = [], [], []
         g1, g2, size = [], [], []
+        m4_1_star, m4_2_star, rho4_star = [], [], []
         flag_psf, flag_star = [], []
         mag, snr, psfex_acc = [], [], []
         ccd_nb = []
@@ -588,9 +595,15 @@ def process(self):
             g1_psf += list(data_j["HSM_G1_PSF"])
             g2_psf += list(data_j["HSM_G2_PSF"])
             size_psf += list(data_j["HSM_T_PSF"])
+            m4_1_psf += list(data_j["HSM_M4_1_PSF"])
+            m4_2_psf += list(data_j["HSM_M4_2_PSF"])
+            rho4_psf += list(data_j["HSM_RHO4_PSF"])
             g1 += list(data_j["HSM_G1_STAR"])
             g2 += list(data_j["HSM_G2_STAR"])
             size += list(data_j["HSM_T_STAR"])
+            m4_1_star += list(data_j["HSM_M4_1_STAR"])
+            m4_2_star += list(data_j["HSM_M4_2_STAR"])
+            rho4_star += list(data_j["HSM_RHO4_STAR"])
 
             # flags
             flag_psf += list(data_j["HSM_FLAG_PSF"])
@@ -636,9 +649,15 @@ def process(self):
             "HSM_G1_PSF": g1_psf,
             "HSM_G2_PSF": g2_psf,
             "HSM_T_PSF": size_psf,
+            "HSM_M4_1_PSF": m4_1_psf,
+            "HSM_M4_2_PSF": m4_2_psf,
+            "HSM_RHO4_PSF": rho4_psf,
             "HSM_G1_STAR": g1,
             "HSM_G2_STAR": g2,
             "HSM_T_STAR": size,
+            "HSM_M4_1_STAR": m4_1_star,
+            "HSM_M4_2_STAR": m4_2_star,
+            "HSM_RHO4_STAR": rho4_star,
             "HSM_FLAG_PSF": flag_psf,
             "HSM_FLAG_STAR": flag_star,
             "MAG": mag,