diff --git a/docs/source/pipeline_canfar.md b/docs/source/pipeline_canfar.md index 7e8bb29c5..96add5758 100644 --- a/docs/source/pipeline_canfar.md +++ b/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--.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`, 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-.fits` instead of +> `validation_psf_conv--.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-.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--.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). diff --git a/scripts/python/convert_psf_pix2world.py b/scripts/python/collate_star_cat.py similarity index 79% rename from scripts/python/convert_psf_pix2world.py rename to scripts/python/collate_star_cat.py index dea99656b..f7ab25b41 100755 --- a/scripts/python/convert_psf_pix2world.py +++ b/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``): input runs live under +``/P/output`` and outputs are named +``validation_psf_conv--.fits``. ``v2.0`` removes the patch concept: +input runs live under a single ``/output`` root and outputs drop the patch +token (``validation_psf_conv-.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`` 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" - + " /P/output;" - " default is {}" + "input base dir; for v1.x runs are expected in" + + " /P/output, for v2.0 (patch-less) in" + + " /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): diff --git a/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py b/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py index bfd9f3417..a01ca095d 100644 --- a/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py +++ b/src/shapepipe/modules/mccd_package/mccd_interpolation_script.py @@ -14,6 +14,7 @@ from sqlitedict import SqliteDict from shapepipe.pipeline import file_io +from shapepipe.modules.psfex_interp_package.psfex_interp import local_wcs_list try: import galsim.hsm as hsm @@ -204,18 +205,34 @@ def _get_galaxy_positions(self): raise KeyError(pos_param_err) galcat.close() - def _get_psfshapes(self): + def _get_psfshapes(self, wcs_list=None): """Get PSF Shapes. Compute shapes of PSF at galaxy positions using HSM. + Parameters + ---------- + wcs_list : list of galsim.BaseWCS, optional + Local WCS at each object position. When provided, HSM adaptive + moments are measured in world coordinates (``use_sky_coords=True``), + so the returned shapes and sizes are already in sky coordinates. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") - psf_moms = [ - hsm.FindAdaptiveMom(Image(psf), strict=False) for psf in self.interp_PSFs - ] + if wcs_list is not None: + psf_moms = [ + hsm.FindAdaptiveMom( + Image(psf, wcs=wcs), strict=False, use_sky_coords=True + ) + for psf, wcs in zip(self.interp_PSFs, wcs_list) + ] + else: + psf_moms = [ + hsm.FindAdaptiveMom(Image(psf), strict=False) + for psf in self.interp_PSFs + ] self.psf_shapes = np.array( [ @@ -253,7 +270,7 @@ def _write_output(self): data = {"VIGNET": self.interp_PSFs} output.save_as_fits(data, sex_cat_path=self._galcat_path) - def _get_starshapes(self, star_vign): + def _get_starshapes(self, star_vign, wcs_list=None): """Get Star Shapes. Compute shapes of stars at stars positions using HSM. @@ -262,6 +279,10 @@ def _get_starshapes(self, star_vign): ---------- star_vign : numpy.ndarray Array containing the star's vignets + wcs_list : list of galsim.BaseWCS, optional + Local WCS at each star position; see :meth:`_get_psfshapes`. When + provided, moments are measured in world coordinates + (``use_sky_coords=True``). """ if import_fail: @@ -270,10 +291,23 @@ def _get_starshapes(self, star_vign): masks = np.zeros_like(star_vign) masks[np.where(star_vign == -1e30)] = 1 - star_moms = [ - hsm.FindAdaptiveMom(Image(star), badpix=Image(mask), strict=False) - for star, mask in zip(star_vign, masks) - ] + if wcs_list is not None: + star_moms = [ + hsm.FindAdaptiveMom( + Image(star, wcs=wcs), + badpix=Image(mask), + strict=False, + use_sky_coords=True, + ) + for star, mask, wcs in zip(star_vign, masks, wcs_list) + ] + else: + star_moms = [ + hsm.FindAdaptiveMom( + Image(star), badpix=Image(mask), strict=False + ) + for star, mask in zip(star_vign, masks) + ] self.star_shapes = np.array( [ @@ -451,6 +485,15 @@ def _interpolate_me(self): 0, ) ).T + # 1-indexed (FITS) positions for local_wcs_list, which wraps + # galsim.AstropyWCS; gal_pos (origin 0) stays for MCCD interp. + gal_pos_wcs = np.array( + self._f_wcs_file[exp_name][ccd]["WCS"].all_world2pix( + self.gal_pos[:, 0][ind_obj], + self.gal_pos[:, 1][ind_obj], + 1, + ) + ).T self.interp_PSFs = interp_MCCD(mccd_model_path, gal_pos, ccd) @@ -497,7 +540,11 @@ def _interpolate_me(self): array_id = np.concatenate((array_id, np.copy(obj_id))) if self._compute_shape: - self._get_psfshapes() + self._get_psfshapes( + local_wcs_list( + self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos_wcs + ) + ) if array_shape is None: array_shape = np.copy(self.psf_shapes) else: diff --git a/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py b/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py index 943433e8f..f0aea4f65 100644 --- a/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py +++ b/src/shapepipe/modules/mccd_package/shapepipe_auxiliary_mccd.py @@ -16,6 +16,7 @@ from astropy.io import fits from cs_util import size as cs_size +from shapepipe.modules.psfex_interp_package.psfex_interp import local_wcs_list from shapepipe.pipeline import file_io NOT_ENOUGH_STARS = "Not enough stars to train the model." @@ -403,6 +404,10 @@ def mccd_interpolation_pipeline( x_pos = galcat[2].data[pos_params[0]] y_pos = galcat[2].data[pos_params[1]] interp_pos = np.array([x_pos, y_pos]).T + # CCD image WCS, stored as a card list in the FITS_LDAC LDAC_IMHEAD HDU + wcs_header = fits.Header.fromstring( + "\n".join(galcat[1].data[0][0]), sep="\n" + ) # Close catalog galcat.close() @@ -411,9 +416,18 @@ def mccd_interpolation_pipeline( if interp_PSFs is not None: if get_shapes: + # Local WCS at each object position lets HSM measure adaptive + # moments directly in world coordinates (use_sky_coords=True). + wcs_list = local_wcs_list( + galsim.AstropyWCS(header=wcs_header), interp_pos + ) PSF_moms = [ - galsim.hsm.FindAdaptiveMom(galsim.Image(psf), strict=False) - for psf in interp_PSFs + galsim.hsm.FindAdaptiveMom( + galsim.Image(psf, wcs=wcs), + strict=False, + use_sky_coords=True, + ) + for psf, wcs in zip(interp_PSFs, wcs_list) ] PSF_shapes = np.array( diff --git a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py index c14aeb136..b6da6623f 100644 --- a/src/shapepipe/modules/psfex_interp_package/psfex_interp.py +++ b/src/shapepipe/modules/psfex_interp_package/psfex_interp.py @@ -17,6 +17,7 @@ from shapepipe.pipeline import file_io try: + import galsim import galsim.hsm as hsm from galsim import Image except ImportError: @@ -30,6 +31,37 @@ FILE_NOT_FOUND = "File_not_found" +def local_wcs_list(wcs, positions): + """Local WCS List. + + Build a list of local (Jacobian) WCS from an astropy/galsim WCS, one per + object position. Attaching the local WCS to an object-centred stamp lets + ``FindAdaptiveMom(use_sky_coords=True)`` measure moments directly in world + coordinates: at the stamp's ``true_center`` galsim evaluates this local WCS, + which must therefore be the linearisation of the full WCS *at the object's + position in the original image*, not at the stamp centre. + + Parameters + ---------- + wcs : galsim.BaseWCS or astropy.wcs.WCS + Full (celestial) WCS for the CCD. + positions : numpy.ndarray + Object positions ``(N, 2)`` in image (pixel) coordinates. + + Returns + ------- + list of galsim.JacobianWCS + One local WCS per object position. + + """ + if not isinstance(wcs, galsim.BaseWCS): + wcs = galsim.AstropyWCS(wcs=wcs) + return [ + wcs.local(image_pos=galsim.PositionD(float(x), float(y))) + for x, y in positions + ] + + class PSFExInterpolator(object): """The PSFEx Interpolator Class. @@ -145,7 +177,9 @@ def process(self): self._w_log.info(f"Psf model file {self._dotpsf_path} not found.") else: if self._compute_shape: - self._get_psfshapes() + self._get_psfshapes( + local_wcs_list(self._galcat_wcs(), self.gal_pos) + ) self._write_output() @@ -307,18 +341,36 @@ def _interpolate(self): self.gal_pos, ) - def _get_psfshapes(self): + def _get_psfshapes(self, wcs_list=None): """Get PSF Shapes. Compute shapes of PSF at galaxy positions using HSM. + Parameters + ---------- + wcs_list : list of galsim.BaseWCS, optional + Local WCS at each object position (see :func:`local_wcs_list`). When + provided, HSM adaptive moments are measured directly in world + coordinates (``use_sky_coords=True``), so the ``g1``/``g2``/``sigma`` + returned are already in sky coordinates and need no post-hoc + rotation. When ``None``, moments are measured in the pixel frame. + """ if import_fail: raise ImportError("Galsim is required to get shapes information") - psf_moms = [ - hsm.FindAdaptiveMom(Image(psf), strict=False) for psf in self.interp_PSFs - ] + if wcs_list is not None: + psf_moms = [ + hsm.FindAdaptiveMom( + Image(psf, wcs=wcs), strict=False, use_sky_coords=True + ) + for psf, wcs in zip(self.interp_PSFs, wcs_list) + ] + else: + psf_moms = [ + hsm.FindAdaptiveMom(Image(psf), strict=False) + for psf in self.interp_PSFs + ] self.psf_shapes = np.array( [ @@ -407,21 +459,30 @@ def process_validation(self, psfex_cat_path): star_dict["SNR"] = np.copy(star_cat.get_data()["SNR_WIN"]) star_cat.close() - self._get_psfshapes() - self._get_starshapes(star_vign) + wcs_list = local_wcs_list( + self._galcat_wcs(), + np.column_stack((star_dict["X"], star_dict["Y"])), + ) + + self._get_psfshapes(wcs_list) + self._get_starshapes(star_vign, wcs_list) psfex_cat_dict = self._get_psfexcatdict(psfex_cat_path) self._write_output_validation(star_dict, psfex_cat_dict) - def _get_starshapes(self, star_vign): + def _get_starshapes(self, star_vign, wcs_list=None): """Get Star Shapes. Compute shapes of stars at stars positions using HSM. Parameters ---------- - numpy.ndarray + star_vign : numpy.ndarray Array containing the star's vignets. + wcs_list : list of galsim.BaseWCS, optional + Local WCS at each star position; see :meth:`_get_psfshapes`. When + provided, moments are measured in world coordinates + (``use_sky_coords=True``). """ if import_fail: @@ -430,10 +491,23 @@ def _get_starshapes(self, star_vign): masks = np.zeros_like(star_vign) masks[np.where(star_vign == -1e30)] = 1 - star_moms = [ - hsm.FindAdaptiveMom(Image(star), badpix=Image(mask), strict=False) - for star, mask in zip(star_vign, masks) - ] + if wcs_list is not None: + star_moms = [ + hsm.FindAdaptiveMom( + Image(star, wcs=wcs), + badpix=Image(mask), + strict=False, + use_sky_coords=True, + ) + for star, mask, wcs in zip(star_vign, masks, wcs_list) + ] + else: + star_moms = [ + hsm.FindAdaptiveMom( + Image(star), badpix=Image(mask), strict=False + ) + for star, mask in zip(star_vign, masks) + ] self.star_shapes = np.array( [ @@ -447,6 +521,26 @@ def _get_starshapes(self, star_vign): ] ) + def _galcat_wcs(self): + """Galaxy Catalogue WCS. + + Read the CCD image WCS from the SExtractor FITS_LDAC galaxy catalogue. + The original image header is stored as a card list in the ``LDAC_IMHEAD`` + HDU (HDU 1); this is the same header the (now retired) pixel-to-world + conversion reconstructed to rotate shapes. + + Returns + ------- + galsim.AstropyWCS + WCS of the CCD the stars/PSFs live on. + + """ + with fits.open(self._galcat_path, memmap=False) as hdu_list: + header = fits.Header.fromstring( + "\n".join(hdu_list[1].data[0][0]), sep="\n" + ) + return galsim.AstropyWCS(header=header) + def _get_psfexcatdict(self, psfex_cat_path): """Get PSFEx Catalogue Dictionary. @@ -656,6 +750,15 @@ def _interpolate_me(self): 0, ) ).T + # 1-indexed (FITS) positions for local_wcs_list, which wraps + # galsim.AstropyWCS; gal_pos (origin 0) stays for PSFEx interp. + gal_pos_wcs = np.array( + self._f_wcs_file[exp_name][ccd]["WCS"].all_world2pix( + self.gal_pos[:, 0][ind_obj], + self.gal_pos[:, 1][ind_obj], + 1, + ) + ).T self.interp_PSFs = self.interpsfex( dot_psf_path, @@ -700,7 +803,11 @@ def _interpolate_me(self): array_id = np.concatenate((array_id, np.copy(obj_id))) if self._compute_shape: - self._get_psfshapes() + self._get_psfshapes( + local_wcs_list( + self._f_wcs_file[exp_name][ccd]["WCS"], gal_pos_wcs + ) + ) if array_shape is None: array_shape = np.copy(self.psf_shapes) else: diff --git a/tests/module/test_collate_star_cat.py b/tests/module/test_collate_star_cat.py new file mode 100644 index 000000000..20d52ce23 --- /dev/null +++ b/tests/module/test_collate_star_cat.py @@ -0,0 +1,70 @@ +"""UNIT TESTS FOR STAR-CATALOGUE COLLATION PATHS. + +Pin the patch vs patch-less (v2.0) path and filename convention of +``scripts/python/collate_star_cat.py``. Runs up to v1.6 carry a ``P`` +token in both the input run directory and the output filename; v2.0 is +patch-less (``patch is None``) and drops that token, reading from a single +``/output`` root and writing ``validation_psf_conv-.fits`` — the +name still matched by the downstream ``validation_psf_conv-*`` glob. +""" + +import importlib.util +from pathlib import Path + +import pytest + +# The collation script lives under scripts/python (not an importable package), +# so load it by path. +_SCRIPT = ( + Path(__file__).resolve().parents[2] + / "scripts" + / "python" + / "collate_star_cat.py" +) +_spec = importlib.util.spec_from_file_location("collate_star_cat", _SCRIPT) +collate_star_cat = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(collate_star_cat) + + +@pytest.mark.parametrize( + "patch, exp_input, exp_output", + [ + ("3", "in/P3/output/", "out/P3"), + (None, "in/output/", "out"), + ], +) +def test_collate_paths(patch, exp_input, exp_output): + """v1.x carries the P token; v2.0 (patch None) drops it.""" + assert collate_star_cat.collate_paths("in", "out", patch) == ( + exp_input, + exp_output, + ) + + +@pytest.mark.parametrize( + "patch, expected", + [ + ("3", "validation_psf_conv-3-0.fits"), + (None, "validation_psf_conv-0.fits"), + ], +) +def test_output_filename(patch, expected): + """The patch token is present for v1.x and absent for v2.0.""" + assert collate_star_cat.output_filename("validation_psf", patch, 0) == expected + + +def test_output_filename_matches_downstream_glob(): + """Both layouts stay under the downstream ``validation_psf_conv-*`` glob.""" + for patch in ("1", None): + assert collate_star_cat.output_filename( + "validation_psf", patch, 5 + ).startswith("validation_psf_conv-") + + +@pytest.mark.parametrize("bad", ["v2", "2.0", "v1.7", ""]) +def test_invalid_version_raises(bad): + """A mistyped -V is rejected rather than falling through to v1.x.""" + obj = collate_star_cat.Convert() + obj._params["version_cat"] = bad + with pytest.raises(ValueError): + obj.run() diff --git a/tests/module/test_hsm_sky_coords.py b/tests/module/test_hsm_sky_coords.py new file mode 100644 index 000000000..3d19c0af6 --- /dev/null +++ b/tests/module/test_hsm_sky_coords.py @@ -0,0 +1,128 @@ +"""UNIT TESTS FOR HSM SKY-COORDINATE MOMENTS. + +Pin the equivalence that motivates measuring PSF/star HSM adaptive moments +directly in world coordinates (``FindAdaptiveMom(use_sky_coords=True)``) instead +of measuring in the pixel frame and rotating afterwards with a WCS Jacobian. + +The reference ``_transform_shape`` below is a verbatim copy of the shape-rotation +math that ``scripts/python/collate_star_cat.py`` used to apply (``transform_shape``, +removed in this change). The test draws an elliptical Gaussian on a stamp with a +nontrivial local WCS (rotation + shear + scale, matching a real CD matrix), then +checks that: + + pixel-frame moments → _transform_shape(local WCS Jacobian) + +agrees with the sky-coordinate moments to numerical precision. + +The local WCS here has positive determinant (no parity flip): production CFHT +runs decomposed with ``flip=False`` (the old ``transform_shape`` flip branch, +which is not a correct galsim reflection, never fired), so the faithful +equivalence to pin is the no-flip one. +""" + +import galsim +import numpy as np +import numpy.testing as npt +import pytest + + +def _transform_shape(mom_list, jac): + """Reference: the old pixel-frame → world shape rotation (verbatim). + + Copied from the removed ``collate_star_cat.transform_shape``. ``mom_list`` + is ``[g1, g2, sigma]`` measured in the pixel frame; ``jac`` is the local WCS. + """ + scale, shear, theta, flip = jac.getDecomposition() + + sig_tmp = mom_list[2] * scale + shape = galsim.Shear(g1=mom_list[0], g2=mom_list[1]) + if 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 + + +# A CD-matrix-like local WCS: 0.187 arcsec/pix scale, ~28 deg rotation and a +# small shear, positive determinant (no parity flip). +_SCALE = 0.187 +_THETA = 28.0 * np.pi / 180.0 +_c, _s = np.cos(_THETA), np.sin(_THETA) +_SHEAR = galsim.Shear(g1=0.06, g2=-0.04) +_A = _SHEAR.getMatrix() * _SCALE @ np.array([[_c, -_s], [_s, _c]]) +LOCAL_WCS = galsim.JacobianWCS(_A[0, 0], _A[0, 1], _A[1, 0], _A[1, 1]) + + +@pytest.mark.parametrize( + "e1, e2, sigma_pix", + [ + (0.3, 0.1, 3.0), + (-0.2, 0.25, 2.5), + (0.15, -0.35, 4.0), + (0.0, 0.0, 3.5), + ], +) +def test_sky_coords_moments_match_jacobian_rotation(e1, e2, sigma_pix): + """use_sky_coords=True reproduces pixel-frame moments + Jacobian rotation. + + The local WCS has positive determinant, so ``getDecomposition`` returns + ``flip=False`` and the reference path is a pure scale+rotation+shear — the + transformation ``use_sky_coords`` performs internally. + """ + # Sanity: the reference decomposition sees no parity flip. + assert LOCAL_WCS.getDecomposition()[3] is False + + # An elliptical Gaussian PSF, drawn on a unit-pixel-scale stamp; sigma is in + # pixels because the pixel scale is 1. + stamp = ( + galsim.Gaussian(sigma=sigma_pix) + .shear(e1=e1, e2=e2) + .drawImage(nx=51, ny=51, scale=1.0, method="no_pixel") + ) + + # Pixel-frame measurement, then the old Jacobian rotation. + moms_pix = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array), strict=False + ) + g1_ref, g2_ref, sigma_ref = _transform_shape( + [ + moms_pix.observed_shape.g1, + moms_pix.observed_shape.g2, + moms_pix.moments_sigma, + ], + LOCAL_WCS, + ) + + # Sky-coordinate measurement: same stamp, local WCS attached. + moms_sky = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array, wcs=LOCAL_WCS), + strict=False, + use_sky_coords=True, + ) + + npt.assert_allclose(moms_sky.observed_shape.g1, g1_ref, atol=1e-9, rtol=1e-6) + npt.assert_allclose(moms_sky.observed_shape.g2, g2_ref, atol=1e-9, rtol=1e-6) + npt.assert_allclose(moms_sky.moments_sigma, sigma_ref, rtol=1e-6) + + +def test_sky_coords_sigma_is_scaled_to_world_units(): + """The sky-frame size is the pixel size times the WCS scale (arcsec/pix). + + Pins the deliberate unit change: SIGMA_*_HSM leaves the measurement in world + (arcsec) units rather than pixels once use_sky_coords is on. + """ + stamp = ( + galsim.Gaussian(sigma=3.0) + .drawImage(nx=51, ny=51, scale=1.0, method="no_pixel") + ) + sigma_pix = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array), strict=False + ).moments_sigma + sigma_sky = galsim.hsm.FindAdaptiveMom( + galsim.Image(stamp.array, wcs=LOCAL_WCS), + strict=False, + use_sky_coords=True, + ).moments_sigma + + npt.assert_allclose(sigma_sky, sigma_pix * _SCALE, rtol=1e-6)