Write the per-tile final_cat as HDF5, one dataset per column

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).
@@ -2084,7 +2084,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:

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
@@ -21,11 +22,12 @@
 import tqdm
 import glob
 import h5py
-from astropy.io import fits
 
 from cs_util import args as cs_args
 from cs_util import logging
 
+from shapepipe.utilities.final_cat import read_final_cat
+
 
 def params_from_run_config(params, defaults):
     """Fill unset paths from a workflow run config.
@@ -374,9 +376,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 or FITS, decided per file
+    (``shapepipe.utilities.final_cat.read_final_cat``); a FITS table is read
+    at HDU ``params["hdu_num"]`` (default 1). 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,33 +392,17 @@ 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:
-            data = hdu_list[params["hdu_num"]].data
-        except:
-            print(f"Error with ID {id}, file{fits_file}")
-            raise
+    data = read_final_cat(
+        cat_file, params["param_list"], hdu=params.get("hdu_num", 1)
+    )
 
     # 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]
-    if missing:
-        raise KeyError(
-            f"{fits_file}: missing {len(missing)} of the "
-            f"{len(params['param_list'])} requested column(s): "
-            f"{' '.join(missing)}"
-        )
+        params["param_list"] = list(data.dtype.names)
+
     extracted_data = {col: data[col] for col in params["param_list"]}
-    dtype = data.dtype
 
-    return extracted_data, dtype
+    return extracted_data, data.dtype
 
 
 def copy_data(param_list, extracted_data, dtype):
@@ -496,26 +487,36 @@ 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
-    <products>/tiles/<shard>/<tile>/final_cat-<tile>.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.
+    <products>/tiles/<shard>/<tile>/final_cat-<tile>.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)]
     if not all_matches:
         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 +578,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)
 

scripts/python/merge_final_cat.py

@@ -5,13 +5,13 @@
 """Script merge_final_cat.py
 
 Merge all final catalogues, created by ShapePipe module
-``make_catalogue_runner``, into a joined numpy binary file.
+``make_cat_runner`` (HDF5 or FITS, decided per file), into a joined numpy
+binary file.
 
 :Authors: Axel Guinot, Martin Kilbinger
 
 """
 
-from astropy.io import fits
 import numpy as np
 import os
 import sys
@@ -21,6 +21,7 @@
 from optparse import OptionParser
 
 from shapepipe.utilities import cfis
+from shapepipe.utilities.final_cat import read_final_cat
 
 
 class param:
@@ -118,7 +119,7 @@ def parse_options(p_def):
         dest="hdu_num",
         type="int",
         default=p_def.hdu_num,
-        help=f"input HDU number, default='{p_def.hdu_num}'",
+        help=f"input HDU number of FITS catalogues, default='{p_def.hdu_num}'",
     )
 
     parser.add_option(
@@ -247,43 +248,6 @@ def read_param_file(path, verbose=False):
     return param_list
 
 
-def get_data(path, hdu_num, param_list):
-    """Get Data.
-
-    Return data of selected columns from FITS file.
-
-    Parameters
-    ----------
-    path: str
-        input file name
-    hdu_num: int
-        HDU number
-    param_list: list of str
-        parameters to be extracted. If none, copy
-        all columns
-
-    Returns
-    -------
-    numpy array
-        data columns
-
-    """
-    hdu_list = fits.open(path)
-    hdu = hdu_list[hdu_num]
-
-    if param_list:
-        cols = []
-        for p in param_list:
-            cols.append(hdu.columns[p])
-        coldefs = fits.ColDefs(cols)
-        hdu_new = fits.BinTableHDU.from_columns(coldefs)
-        d = hdu_new.data
-    else:
-        d = hdu.data
-
-    return d
-
-
 def main(argv=None):
 
 
@@ -317,9 +281,9 @@ def main(argv=None):
         tile_ID_list = cfis.read_list(param.tile_ID_list_path)
 
     if param.verbose:
-        print("Find input catalogue FITS files")
+        print("Find input catalogue files")
     l = os.listdir(path=path)
-    ext = "fits"
+    ext = (".hdf5", ".fits")
     lpath = []
     for this_l in l:
 
@@ -348,7 +312,7 @@ def main(argv=None):
     count = 0
 
     # Determine number of columns and keys from first catalogue file
-    d_tmp = get_data(lpath[0], param.hdu_num, param.param_list)
+    d_tmp = read_final_cat(lpath[0], param.param_list, hdu=param.hdu_num)
     d = np.zeros(d_tmp.shape, dtype=d_tmp.dtype)
     for key in d_tmp.dtype.names:
         d[key] = d_tmp[key]
@@ -360,7 +324,7 @@ def main(argv=None):
     for fname in lpath[1:]:
 
         try:
-            d_tmp = get_data(fname, param.hdu_num, param.param_list)
+            d_tmp = read_final_cat(fname, param.param_list, hdu=param.hdu_num)
             dd = np.zeros(d_tmp.shape, dtype=d.dtype)
 
             for key in d_tmp.dtype.names: