diff --git a/astra.yaml b/astra.yaml index 382669d62..0570042a3 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). @@ -2024,7 +2024,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/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/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/unit/test_final_cat_merge_invariants.py b/tests/unit/test_final_cat_merge_invariants.py index 7dcfc6183..033a40d6e 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/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 7b839ba7e..7805a10b3 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": "ea546ec59bcd13c0f8c9ee2c7fa5dde4773975277eed8f46ab63e6d124db0535", + "tile_make_cat": "3cbf638aa39f1574d71bee306b42030246a3a8a9db0749877d7289bead09534a", "tile_merge_cats": "ff21216ea804dccc2d2c290d2b2499d5d05f0c34c0a56993c233f43fe3c06bdb", "tile_merge_headers": "7a344849d62936e2f5598dc8731a2c4947eff2c4c7218b1a731dd2a9577e7111", "tile_ngmix": "5192e65a72b3b6186b29d7bbecfe48751c7bc9ec68a61083ef4092db9f86f02f", @@ -24,7 +24,7 @@ "tile_vignets": "ab52bf2c6ede04915c77f30a44be0cf707c4609ddf8ee768bfaf598cfff31481" }, "schema": 1, - "sha256": "1c792868397e4a2f1144d68d153c4a0a23e9ffd07ecdb95b98896d1b990c8ce1", + "sha256": "d485e25481534d73b2a5c3a40a87f65ac5111e8353968ae6169b7bd1d88fcd97", "unit_pre": { "exp_get_images": "8dec850af212879f225fcf27a5f1281e1a075264c7b97d38c2214395d360168c", "exp_psf": "f2358ddf7385918dc5033d10b37f6dc97a15d02b071a3ea0a4619a5f7e6f5bec", diff --git a/workflow/README.md b/workflow/README.md index 4f0b135c3..89434138b 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 ac9d0760c..180ed29b8 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/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 diff --git a/workflow/rules/tile.smk b/workflow/rules/tile.smk index f607e3617..3dd4e4bdd 100644 --- a/workflow/rules/tile.smk +++ b/workflow/rules/tile.smk @@ -877,7 +877,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: @@ -899,7 +899,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") 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 0d4594ec4..e4a0b4d7c 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, @@ -125,7 +125,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: @@ -151,7 +151,6 @@ def main() -> None: help="names the campaign's group in the output file") p.add_argument("--param-file", required=True, type=Path, help="the input type's final_cat.param — the column list") - 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") @@ -163,7 +162,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: