From d0bf91fde9993fa835b24e4233e2009dd6848d45 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 2 Oct 2026 12:23:29 +0200 Subject: [PATCH 1/3] make_cat: assemble the final catalogue in memory, write HDF5 once make_cat collects every stage's columns (SExtractor, TILE_ID and TILE_UNIQUE_ID, ngmix, per-epoch PSF slots, MASK_*) in one ordered dict and writes final_cat.hdf5 at the end: one lzf-compressed dataset per column, vector columns as 2-D datasets, native dtypes and byte order. With a single write, building on node-local disk has nothing left to buy (the write of a 9,226 x 174 tile takes 0.17 s on candide's NFS), so WORK_DIR and its stale-file handling go. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01JkM3cmXPuuSrNAf7vJde9e --- .../modules/make_cat_package/make_cat.py | 164 ++++++--------- src/shapepipe/modules/make_cat_runner.py | 60 ++---- tests/module/test_make_cat.py | 187 +++++++----------- tests/module/test_make_cat_mask_ext.py | 67 ++----- tests/module/test_read_ext_sexcat.py | 13 +- workflow/config/cfis/config_tile_Mc.ini | 5 - 6 files changed, 164 insertions(+), 332 deletions(-) diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index 8a6f9548a..4b5e6f164 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -9,6 +9,7 @@ import os import re +import h5py import numpy as np from astropy import coordinates as coords from astropy import units as u @@ -42,109 +43,77 @@ def get_output_name(output_dir, file_number_string): output path name """ - return f"{output_dir}/final_cat{file_number_string}.fits" + return f"{output_dir}/final_cat{file_number_string}.hdf5" -def prepare_final_cat_file(output_path, file_number_string): - """Prepare Final Catalogue File. +def write_final_cat(output_path, columns): + """Write Final Catalogue. - Create a ``FITSCatalogue`` object for the current file. + Write the assembled catalogue to HDF5 in one pass: one dataset per + column, in ``columns`` order, each lzf-compressed. A vector column is a + 2-D dataset with one row per object. Parameters ---------- output_path : str - Output file path - file_number_string : str - String with current file numbering - - Returns - ------- - file_io.FITSCatalogue - Output FITS file + Output file path; an existing file is replaced + columns : dict + Column name to array, every array with one entry per object """ - - output_name = get_output_name(output_path, file_number_string) - - return file_io.FITSCatalogue( - output_name, - open_mode=file_io.BaseCatalogue.OpenMode.ReadWrite, - ) - - -def remove_field_name(arr, name): - """Remove Field Name. - - Remove a column of a structured array from the given name. - - Parameters - ---------- - arr : numpy.ndarray - A numpy strucured array - name : str - Name of the field to remove - - Returns - ------- - numpy.ndarray - The structured array with the field removed - - """ - names = list(arr.dtype.names) - if name in names: - names.remove(name) - arr2 = arr[names] - return arr2 + with h5py.File(output_path, "w", track_order=True) as cat: + for name, values in columns.items(): + values = np.asarray(values) + cat.create_dataset( + name, + data=values.astype(values.dtype.newbyteorder("=")), + compression="lzf", + ) -def save_sextractor_data(final_cat_file, sexcat_path, remove_vignet=True): - """Save SExtractor Data. +def read_sextractor_data(sexcat_path, remove_vignet=True): + """Read SExtractor Data. - Save the SExtractor catalogue into the final one, adding the tile as - ``TILE_ID`` (float ``RRR.DDD``) and the survey-wide object ID - ``TILE_UNIQUE_ID`` (:func:`shapepipe.utilities.cfis.get_tile_unique_id` - of the tile and ``NUMBER``). The tile is read from the catalogue's file - name, e.g. ``sexcat-301-279.fits``. + Read the SExtractor catalogue as the final catalogue's first columns, + adding the tile as ``TILE_ID`` (float ``RRR.DDD``) and the survey-wide + object ID ``TILE_UNIQUE_ID`` + (:func:`shapepipe.utilities.cfis.get_tile_unique_id` of the tile and + ``NUMBER``). The tile is read from the catalogue's file name, e.g. + ``sexcat-301-279.fits``. Parameters ---------- - final_cat_file : file_io.FITSCatalogue - Final catalogue sexcat_path : str - Path to SExtractor catalogue to save + Path to SExtractor catalogue remove_vignet : bool - If ``True`` will not save the ``VIGNET`` field into the final catalogue + If ``True`` will not keep the ``VIGNET`` field Returns ------- - int - Number of objects saved + dict + Column name to array, in the SExtractor catalogue's column order @sc [decision:catalogue_assembly.tile_overlap_handling] """ sexcat_file = file_io.FITSCatalogue(sexcat_path, SEx_catalogue=True) sexcat_file.open() data = np.copy(sexcat_file.get_data()) - if remove_vignet: - data = remove_field_name(data, "VIGNET") - cat_size = len(data) + sexcat_file.close() + + columns = { + name: data[name] + for name in data.dtype.names + if not (remove_vignet and name == "VIGNET") + } tile_name = os.path.basename(sexcat_path) nix, niy = cfis.get_tile_number(tile_name) - tile_id_array = np.full(cat_size, float(f"{nix}.{niy}")) - unique_id = cfis.get_tile_unique_id( + columns["TILE_ID"] = np.full(len(data), float(f"{nix}.{niy}")) + columns["TILE_UNIQUE_ID"] = cfis.get_tile_unique_id( cfis.get_tile_id(tile_name), data["NUMBER"] ) - final_cat_file.save_as_fits(data, ext_name="RESULTS") - final_cat_file.open() - final_cat_file.add_cols( - {"TILE_ID": tile_id_array, "TILE_UNIQUE_ID": unique_id} - ) - - sexcat_file.close() - - return cat_size + return columns def parse_mask_ext_paths(paths_str): @@ -172,17 +141,17 @@ def parse_mask_ext_paths(paths_str): return band_paths -def save_mask_ext_data(final_cat_file, band_paths, w_log): +def save_mask_ext_data(final_cat, band_paths, w_log): """Save External Mask Data. Query per-band external healsparse masks at each object's world position and write one ``MASK_`` column per band into the final catalogue. Object positions are read from the SExtractor windowed world coordinates - (``XWIN_WORLD`` = RA, ``YWIN_WORLD`` = Dec, both in degrees) carried in the - ``RESULTS`` extension. Objects falling outside a map's coverage receive - that map's sentinel value (``healsparse.HealSparseMap.get_values_pos`` - returns the map's sentinel — ``-1`` for integer maps — verbatim), which is - the documented off-map flag. + (``XWIN_WORLD`` = RA, ``YWIN_WORLD`` = Dec, both in degrees). Objects + falling outside a map's coverage receive that map's sentinel value + (``healsparse.HealSparseMap.get_values_pos`` returns the map's sentinel — + ``-1`` for integer maps — verbatim), which is the documented off-map + flag. The lookup itself is ``shapepipe.utilities.mask_query.query_map``, shared with the ``mask_query`` module: one primitive, two consumers. @@ -193,36 +162,31 @@ def save_mask_ext_data(final_cat_file, band_paths, w_log): Parameters ---------- - final_cat_file : file_io.FITSCatalogue - Final catalogue + final_cat : dict + Final catalogue columns, updated in place band_paths : dict Mapping from band name to healsparse map path w_log : logging.Logger Logging instance """ - final_cat_file.open() - ra = np.copy(final_cat_file.get_data()["XWIN_WORLD"]) - dec = np.copy(final_cat_file.get_data()["YWIN_WORLD"]) + ra = final_cat["XWIN_WORLD"] + dec = final_cat["YWIN_WORLD"] - mask_cols = {} for band, path in band_paths.items(): w_log.info(f"Query external mask for band {band}: {path}") - mask_cols[f"MASK_{band}"] = mask_query.query_map(path, ra, dec) - final_cat_file.add_cols(mask_cols) - - final_cat_file.close() + final_cat[f"MASK_{band}"] = mask_query.query_map(path, ra, dec) class SaveCatalogue: """Save Catalogue. - Class to save catalogue. + Add shape-measurement and PSF columns to the final catalogue. Parameters ---------- - final_cat_file : str - Final catalogue file name + final_cat : dict + Final catalogue columns, updated in place; must hold ``NUMBER`` cat_size_target : int target catalogue size w_log : logging.Logger @@ -230,9 +194,9 @@ class SaveCatalogue: """ - def __init__(self, final_cat_file, cat_size_target, w_log): + def __init__(self, final_cat, cat_size_target, w_log): - self._final_cat_file = final_cat_file + self._final_cat = final_cat self._cat_size_target = cat_size_target self._w_log = w_log @@ -264,9 +228,7 @@ def process( """ self._output_dict = {} - - self._final_cat_file.open() - self._obj_id = np.copy(self._final_cat_file.get_data()["NUMBER"]) + self._obj_id = self._final_cat["NUMBER"] err_msg = None if mode == "ngmix": @@ -281,9 +243,7 @@ def process( ) if err_msg is None: - self._final_cat_file.add_cols(self._output_dict) - - self._final_cat_file.close() + self._final_cat.update(self._output_dict) return err_msg @@ -522,9 +482,9 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): # Original image PSF (average_original_psf) and metacal # reconvolution kernel (average_multiepoch_psf). Both PSF - # families share ONE write template, so the FITS column - # name (``{obj}``) and the res-key it reads (``{family}``) - # are generated from the same pair and cannot drift apart. + # families share ONE write template, so the column name + # (``{obj}``) and the res-key it reads (``{family}``) are + # generated from the same pair and cannot drift apart. for family, obj in ( ("orig", "PSF_ORIG"), ("reconv", "PSF_RECONV"), @@ -612,7 +572,7 @@ def _save_psf_data(self, galaxy_psf_path, n_epoch_slots=None): """ galaxy_psf_cat = SqliteDict(galaxy_psf_path) - n_epoch = self._final_cat_file.get_data()["N_EPOCH"] + n_epoch = self._final_cat["N_EPOCH"] if n_epoch_slots is None: n_slots = np.max(n_epoch) + 1 else: diff --git a/src/shapepipe/modules/make_cat_runner.py b/src/shapepipe/modules/make_cat_runner.py index 468e7e4e1..540cb8afb 100644 --- a/src/shapepipe/modules/make_cat_runner.py +++ b/src/shapepipe/modules/make_cat_runner.py @@ -6,15 +6,12 @@ """ -import os -import shutil - from shapepipe.modules.make_cat_package import make_cat from shapepipe.modules.module_decorator import module_runner @module_runner( - version="1.1", + version="2.0", input_module=[ "sextractor_runner", "psfex_interp_runner", @@ -26,7 +23,7 @@ "ngmix", ], file_ext=[".fits", ".sqlite", ".fits"], - depends=["numpy", "sqlitedict"], + depends=["numpy", "h5py", "sqlitedict"], ) def make_cat_runner( input_file_list, @@ -63,50 +60,24 @@ def make_cat_runner( else: n_epoch_slots = None - # The catalogue is built in WORK_DIR when set (e.g. node-local disk: - # each save stage rewrites the whole file) and moved to the run's output - # directory once complete. - if config.has_option(module_config_sec, "WORK_DIR"): - work_dir = config.getexpanded(module_config_sec, "WORK_DIR") - os.makedirs(work_dir, exist_ok=True) - else: - work_dir = run_dirs["output"] - work_path = make_cat.get_output_name(work_dir, file_number_string) - # save_as_fits appends to an existing file: a work file left by an - # earlier attempt must not survive into this one. - if os.path.exists(work_path): - os.remove(work_path) - - # Set final output file - final_cat_file = make_cat.prepare_final_cat_file( - work_dir, - file_number_string, - ) - - # Save SExtractor data + # The catalogue is assembled in memory, column by column, and written once. w_log.info("Save SExtractor data") - cat_size_sextractor = make_cat.save_sextractor_data( - final_cat_file, tile_sexcat_path - ) + final_cat = make_cat.read_sextractor_data(tile_sexcat_path) - # Save shape data - sc_inst = make_cat.SaveCatalogue(final_cat_file, cat_size_sextractor, w_log) + sc_inst = make_cat.SaveCatalogue( + final_cat, len(final_cat["NUMBER"]), w_log + ) w_log.info("Save shape measurement data") for shape_type in shape_type_list: w_log.info(f"Save {shape_type.lower()} data") err_msg = sc_inst.process(shape_type.lower(), shape1_cat_path) - - - # If error message: delete (incomplete) output file and raise error + # An incomplete catalogue is never written. if err_msg is not None: - os.remove(work_path) - #raise ValueError(err_msg) w_log.info(err_msg) + return None, None if save_psf: - err_msg = sc_inst.process( - "psf", galaxy_psf_path, n_epoch_slots=n_epoch_slots - ) + sc_inst.process("psf", galaxy_psf_path, n_epoch_slots=n_epoch_slots) # Optional per-band external healsparse mask lookup (UNIONS-WL/spherex#38): # add one MASK_ column per band, queried at each object's world @@ -116,12 +87,11 @@ def make_cat_runner( config.getexpanded(module_config_sec, "MASK_EXT_PATHS") ) w_log.info("Save external mask data") - make_cat.save_mask_ext_data(final_cat_file, band_paths, w_log) + make_cat.save_mask_ext_data(final_cat, band_paths, w_log) - if work_dir != run_dirs["output"] and os.path.exists(work_path): - shutil.move( - work_path, - make_cat.get_output_name(run_dirs["output"], file_number_string), - ) + make_cat.write_final_cat( + make_cat.get_output_name(run_dirs["output"], file_number_string), + final_cat, + ) return None, None diff --git a/tests/module/test_make_cat.py b/tests/module/test_make_cat.py index e29109d19..b5b48548b 100644 --- a/tests/module/test_make_cat.py +++ b/tests/module/test_make_cat.py @@ -17,6 +17,7 @@ import functools import operator +import h5py import numpy as np import numpy.testing as npt import pytest @@ -28,7 +29,6 @@ from shapepipe.modules.make_cat_package.make_cat import SaveCatalogue from shapepipe.modules.make_cat_runner import make_cat_runner from shapepipe.modules.ngmix_package.ngmix import Ngmix -from shapepipe.pipeline import file_io from shapepipe.pipeline.config import CustomParser from shapepipe.utilities import cfis @@ -396,22 +396,12 @@ def _write_galaxy_psf_cat(path, per_obj): db.close() -class _FinalCatStub: - """Stand-in for the FITSCatalogue ``_save_psf_data`` reads N_EPOCH from.""" - - def __init__(self, n_epoch): - self._n_epoch = np.asarray(n_epoch) - - def get_data(self): - return {"N_EPOCH": self._n_epoch} - - def _run_save_psf(galaxy_psf_path, obj_id, n_epoch, n_epoch_slots=None): """Drive ``_save_psf_data`` and return its populated output dict.""" inst = object.__new__(SaveCatalogue) inst._obj_id = np.asarray(obj_id) inst._output_dict = {} - inst._final_cat_file = _FinalCatStub(n_epoch) + inst._final_cat = {"N_EPOCH": np.asarray(n_epoch)} inst._save_psf_data(str(galaxy_psf_path), n_epoch_slots=n_epoch_slots) return inst._output_dict @@ -522,15 +512,7 @@ def _write_sex_like_cat(path, data): def _numbered_data(obj_ids): - """A ``NUMBER`` + dummy second field structured array, one row per id. - - A structured array with a single field gets collapsed by - ``file_io.FITSCatalogue._save_to_fits`` into one row holding a - vector-valued column (``len(names) == 1`` triggers a - ``data = np.array([data])`` wrap), so every SExtractor-like fixture - that ``save_sextractor_data`` re-saves through ``save_as_fits`` needs a - second field to keep one row per object. - """ + """A minimal SExtractor-like table: ``NUMBER`` plus one position field.""" return np.array( [(oid, 0.0) for oid in obj_ids], dtype=[("NUMBER", "i8"), ("X_IMAGE", "f8")], @@ -564,39 +546,36 @@ def test_make_cat_runner_ships_every_detection_unclassified(tmp_path): ) assert result == (None, None) - final_cat = file_io.FITSCatalogue( - str(make_cat.get_output_name(str(tmp_path), "-350-100")) - ) - final_cat.open() - data = final_cat.get_data() - final_cat.close() - npt.assert_array_equal(data["NUMBER"], obj_ids) - assert not [name for name in data.dtype.names if "SPREAD" in name] + with h5py.File(tmp_path / "final_cat-350-100.hdf5", "r") as cat: + npt.assert_array_equal(cat["NUMBER"][()], obj_ids) + assert not [name for name in cat if "SPREAD" in name] + # No MASK_EXT_PATHS, no mask columns. + assert not [name for name in cat if name.startswith("MASK_")] -def test_make_cat_runner_work_dir_publishes_same_catalogue( - tmp_path, monkeypatch -): - """WORK_DIR moves the build elsewhere; the published catalogue is the same. +def test_make_cat_runner_writes_one_hdf5_dataset_per_column(tmp_path): + """All save stages land in one hdf5 file, one lzf dataset per column. - All three save stages run (ngmix, per-epoch PSF slots, one external mask - band). The catalogue built in WORK_DIR is moved to the run's output - directory byte-identical to one built there directly, nothing is left in - WORK_DIR, and a work file a previous attempt left behind does not leak - into it (``save_as_fits`` appends to an existing file). + Runs every stage (ngmix, per-epoch PSF slots, one external mask band). + Columns keep the order the stages add them, a vector column is a 2-D + dataset with one row per object, and each column keeps the dtype its + stage built, in native byte order (the SExtractor input is big-endian + FITS). """ healsparse = pytest.importorskip("healsparse") obj_ids = [1, 2, 3] ra = np.array([10.0, 10.1, 200.0]) dec = np.array([20.0, 20.1, -40.0]) + flux_aper = np.arange(9, dtype=np.float32).reshape(3, 3) sexcat = np.array( - list(zip(obj_ids, [1, 2, 0], ra, dec)), + list(zip(obj_ids, [1, 2, 0], ra, dec, flux_aper)), dtype=[ ("NUMBER", "i8"), ("N_EPOCH", "i8"), ("XWIN_WORLD", "f8"), ("YWIN_WORLD", "f8"), + ("FLUX_APER", "f4", (3,)), ], ) tile_sexcat_path = tmp_path / "tile_sexcat-350-100.fits" @@ -621,48 +600,45 @@ def test_make_cat_runner_work_dir_publishes_same_catalogue( ) mask_path = tmp_path / "mask_r.hsp" mask.write(str(mask_path)) - inputs = [str(tile_sexcat_path), str(galaxy_psf_path), str(ngmix_path)] - stages = { + + config = CustomParser() + config.read_dict({"MAKE_CAT_RUNNER": { "SHAPE_MEASUREMENT_TYPE": "ngmix", "SAVE_PSF_DATA": "True", "N_EPOCH_SLOTS": "3", "MASK_EXT_PATHS": f"r:{mask_path}", - } - - published = {} - for label, work_dir in ( - ("direct", {}), - ("staged", {"WORK_DIR": "$SP_TEST_LOCAL/make_cat"}), - ): - out_dir = tmp_path / label - out_dir.mkdir() - config = CustomParser() - config.read_dict({"MAKE_CAT_RUNNER": {**stages, **work_dir}}) - if work_dir: - local = tmp_path / "local" - monkeypatch.setenv("SP_TEST_LOCAL", str(local)) - (local / "make_cat").mkdir(parents=True) - stale = make_cat.get_output_name(str(local / "make_cat"), "-350-100") - _write_sex_like_cat(stale, _numbered_data([7, 8])) - assert make_cat_runner( - inputs, {"output": str(out_dir)}, "-350-100", config, - "MAKE_CAT_RUNNER", _NullLogger(), - ) == (None, None) - path = make_cat.get_output_name(str(out_dir), "-350-100") - with open(path, "rb") as f: - published[label] = f.read() - - assert published["staged"] == published["direct"] - assert list((tmp_path / "local" / "make_cat").iterdir()) == [] - - with fits.open(make_cat.get_output_name(str(tmp_path / "staged"), "-350-100")) as hdul: - data = hdul[1].data - names = data.columns.names - npt.assert_array_equal(data["NUMBER"], obj_ids) - for col in ("TILE_ID", "NGMIX_MCAL_FLAGS", "HSM_G1_PSF_3", "EXP_ID_2"): - assert col in names, col - npt.assert_allclose(data["HSM_G1_PSF_1"], [0.01, 0.05, -10.0]) - npt.assert_array_equal(data["MASK_r"], [64, 64, -1]) + }}) + out_dir = tmp_path / "output" + out_dir.mkdir() + assert make_cat_runner( + [str(tile_sexcat_path), str(galaxy_psf_path), str(ngmix_path)], + {"output": str(out_dir)}, "-350-100", config, + "MAKE_CAT_RUNNER", _NullLogger(), + ) == (None, None) + assert [p.name for p in out_dir.iterdir()] == ["final_cat-350-100.hdf5"] + + with h5py.File(out_dir / "final_cat-350-100.hdf5", "r") as cat: + names = list(cat) + assert names[:7] == [ + "NUMBER", "N_EPOCH", "XWIN_WORLD", "YWIN_WORLD", "FLUX_APER", + "TILE_ID", "TILE_UNIQUE_ID", + ] + assert names.index("NGMIX_G1_NOSHEAR") < names.index("HSM_G1_PSF_1") + assert names[-1] == "MASK_r" + for name in names: + assert cat[name].shape[0] == len(obj_ids), name + assert cat[name].compression == "lzf", name + assert cat[name].dtype.byteorder in "=|<", name + assert cat["FLUX_APER"].shape == (3, 3) + npt.assert_array_equal(cat["FLUX_APER"][()], flux_aper) + assert cat["NUMBER"].dtype == np.int64 + assert cat["TILE_UNIQUE_ID"].dtype == np.int64 + assert cat["HSM_FLAG_PSF_1"].dtype == np.int16 + assert cat["EXP_ID_2"].dtype == np.int32 + npt.assert_array_equal(cat["NUMBER"][()], obj_ids) + npt.assert_allclose(cat["HSM_G1_PSF_1"][()], [0.01, 0.05, -10.0]) + npt.assert_array_equal(cat["EXP_ID_2"][()], [-1, 2358123, -1]) + npt.assert_array_equal(cat["MASK_r"][()], [64, 64, -1]) @pytest.mark.parametrize("shear", SHEAR_EXTS) @@ -756,27 +732,6 @@ def test_galaxy_cut_admits_only_measured_objects(tmp_path): } -class _ProcessCatStub(_FinalCatStub): - """FITSCatalogue stand-in for ``SaveCatalogue.process``; records add_cols.""" - - def __init__(self, obj_id, n_epoch): - super().__init__(n_epoch) - self._number = np.asarray(obj_id) - self.cols = {} - - def open(self): - pass - - def close(self): - pass - - def get_data(self): - return {"NUMBER": self._number, "N_EPOCH": self._n_epoch} - - def add_cols(self, columns): - self.cols.update(columns) - - def _slot_numbers(out, family): """The slot numbers ``n`` present in ``out`` for ``_n`` columns.""" prefix = f"{family}_" @@ -808,10 +763,9 @@ def test_save_psf_data_fixed_slots_pad_every_family(tmp_path): # Driven through ``process``, the entry point the runner calls. n_slots = 7 - cat = _ProcessCatStub([101, 202, 303], n_epoch=[1, 2, 0]) - sc = SaveCatalogue(cat, 3, _NullLogger()) + out = {"NUMBER": np.array([101, 202, 303]), "N_EPOCH": np.array([1, 2, 0])} + sc = SaveCatalogue(out, 3, _NullLogger()) assert sc.process("psf", str(galaxy_psf_path), n_epoch_slots=n_slots) is None - out = cat.cols n_filled = [1, 2, 0] for family, sentinel in _PSF_SLOT_SENTINELS.items(): @@ -978,7 +932,7 @@ def _write_sexcat(path, number): ]).writeto(path, overwrite=True) -def test_save_sextractor_data_writes_tile_unique_id(tmp_path): +def test_read_sextractor_data_adds_tile_unique_id(tmp_path): """SExtractor-mode final catalogue carries tile_id * 10**6 + NUMBER. ``NUMBER`` is deliberately gapped and unsorted: the ID is built from the @@ -988,28 +942,19 @@ def test_save_sextractor_data_writes_tile_unique_id(tmp_path): sexcat = tmp_path / "sexcat-301-279.fits" _write_sexcat(sexcat, number) - final_cat = make_cat.prepare_final_cat_file(str(tmp_path), "-301-279") - n_obj = make_cat.save_sextractor_data(final_cat, str(sexcat)) - final_cat.close() - - assert n_obj == len(number) - with fits.open(tmp_path / "final_cat-301-279.fits") as hdul: - data = hdul["RESULTS"].data - assert "VIGNET" not in data.names - npt.assert_array_equal(data["NUMBER"], number) - assert data["TILE_UNIQUE_ID"].dtype.newbyteorder("=") == np.int64 - npt.assert_array_equal( - data["TILE_UNIQUE_ID"], 301279 * 10**6 + number - ) - npt.assert_allclose(data["TILE_ID"], 301.279) + data = make_cat.read_sextractor_data(str(sexcat)) + + assert list(data) == ["NUMBER", "XWIN_WORLD", "TILE_ID", "TILE_UNIQUE_ID"] + npt.assert_array_equal(data["NUMBER"], number) + assert data["TILE_UNIQUE_ID"].dtype.newbyteorder("=") == np.int64 + npt.assert_array_equal(data["TILE_UNIQUE_ID"], 301279 * 10**6 + number) + npt.assert_allclose(data["TILE_ID"], 301.279) -def test_save_sextractor_data_refuses_number_beyond_id_range(tmp_path): - """A NUMBER that would overflow into the tile digits raises, writes nothing.""" +def test_read_sextractor_data_refuses_number_beyond_id_range(tmp_path): + """A NUMBER that would overflow into the tile digits raises.""" sexcat = tmp_path / "sexcat-301-279.fits" _write_sexcat(sexcat, np.array([1, 10**6])) - final_cat = make_cat.prepare_final_cat_file(str(tmp_path), "-301-279") with pytest.raises(cfis.CfisError): - make_cat.save_sextractor_data(final_cat, str(sexcat)) - assert not (tmp_path / "final_cat-301-279.fits").exists() + make_cat.read_sextractor_data(str(sexcat)) diff --git a/tests/module/test_make_cat_mask_ext.py b/tests/module/test_make_cat_mask_ext.py index 91d0f5154..96050747f 100644 --- a/tests/module/test_make_cat_mask_ext.py +++ b/tests/module/test_make_cat_mask_ext.py @@ -2,12 +2,11 @@ Exercises the optional per-band external healsparse mask lookup added to ``make_cat`` (PR #847 §4, the ShapePipe end of UNIONS-WL/spherex#38). A small -synthetic ``final_cat`` FITS carrying known ``XWIN_WORLD`` / ``YWIN_WORLD`` +synthetic final catalogue carrying known ``XWIN_WORLD`` / ``YWIN_WORLD`` object positions is queried against synthetic healsparse maps of known value, locking in: (1) the ``MASK_`` column name and per-object values, (2) the off-map sentinel (``-1`` for integer maps) written verbatim for objects outside -coverage, (3) multi-band handling, and (4) that absent config leaves the -catalogue untouched. +coverage, and (3) multi-band handling. """ import numpy as np @@ -17,7 +16,6 @@ healsparse = pytest.importorskip("healsparse") from shapepipe.modules.make_cat_package import make_cat -from shapepipe.pipeline import file_io class _NullLogger: @@ -48,26 +46,13 @@ def _make_map(value, dtype=np.int16, sentinel=-1): return smap -def _write_final_cat(path): - """Write a synthetic final_cat FITS with a RESULTS ext of known positions.""" - data = np.empty( - len(RA), - dtype=[ - ("NUMBER", "i4"), - ("XWIN_WORLD", "f8"), - ("YWIN_WORLD", "f8"), - ], - ) - data["NUMBER"] = np.arange(len(RA)) - data["XWIN_WORLD"] = RA - data["YWIN_WORLD"] = DEC - - cat = file_io.FITSCatalogue( - str(path), - open_mode=file_io.BaseCatalogue.OpenMode.ReadWrite, - ) - cat.save_as_fits(data, ext_name="RESULTS") - return cat +def _final_cat(): + """Synthetic final-catalogue columns with known positions.""" + return { + "NUMBER": np.arange(len(RA)), + "XWIN_WORLD": RA, + "YWIN_WORLD": DEC, + } def test_parse_mask_ext_paths(): @@ -91,39 +76,19 @@ def test_mask_ext_columns(tmp_path): u_map.write(str(u_path)) g_map.write(str(g_path)) - cat_path = tmp_path / "final_cat-000.fits" - _write_final_cat(cat_path) - - cat = file_io.FITSCatalogue( - str(cat_path), - open_mode=file_io.BaseCatalogue.OpenMode.ReadWrite, - ) + cat = _final_cat() make_cat.save_mask_ext_data( cat, {"u": str(u_path), "g": str(g_path)}, _NullLogger(), ) - cat.open() - data = cat.get_data() # On-map objects (first three) carry the map value; the off-map object # (last) carries the map's -1 sentinel. - npt.assert_array_equal(data["MASK_u"], [16, 16, 16, -1]) - npt.assert_array_equal(data["MASK_g"], [32, 32, 32, -1]) + npt.assert_array_equal(cat["MASK_u"], [16, 16, 16, -1]) + npt.assert_array_equal(cat["MASK_g"], [32, 32, 32, -1]) # Integer dtype preserved from the map. - assert np.issubdtype(data["MASK_u"].dtype, np.integer) - cat.close() - - -def test_mask_ext_absent_is_noop(tmp_path): - """Not calling the lookup leaves the catalogue columns unchanged.""" - cat_path = tmp_path / "final_cat-001.fits" - _write_final_cat(cat_path) - - cat = file_io.FITSCatalogue(str(cat_path)) - cat.open() - cols = set(cat.get_data().dtype.names) - cat.close() - - assert cols == {"NUMBER", "XWIN_WORLD", "YWIN_WORLD"} - assert not any(name.startswith("MASK_") for name in cols) + assert np.issubdtype(cat["MASK_u"].dtype, np.integer) + assert list(cat) == [ + "NUMBER", "XWIN_WORLD", "YWIN_WORLD", "MASK_u", "MASK_g" + ] diff --git a/tests/module/test_read_ext_sexcat.py b/tests/module/test_read_ext_sexcat.py index fa406b9a6..7ae89a73c 100644 --- a/tests/module/test_read_ext_sexcat.py +++ b/tests/module/test_read_ext_sexcat.py @@ -6,7 +6,7 @@ extension carrying the tile header, the SExtractor column aliases, one ``VIGNET`` stamp per object cut from the image, and the input ``NUMBER`` kept as is. It follows the catalogue through -``make_cat.save_sextractor_data``, which builds ``TILE_UNIQUE_ID``. The rest +``make_cat.read_sextractor_data``, which builds ``TILE_UNIQUE_ID``. The rest covers the segmentation map: relabelling it to the catalogue's ``NUMBER`` and setting neighbours' ``VIGNET`` pixels to -1e30, as SExtractor does, which is all ngmix reads to mask neighbours. @@ -111,14 +111,11 @@ def test_vignets_are_cut_from_the_image_and_padded_as_sextractor(ldac): assert vignets[2, STAMP // 2, STAMP // 2] == 29 * 1000 + 39 -def test_tile_unique_id_reaches_the_final_catalogue(ldac, tmp_path): +def test_tile_unique_id_reaches_the_final_catalogue(ldac): """make_cat builds the ID from the tile and the catalogue's own NUMBER.""" - final = make_cat.prepare_final_cat_file(str(tmp_path), "-301-279") - n_obj = make_cat.save_sextractor_data(final, str(ldac)) - assert n_obj == len(OBJECTS) - with fits.open(tmp_path / "final_cat-301-279.fits") as hdul: - data = hdul["RESULTS"].data - assert "VIGNET" not in data.names + data = make_cat.read_sextractor_data(str(ldac)) + assert len(data["NUMBER"]) == len(OBJECTS) + assert "VIGNET" not in data npt.assert_array_equal( data["TILE_UNIQUE_ID"], 301279 * 10**6 + np.array([1, 2, 7]) ) diff --git a/workflow/config/cfis/config_tile_Mc.ini b/workflow/config/cfis/config_tile_Mc.ini index addb88861..636803f48 100644 --- a/workflow/config/cfis/config_tile_Mc.ini +++ b/workflow/config/cfis/config_tile_Mc.ini @@ -80,11 +80,6 @@ NUMBERING_SCHEME = -000-000 SHAPE_MEASUREMENT_TYPE = ngmix -# The catalogue is built here and moved to OUTPUT_DIR once complete. Node-local -# ($SP_LOCAL, exported by tile.smk's tile_local() prologue): each save stage -# rewrites the whole file, which is tens of times slower on NFS scratch. -WORK_DIR = $SP_LOCAL/make_cat - # Save per-epoch PSF shapes and epoch identity (HSM_*_PSF_n, EXP_ID_n, # CCD_n). The campaign merge needs one column schema across all tiles, so # the slot count is fixed rather than per-tile max(N_EPOCH)+1. It must From bf9ee8aa2d1a96eb161be7399003b789dee0e059 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 2 Oct 2026 12:23:29 +0200 Subject: [PATCH 2/3] Read per-tile HDF5 final catalogues in the campaign merge create_final_cat.read_data reads a per-tile HDF5 catalogue (one dataset per column) into the same structured array the merge already writes, and still reads FITS catalogues made before the format change, so hand merges over existing runs keep working. find_final_cat looks for .hdf5 before .fits in both layouts. merge_final_cat reads tiles///final_cat-.hdf5 and drops the unused --hdu. The merge invariants now write their tile catalogues with make_cat's own writer, cover a vector column, and check that read_data gives one result from the HDF5 and FITS forms of a catalogue. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01JkM3cmXPuuSrNAf7vJde9e --- scripts/python/create_final_cat.py | 72 +++++++++------- tests/unit/test_final_cat_merge_invariants.py | 85 ++++++++++++++----- workflow/scripts/hdf5_reconcile.py | 4 +- workflow/scripts/merge_final_cat.py | 11 ++- 4 files changed, 113 insertions(+), 59 deletions(-) diff --git a/scripts/python/create_final_cat.py b/scripts/python/create_final_cat.py index 8e989b627..4154621c0 100755 --- a/scripts/python/create_final_cat.py +++ b/scripts/python/create_final_cat.py @@ -2,8 +2,9 @@ """Script create_final_cat.py -Create and update hdf5 file of all final ShapePipe output FITS files, runs of -ShapePipe module ``make_catalogue_runner``. Supercedes `merge_final_cat.py`. +Create and update hdf5 file of all per-tile final catalogues, the output of +ShapePipe module ``make_cat_runner`` (hdf5, one dataset per column; FITS for +catalogues written before that format). Supercedes `merge_final_cat.py`. Usage: in parent dir of patches: create_final_cat.py -p ~/shapepipe/workflow/config/cfis/final_cat.param -i . -P 7 -v -m final_cat_P7.hdf5 @@ -374,9 +375,14 @@ def get_patch_group(hdf5_file, patch, verbose=False): return patch_group -def read_data(fits_file, params): +def read_data(cat_file, params): """Read the parameter list's columns out of one catalogue. + The catalogue is hdf5, one dataset per column (a vector column is 2-D), + or, for one with a ``.fits`` suffix, a FITS table at HDU + ``params["hdu_num"]``. Returns the requested columns and a structured + dtype describing them. + @sc [label:schema] read-data-raises-on-missing-column A requested column the catalogue lacks raises `KeyError` naming it; it is never skipped or filled. `copy_data` keeps only columns present in the @@ -385,31 +391,29 @@ def read_data(fits_file, params): silently narrower, with that slot's exposure identity gone. Enforced by tests/unit/test_final_cat_merge_invariants.py. """ - with fits.open(fits_file) as hdu_list: - try: + if cat_file.endswith(".fits"): + with fits.open(cat_file) as hdu_list: data = hdu_list[params["hdu_num"]].data - except: - print(f"Error with ID {id}, file{fits_file}") - raise + dtype = data.dtype + else: + with h5py.File(cat_file, "r") as cat: + dtype = np.dtype([(col, ds.dtype, ds.shape[1:]) + for col, ds in cat.items()]) + wanted = params["param_list"] or dtype.names + data = {col: cat[col][()] for col in wanted if col in cat} # If columns not given on input: read all column names if params["param_list"] is None: - params["param_list"] = [col for col in data.keys()] - - # RAISE, do not print and fall through. The bare `except:` this replaces - # left extracted_data and dtype unbound, so the caller's own error was an - # UnboundLocalError from the return statement below, naming neither the - # file nor the column that was actually missing. - present = set(data.dtype.names or ()) - missing = [col for col in params["param_list"] if col not in present] + params["param_list"] = list(dtype.names) + + missing = [col for col in params["param_list"] if col not in dtype.names] if missing: raise KeyError( - f"{fits_file}: missing {len(missing)} of the " + f"{cat_file}: missing {len(missing)} of the " f"{len(params['param_list'])} requested column(s): " f"{' '.join(missing)}" ) extracted_data = {col: data[col] for col in params["param_list"]} - dtype = data.dtype return extracted_data, dtype @@ -496,17 +500,23 @@ def collect_tile_ids_image_sims(patch_path): return result +# make_cat_runner writes hdf5; catalogues made before that are FITS. +FINAL_CAT_SUFFIXES = (".hdf5", ".fits") + + def find_final_cat(id, id_path, run_prefix): """Path of this tile's final catalogue, or None. Flat layout first: the unified workflow copies one catalogue per tile to - /tiles///final_cat-.fits, in DOT form and with - no run sub-tree. Then the legacy layout, where the catalogue sits under the - tile's own shapepipe run dir in DASH form, newest run wins. + /tiles///final_cat-.hdf5, in DOT form and + with no run sub-tree. Then the run-dir layout, where the catalogue sits + under the tile's own shapepipe run dir in DASH form, newest run wins. + Either layout may hold an hdf5 or a FITS catalogue; hdf5 is preferred. """ - flat = os.path.join(id_path, f"final_cat-{id}.fits") - if os.path.exists(flat): - return flat + for suffix in FINAL_CAT_SUFFIXES: + flat = os.path.join(id_path, f"final_cat-{id}{suffix}") + if os.path.exists(flat): + return flat base_pattern = os.path.join(id_path, "output", run_prefix) all_matches = [d for d in glob.glob(base_pattern) if os.path.isdir(d)] @@ -514,8 +524,12 @@ def find_final_cat(id, id_path, run_prefix): return None newest_dir = max(all_matches, key=os.path.getmtime) id_dash = re.sub(r"\.", "-", id) - legacy = f"{newest_dir}/make_cat_runner/output/final_cat-{id_dash}.fits" - return legacy if os.path.exists(legacy) else None + for suffix in FINAL_CAT_SUFFIXES: + in_run = (f"{newest_dir}/make_cat_runner/output/" + f"final_cat-{id_dash}{suffix}") + if os.path.exists(in_run): + return in_run + return None def process(params): @@ -577,13 +591,13 @@ def process(params): print(f"Skipping {id} (already processed)") continue - fits_file = find_final_cat(id, id_path, run_prefix) - if fits_file is None: + cat_file = find_final_cat(id, id_path, run_prefix) + if cat_file is None: if params["verbose"]: print(f"Final cat for {id} not found, continuing") continue - extracted_data, dtype = read_data(fits_file, params) + extracted_data, dtype = read_data(cat_file, params) structured_data = copy_data(params["param_list"], extracted_data, dtype) diff --git a/tests/unit/test_final_cat_merge_invariants.py b/tests/unit/test_final_cat_merge_invariants.py index 21eb6adad..7055ec694 100644 --- a/tests/unit/test_final_cat_merge_invariants.py +++ b/tests/unit/test_final_cat_merge_invariants.py @@ -15,6 +15,9 @@ tile catalogue. Objects are matched on ``NUMBER``, so a merge may reorder rows but not move one field of a row without the others. +Tile catalogues are written by make_cat's own writer (``write_final_cat``: +hdf5, one dataset per column), so the producer and the merge meet here. + Failure modes each test was checked against, by editing the code under test: a listed column dropped or an unlisted one kept (columns); one family's slots swapped, or one family's rows reordered (per-epoch); never-fit rows dropped or @@ -22,17 +25,18 @@ (missing column). """ +import importlib.util import sqlite3 import subprocess import sys from pathlib import Path +import h5py import numpy as np import pytest -from numpy.lib.recfunctions import repack_fields +from astropy.io import fits -h5py = pytest.importorskip("h5py") -fits = pytest.importorskip("astropy.io.fits") +from shapepipe.modules.make_cat_package.make_cat import write_final_cat REPO_ROOT = Path(__file__).resolve().parents[2] SCRIPT = REPO_ROOT / "workflow" / "scripts" / "merge_final_cat.py" @@ -48,10 +52,13 @@ ) EPOCH_DTYPE = {"EXP_ID": "i4", "CCD": "i4", "HSM_G1_PSF": "f8", "HSM_G2_PSF": "f8"} +# A vector column: a 2-D dataset in the tile file, a subarray field merged. +VECTORS = (("FLUX_APER", "f4", (3,)),) # In each tile but not in the param file: the merge must leave these behind. UNLISTED = (("MAG_UNLISTED", "f4"), ("EXP_ID_4", "i4"), ("CCD_4", "i4")) PARAM_LIST = [name for name, _ in SCALARS] + [ + name for name, *_ in VECTORS] + [ f"{fam}_{n}" for n in range(1, N_SLOTS + 1) for fam in EPOCH_FAMILIES] @@ -62,7 +69,7 @@ def _slot(fam, n): def _tile_catalogue(seed, rows=40, n_unfit=5): """One tile's catalogue: listed columns, unlisted ones, never-fit rows.""" rng = np.random.default_rng(seed) - columns = list(SCALARS) + [ + columns = list(SCALARS) + list(VECTORS) + [ (_slot(fam, n), EPOCH_DTYPE[fam]) for n in range(1, N_SLOTS + 1) for fam in EPOCH_FAMILIES ] + list(UNLISTED) @@ -71,6 +78,7 @@ def _tile_catalogue(seed, rows=40, n_unfit=5): cat = np.zeros(rows, dtype=columns) cat["NUMBER"] = rng.permutation(rows) + 1 cat["XWIN_WORLD"] = rng.uniform(0, 360, rows) + cat["FLUX_APER"] = rng.uniform(0, 100, (rows, 3)) cat["MAG_UNLISTED"] = rng.uniform(18, 25, rows) cat["EXP_ID_4"] = cat["CCD_4"] = -1 n_epoch = rng.integers(1, N_SLOTS + 1, rows) @@ -103,13 +111,12 @@ def _campaign(root: Path, drop=None): sources = {} for i, tile in enumerate(TILES): cat = _tile_catalogue(seed=1000 + i) - if drop and tile == TILES[-1]: - keep = [c for c in cat.dtype.names if c != drop] - cat = repack_fields(cat[keep]) - path = products / "tiles" / tile[:2] / tile / f"final_cat-{tile}.fits" + columns = {name: cat[name] for name in cat.dtype.names + if not (drop and tile == TILES[-1] and name == drop)} + path = products / "tiles" / tile[:2] / tile / f"final_cat-{tile}.hdf5" path.parent.mkdir(parents=True) - fits.HDUList([fits.PrimaryHDU(), fits.BinTableHDU(cat)]).writeto(path) - sources[tile] = fits.getdata(path, 1) + write_final_cat(str(path), columns) + sources[tile] = cat tile_list = root / "tiles.txt" tile_list.write_text("\n".join(TILES) + "\n") index = root / "index.sqlite" @@ -167,6 +174,14 @@ def test_per_epoch_tuples_survive_the_merge(merged): assert np.array_equal(got[col], want[col]), (tile, col) +def test_vector_columns_merge_as_subarrays(merged): + out, sources = merged + for tile, data in out.items(): + assert data.dtype["FLUX_APER"].shape == (3,), tile + got, want = _by_number(data), _by_number(sources[tile]) + assert np.array_equal(got["FLUX_APER"], want["FLUX_APER"]), tile + + def test_never_fit_rows_pass_through(merged): out, sources = merged for tile, data in out.items(): @@ -191,25 +206,22 @@ def test_a_tile_missing_a_listed_column_stops_the_merge(tmp_path): def _rewrite_tile(tile_path, retype=None): """Rewrite one tile's catalogue, optionally narrowing one column to f4.""" - cat = fits.getdata(tile_path, 1) - arr = np.array(cat) + with h5py.File(tile_path, "r") as cat: + columns = {name: cat[name][()] for name in cat} if retype: - dtype = [(n, "f4" if n == retype else arr.dtype[n]) - for n in arr.dtype.names] - arr = arr.astype(dtype) - fits.HDUList([fits.PrimaryHDU(), fits.BinTableHDU(arr)]).writeto( - tile_path, overwrite=True) + columns[retype] = columns[retype].astype("f4") + write_final_cat(str(tile_path), columns) def test_one_column_type_per_campaign(tmp_path): - """A tile rewritten with the same types refreshes alone (FITS byte order - is not a type change); a tile whose column changes type beside tiles that - kept the old one is refused, naming the column and both dtypes, and the - published catalogue is left as it was.""" + """A tile rewritten with the same types refreshes alone; a tile whose + column changes type beside tiles that kept the old one is refused, naming + the column and both dtypes, and the published catalogue is left as it + was.""" argv, output, _ = _campaign(tmp_path) assert subprocess.run(argv, capture_output=True).returncode == 0 tile = lambda t: (output.parent / "tiles" / t[:2] / t - / f"final_cat-{t}.fits") + / f"final_cat-{t}.hdf5") _rewrite_tile(tile(TILES[0])) run = subprocess.run(argv, capture_output=True, text=True) @@ -223,3 +235,32 @@ def test_one_column_type_per_campaign(tmp_path): assert "NGMIX_T_NOSHEAR" in run.stderr assert "float32" in run.stderr and "float64" in run.stderr, run.stderr assert output.read_bytes() == before + + +def _create_final_cat(): + path = REPO_ROOT / "scripts" / "python" / "create_final_cat.py" + spec = importlib.util.spec_from_file_location("create_final_cat", path) + mod = importlib.util.module_from_spec(spec) + spec.loader.exec_module(mod) + return mod + + +def test_read_data_gives_one_result_from_hdf5_and_fits(tmp_path): + """A hand merge over catalogues made before the hdf5 format (FITS) reads + them into the same structured array as their hdf5 form.""" + cfc = _create_final_cat() + cat = _tile_catalogue(seed=7) + h5 = tmp_path / "final_cat-210-282.hdf5" + write_final_cat(str(h5), {n: cat[n] for n in cat.dtype.names}) + fts = tmp_path / "final_cat-210-282.fits" + fits.HDUList([fits.PrimaryHDU(), fits.BinTableHDU(cat)]).writeto(fts) + + merged = [] + for path in (h5, fts): + params = {"hdu_num": 1, "param_list": list(PARAM_LIST)} + extracted, dtype = cfc.read_data(str(path), params) + merged.append(cfc.copy_data(PARAM_LIST, extracted, dtype)) + from_h5, from_fits = merged + assert from_h5.dtype.names == from_fits.dtype.names == tuple(PARAM_LIST) + for name in PARAM_LIST: + assert np.array_equal(from_h5[name], from_fits[name]), name diff --git a/workflow/scripts/hdf5_reconcile.py b/workflow/scripts/hdf5_reconcile.py index 4fa096d90..e42032759 100644 --- a/workflow/scripts/hdf5_reconcile.py +++ b/workflow/scripts/hdf5_reconcile.py @@ -106,8 +106,8 @@ def stamp(path: Path) -> tuple: def column_types(dtype) -> dict: """``{column: (kind, itemsize, shape)}``: a dataset's schema, as compared - across units. Byte order is left out: FITS sources are big-endian, hdf5 - may hand them back either way, and neither changes a value.""" + across units. Byte order is left out: a source may store either, and + neither changes a value.""" return {n: (dtype[n].base.kind, dtype[n].base.itemsize, dtype[n].shape) for n in dtype.names} diff --git a/workflow/scripts/merge_final_cat.py b/workflow/scripts/merge_final_cat.py index afa05373a..0ddbbc44c 100644 --- a/workflow/scripts/merge_final_cat.py +++ b/workflow/scripts/merge_final_cat.py @@ -18,7 +18,7 @@ WHAT IT REUSES, AND WHAT IT DOES NOT. The column extraction is ``create_final_cat.py``'s — ``read_param_file`` for the parameter list, ``read_data`` and ``copy_data`` for pulling those columns out of one catalogue -with their FITS dtypes — so the column grammar keeps exactly one definition. +with their stored dtypes — so the column grammar keeps exactly one definition. Those three are REPRODUCIBLE FUNCTIONS, and this PR is what made them so: the parameter list comes back ordered rather than through a set, ``copy_data`` allocates the requested columns alone rather than leaving every other column of @@ -29,8 +29,8 @@ Its ``process()`` is NOT used and neither is any of its discovery: that function walks a directory tree the workflow does not have and never will, and it groups by a unit ShapePipe v2 no longer has. This script walks the workflow's own -products tree instead (``tiles/<2-char prefix>//final_cat-.fits``) and -writes the hdf5 itself. +products tree instead (``tiles/<2-char prefix>//final_cat-.hdf5``, +one dataset per column) and writes the merged hdf5 itself. WHERE ``create_final_cat.py`` IS FOUND. Beside this workflow, at ``/scripts/python/create_final_cat.py`` — resolved relative to THIS file, @@ -166,7 +166,7 @@ def catalogues(products_dir: Path, tile_list: Path, index_db: Path) -> list: out, missing = [], [] for tile in sorted(build_index.campaign_tiles(tile_list, index_db)): path = (products_dir / "tiles" / tile[:2] / tile - / f"final_cat-{tile}.fits") + / f"final_cat-{tile}.hdf5") if path.exists(): out.append((tile, path)) else: @@ -195,7 +195,6 @@ def main() -> None: p.add_argument("--tile-detection", required=True, choices=TILE_DETECTIONS, help="where the tile detections came from (config " "tile_detection); selects the columns requested") - p.add_argument("--hdu", type=int, default=1) p.add_argument("--snapshot-json", type=Path, default=None, help="sp run's code snapshot (bin/sp's " "$STATE_DIR/code/snapshot.json); absent outside sp run") @@ -209,7 +208,7 @@ def main() -> None: sys.exit(f"merge_final_cat: no columns read from {args.param_file}") # read_data/copy_data read their knobs out of this dict, exactly as # create_final_cat.py's own main() builds it. - params = {"hdu_num": args.hdu, "param_list": param_list, "verbose": False} + params = {"param_list": param_list, "verbose": False} tiles = catalogues(args.products_dir, args.tile_list, args.index_db) if not tiles: From ce35edd6fedbe8d01298378d4bfd504527f8362f Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Fri, 2 Oct 2026 12:23:29 +0200 Subject: [PATCH 3/3] Declare the per-tile final_cat as final_cat-.hdf5 tile_make_cat publishes final_cat-.hdf5 and final_cat_merge reads it. The rendered tile_make_cat shell changes, so the params pin moves: this is a campaign boundary. astra.yaml records the format. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01JkM3cmXPuuSrNAf7vJde9e --- astra.yaml | 4 ++-- tests/workflow/harness.py | 2 +- tests/workflow/params_pin.json | 4 ++-- workflow/README.md | 2 +- workflow/Snakefile | 2 +- workflow/rules/tile.smk | 4 ++-- 6 files changed, 9 insertions(+), 9 deletions(-) diff --git a/astra.yaml b/astra.yaml index e566a1253..cd03cf9b4 100644 --- a/astra.yaml +++ b/astra.yaml @@ -69,7 +69,7 @@ inputs: outputs: - id: final_cat type: data - format: fits + format: hdf5 description: >- Per-tile shear catalogue family, the terminal science product (make_cat_runner). @@ -1994,7 +1994,7 @@ analyses: outputs: - id: tile_final_cat type: data - format: fits + format: hdf5 description: The assembled per-tile science catalogue (final_cat family). inputs: [ngmix_chunks] decisions: diff --git a/tests/workflow/harness.py b/tests/workflow/harness.py index 1c8420f26..b82cc0308 100644 --- a/tests/workflow/harness.py +++ b/tests/workflow/harness.py @@ -155,7 +155,7 @@ def tile_manifest(self, tile, stage): def final_cat(self, tile): """Return the expected persistent catalogue for one tile.""" return (self.products_dir / "tiles" / tile[:2] / tile - / f"final_cat-{tile}.fits") + / f"final_cat-{tile}.hdf5") def persist_manifest(self, exp): """Return the expected persistent manifest for one exposure.""" diff --git a/tests/workflow/params_pin.json b/tests/workflow/params_pin.json index eeb61d950..f5674c2e2 100644 --- a/tests/workflow/params_pin.json +++ b/tests/workflow/params_pin.json @@ -16,7 +16,7 @@ "tile_find_exposures": "8704317871744996c44351c2836fcb222d7986a602e9e046a90d684ba7b3c184", "tile_get_catalogue": "7e3f889a955a14b0b917a015c2a8bc90e433a89b1cc5e00b8ca94695c4c93d3c", "tile_get_images": "331a67e747f211ebf4c14b946a7af7f9fc9f55243d69d2791d74aecc3ca228c3", - "tile_make_cat": "a57518b04c11f70bf41b320532fddffdd1115fdac604d7dd3928eb61cca28f23", + "tile_make_cat": "ad1256f27109ec443b6ab870eb04550ff260d0f365cb71864b4c33637a07f3df", "tile_merge_cats": "ff21216ea804dccc2d2c290d2b2499d5d05f0c34c0a56993c233f43fe3c06bdb", "tile_merge_headers": "7a344849d62936e2f5598dc8731a2c4947eff2c4c7218b1a731dd2a9577e7111", "tile_ngmix": "6ed10a3a4d3658ba0303fb0c23ad5100ec8cbdc638ca36c3b25875087b5c3415", @@ -24,7 +24,7 @@ "tile_vignets": "9d4ae0d99c18217f2185f245281312454c8a219ec1628176e08c71a5efc4dc91" }, "schema": 1, - "sha256": "f9772987e517e9827c8d501d310ec323be251a62b7b89aad0725a03f8f4f120e", + "sha256": "83b343ac61586714df2b916845f0ee43e3466ed826a38b5cd65f5e7d65bb363d", "unit_pre": { "exp_get_images": "8dec850af212879f225fcf27a5f1281e1a075264c7b97d38c2214395d360168c", "exp_psf": "f2358ddf7385918dc5033d10b37f6dc97a15d02b071a3ea0a4619a5f7e6f5bec", diff --git a/workflow/README.md b/workflow/README.md index a1654fcb4..bb3ddc871 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -396,7 +396,7 @@ profiles/nibi/config.yaml SLURM executor; apptainer SDM; per-user jobs cap; kee apart, and nothing else would notice if they did — a column added to one writer would just be missing from the other's product. `tests/unit/` `test_star_cat_columns.py` is what holds them together. - `final_cat_merge` collects every ready tile's `final_cat-.fits` into + `final_cat_merge` collects every ready tile's `final_cat-.hdf5` into `/final_cat_.hdf5`: one dataset per tile under a group named for the campaign, the `final_cat.param` columns, an `n_tiles` attribute. That schema is what sp_validation's reader opens, so it is fixed; the column diff --git a/workflow/Snakefile b/workflow/Snakefile index 6aed81cff..31d176098 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -400,7 +400,7 @@ def final_cat(tile): stores that reclamation deleted — the rerun avalanche the cut exists to prevent, arriving by way of the purge instead. """ - return f"{PRODUCTS_DIR}/tiles/{tile[:2]}/{tile}/final_cat-{tile}.fits" + return f"{PRODUCTS_DIR}/tiles/{tile[:2]}/{tile}/final_cat-{tile}.hdf5" def unit_num(unit): """$SP_UNIT_NUM: ShapePipe's image-number convention, dot -> dash, leading diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index f105d23b4..22a61a11a 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -873,7 +873,7 @@ rule tile_make_cat: ms = rules.tile_merge_cats.output.manifest, output: manifest = f"{TILE_DIR}/manifests/tile_make_cat.json", - final_cat = f"{PROD_TILE_DIR}/final_cat-{{tile}}.fits", + final_cat = f"{PROD_TILE_DIR}/final_cat-{{tile}}.hdf5", log: f"{TILE_DIR}/logs/tile_make_cat.json" params: @@ -895,7 +895,7 @@ rule tile_make_cat: sp_shell("tile_make_cat", "config_tile_Mc.ini", post="if [ $rc -eq 0 ]; then\n" ' cp -f "$(ls -1 "$SP_RUN"/output/run_sp_tile_Mc/make_cat_runner' - '/output/final_cat*.fits | head -1)" {output.final_cat}\n' + '/output/final_cat*.hdf5 | head -1)" {output.final_cat}\n' "fi\n")