HSM: measure PSF/star shapes directly in sky coords (use_sky_coords)

docs/source/pipeline_canfar.md

@@ -333,21 +333,43 @@ 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
 
-Convert all input validation PSF files and create directories per patch `P?`.
-Create files `validation_psf_conv-<patchnum>-<idx>.fits` (for the v1.4 setup only one file):
+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
+`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<n>`, 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-<idx>.fits` instead of
+> `validation_psf_conv-<patchnum>-<idx>.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-<idx>.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-<patchnum>-<idx>.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).

scripts/python/convert_psf_pix2world.py → scripts/python/collate_star_cat.py

@@ -1,5 +1,32 @@
 #! /usr/bin/env python3
 
+"""COLLATE STAR CATALOGUES.
+
+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<n>``): input runs live under
+``<base>/P<patch>/output`` and outputs are named
+``validation_psf_conv-<patch>-<idx>.fits``. ``v2.0`` removes the patch concept:
+input runs live under a single ``<base>/output`` root and outputs drop the patch
+token (``validation_psf_conv-<idx>.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 (``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
 import os
 import re
@@ -16,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<patch>`` 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.
 
@@ -378,6 +457,7 @@ def params_default(self):
         self._params = {
             "input_base_dir": ".",
             "output_base_dir": ".",
+            "version_cat": "v2.0",
             "mode": "merge",
             "patches": "",
             "psf": "psfex",
@@ -386,6 +466,7 @@ def params_default(self):
 
         self._short_options = {
             "input_base_dir": "-i",
+            "version_cat": "-V",
             "mode": "-m",
             "psf": "-p",
             "patches": "-P",
@@ -395,15 +476,19 @@ def params_default(self):
 
         self._help_strings = {
             "input_base_dir": (
-                "input base dir, runs are expected in"
-                + " <input_base_dir>/P<patch?>/output;"
-                " default is {}"
+                "input base dir; for v1.x runs are expected in"
+                + " <input_base_dir>/P<patch>/output, for v2.0 (patch-less) in"
+                + " <input_base_dir>/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
@@ -454,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"])
@@ -464,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)
@@ -516,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")
@@ -555,7 +660,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 +668,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 +682,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,8 +715,23 @@ 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):