psf validation (with 4th-order moments) enabled: merge star cat, convert pix

Author: martinkilbinger

docs/source/pipeline_canfar.md

@@ -122,12 +122,17 @@ create_final_cat.py -m final_cat_$patch.hdf5 -i . -p $patch/cfis/final_cat.param
 
 # Star catalogue
 combine_runs.bash  -p $psf -c psf
+
+## The following are not necessary, only to get merged catalogue with non-converted quantities
 shapepipe_run -c $SP_CONFIG/config_Ms_$psf.ini
 shapepipe_run -c $SP_CONFIG/config_Pl_$psf.ini
 
 # With fourth-order moment, additional psf interpol validation runs
 combine_runs.bash  -p $psf -c psf_val
 
+## The following are not necessary, only to get merged catalogue with non-converted quantities
+shapepipe_run -c $SP_CONFIG/config_Msval_$psf.ini
+
 
 # Convert star cat to WCS
 ## Convert all input validation psf files and create directories par patch

pyproject.toml

@@ -58,9 +58,10 @@ include-package-data = true
     "scripts/sh/curl_canfar_local.sh",
     "scripts/sh/functions.sh",
     "scripts/sh/stats_jobs_canfar.sh",
-    "scripts/sh/combine_runs.sh",
+    "scripts/sh/combine_runs.bash",
     "scripts/python/update_runs_log_file.py",
     "scripts/python/summary_run.py",
     "scripts/python/summary_params_pre_v2.py",
     "scripts/python/link_to_exp_for_tile.py",
+    "scripts/python/convert_psf_pix2world.py",
 ]

scripts/python/convert_psf_pix2world.py

@@ -16,45 +16,94 @@
 from cs_util import logging
 
 
-def transform_shape(mom_list, jac):
+def transform_shape(mom_list, jac, rescale=False):
     """Transform Shape.
 
     Transform shape (ellipticity and size) using a Jacobian.
 
     Parameters
     ----------
     mom_list : list
-        input moment measurements; each list element contains
-        first and second ellipticity component and size
+        input "moment" measurements; each list element contains
+        first and second ellipticity components and size
     jac : galsim.JacobianWCS
         Jacobian transformation matrix information
+    rescale : bool, optional
+        rescale ellipticity to modulus if ``True``, to prevent
+        values outside [-1; 1]; default is ``False``
 
     Returns
     -------
     list
         transformed shape parameters, which are
-        first and second ellipticity component and size
+        first and second ellipticity components and size
 
     """
     scale, shear, theta, flip = jac.getDecomposition()
 
     sig_tmp = mom_list[2] * scale
-    shape = galsim.Shear(g1=mom_list[0], g2=mom_list[1])
+
+    if rescale:
+        modulus = np.sqrt(mom_list[0] ** 2 + mom_list[1] ** 2) * 1.1
+    else:
+        modulus = 1
+
+    shape = galsim.Shear(g1=mom_list[0]/modulus, g2=mom_list[1]/modulus)
+
+    if flip:
+        print(f"g1 sign change, flip = {flip}")
+        shape = galsim.Shear(g1=-shape.g1, g2=shape.g2)
+
+    shape = galsim.Shear(g=shape.g, beta=shape.beta + theta)
+    shape = shear + shape
+
+    return shape.g1 * modulus, shape.g2 * modulus, sig_tmp
+
+
+def transform_shape_scale(mom_list, jac):
+    """Transform Shape.
+
+    Transform shape (ellipticity and size) using a Jacobian.
+    Input "ellipticities" can be outside [-1, 1].
+
+    Parameters
+    ----------
+    mom_list : list
+        input "moment" measurements; each list element contains
+        first and second ellipticity components
+    jac : galsim.JacobianWCS
+        Jacobian transformation matrix information
+
+    Returns
+    -------
+    list
+        transformed shape parameters, which are
+        first and second "ellipticity" components
+
+    """
+    _, shear, theta, flip = jac.getDecomposition()
+
+    # Scale to within [-1, 1]
+    modulus = np.sqrt(mom_list[0] ** 2 + mom_list[1] ** 2) * 1.1
+    shape = galsim.Shear(g1=mom_list[0]/modulus, g2=mom_list[1]/modulus)
+
     if flip:
-        # The following output is not observed
-        print("FLIP!")
+        print(f"g1 sign change, flip = {flip}")
         shape = galsim.Shear(g1=-shape.g1, g2=shape.g2)
+
     shape = galsim.Shear(g=shape.g, beta=shape.beta + theta)
     shape = shear + shape
 
-    return shape.g1, shape.g2, sig_tmp
+    # Scale back to original modulus
+
+    return shape.g1 * modulus, shape.g2 * modulus
 
 
 class Loc2Glob(object):
     r"""Change from local to global coordinates.
 
     Class to pass from local coordinates to global coordinates under
-    CFIS (CFHT) MegaCam instrument. The geometrical informcation of the
+    CFIS (CFHT) MegaCam instrument. The geometrical information of the
     instrument is encoded in this function.
 
     Parameters
@@ -378,6 +427,7 @@ def params_default(self):
         self._params = {
             "input_base_dir": ".",
             "output_base_dir": ".",
+            "sub_dir_pattern": "run_sp_combined_psf",
             "mode": "merge",
             "patches": "",
             "psf": "psfex",
@@ -399,6 +449,7 @@ def params_default(self):
                 + " <input_base_dir>/P<patch?>/output;"
                 " default is {}"
             ),
+            "sub_dir_pattern": "run directory pattern, default is {}",
             "mode": (
                 "run mode, allowed are 'merge', 'test'; default is" + " '{}'"
             ),
@@ -415,10 +466,14 @@ def params_default(self):
             ("E1_PSF_HSM", float),
             ("E2_PSF_HSM", float),
             ("SIGMA_PSF_HSM", float),
+            ("M_4_PSF_1", float),
+            ("M_4_PSF_2", float),
             ("FLAG_PSF_HSM", float),
             ("E1_STAR_HSM", float),
             ("E2_STAR_HSM", float),
             ("SIGMA_STAR_HSM", float),
+            ("M_4_STAR_1", float),
+            ("M_4_STAR_2", float),
             ("FLAG_STAR_HSM", float),
             ("CCD_NB", int),
         ]
@@ -435,11 +490,8 @@ def update_params(self):
 
         """
         if self._params["psf"] == "psfex":
-            #self._params["sub_dir_pattern"] = "run_sp_exp_202"
-            self._params["sub_dir_pattern"] = "run_sp_combined_psf"
             self._params["sub_dir_psfint"] = "psfex_interp_runner"
         elif self._params["psf"] == "mccd":
-            self._params["sub_dir_pattern"] = "run_sp_exp_SxSePsf_202"
             self._params["sub_dir_psfint"] = "mccd_fit_val_runner"
             self._params["sub_dir_setools"] = "setools_runner/output/mask"
         else:
@@ -471,7 +523,7 @@ def run(self):
             print("Running patch:", patch)
 
             patch_dir = f"{self._params['input_base_dir']}/P{patch}/output/"
-            subdirs = f"{patch_dir}/{self._params['sub_dir_pattern']}*"
+            subdirs = f"{patch_dir}/{self._params['sub_dir_pattern']}"
             exp_run_dirs = glob.glob(subdirs)
             n_exp_runs = len(exp_run_dirs)
             print(
@@ -567,9 +619,13 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
             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_m41_psf = np.zeros_like(psf_file[mod])
+            new_m42_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_m41_star = np.zeros_like(psf_file[mod])
+            new_m42_star = np.zeros_like(psf_file[mod])
 
             if self._params["psf"] == "psfex":
                 header = fits.Header.fromstring(
@@ -580,10 +636,6 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
                 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"],
@@ -592,6 +644,7 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
                         )
                     except Exception:
                         continue
+
                     g1_psf_tmp, g2_psf_tmp, sig_psf_tmp = transform_shape(
                         [
                             obj["E1_PSF_HSM"],
@@ -600,10 +653,14 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
                         ],
                         jac,
                     )
-
                     new_e1_psf[ind] = g1_psf_tmp
                     new_e2_psf[ind] = g2_psf_tmp
                     new_sig_psf[ind] = sig_psf_tmp
+                    new_m41_psf[ind], new_m42_psf[ind], _ = transform_shape(
+                        [obj["M_4_PSF_1"], obj["M_4_PSF_2"], 0],
+                        jac,
+                        rescale=True,
+                    )
 
                     g1_star_tmp, g2_star_tmp, sig_star_tmp = transform_shape(
                         [
@@ -616,6 +673,12 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
                     new_e1_star[ind] = g1_star_tmp
                     new_e2_star[ind] = g2_star_tmp
                     new_sig_star[ind] = sig_star_tmp
+                    new_m41_star[ind], new_m42_star[ind], _ = transform_shape(
+                        [obj["M_4_STAR_1"], obj["M_4_STAR_2"], 0],
+                        jac,
+                        rescale=True,
+                    )
+
                     k += 1
 
                 exp_cat = np.array(
@@ -631,10 +694,14 @@ def transform_exposures(self, output_dir, patch, idx, exp_run_dir):
                                     new_e1_psf,
                                     new_e2_psf,
                                     new_sig_psf,
+                                    new_m41_psf,
+                                    new_m42_psf,
                                     psf_file["FLAG_PSF_HSM"],
                                     new_e1_star,
                                     new_e2_star,
                                     new_sig_star,
+                                    new_m41_star,
+                                    new_m42_star,
                                     psf_file["FLAG_STAR_HSM"],
                                     np.ones_like(psf_file["RA"], dtype=int)
                                     * ccd_id,

src/shapepipe/modules/merge_starcat_package/merge_starcat.py

@@ -554,8 +554,8 @@ def process(self):
 
         """
         x, y, ra, dec = [], [], [], []
-        g1_psf, g2_psf, size_psf = [], [], []
-        g1, g2, size = [], [], []
+        g1_psf, g2_psf, size_psf, m41_psf, m42_psf = [], [], [], [], []
+        g1, g2, size, m41, m42 = [], [], [], [], []
         flag_psf, flag_star = [], []
         mag, snr, psfex_acc = [], [], []
         ccd_nb = []
@@ -584,8 +584,12 @@ def process(self):
             g1_psf += list(data_j["E1_PSF_HSM"])
             g2_psf += list(data_j["E2_PSF_HSM"])
             size_psf += list(data_j["SIGMA_PSF_HSM"] ** 2)
+            m41_psf += list(data_j["M_4_PSF_1"])
+            m42_psf += list(data_j["M_4_PSF_1"])
             g1 += list(data_j["E1_STAR_HSM"])
             g2 += list(data_j["E2_STAR_HSM"])
+            m41 += list(data_j["M_4_STAR_1"])
+            m42 += list(data_j["M_4_STAR_2"])
             size += list(data_j["SIGMA_STAR_HSM"] ** 2)
 
             # flags
@@ -633,9 +637,13 @@ def process(self):
             "E1_PSF_HSM": g1_psf,
             "E2_PSF_HSM": g2_psf,
             "SIGMA_PSF_HSM": np.sqrt(size_psf),
+            "M_4_PSF_1": m41_psf,
+            "M_4_PSF_2": m42_psf,
             "E1_STAR_HSM": g1,
             "E2_STAR_HSM": g2,
             "SIGMA_STAR_HSM": np.sqrt(size),
+            "M_4_STAR_1": m41,
+            "M_4_STAR_2": m42,
             "FLAG_PSF_HSM": flag_psf,
             "FLAG_STAR_HSM": flag_star,
             "MAG": mag,