Exposure-level HealSparse maps: exposure count and flagged pixels

astra.yaml

@@ -231,7 +231,61 @@ analyses:
           catalogue's MASK_<flag value>_<name> columns.
         inputs: [exposure_flags, sky_masks]
         decisions: [pixel_mask_source, psf_star_mask_veto, sky_mask_application]
+      - id: exposure_maps
+        type: data
+        format: hsp
+        description: >-
+          nexp_<run>.hsp and nflagged_<run>.hsp at the mask ladder's
+          resolution: per sky pixel, the number of exposures whose CCD with a
+          valid PSF model covers it, and the flagged CCD pixels of those
+          exposures that fall in it. Products; nothing in the workflow reads
+          them.
+        inputs: [exposure_flags]
+        decisions: [nexp_map_valid_psf_ccds, defect_map_from_flags]
     decisions:
+      nexp_map_valid_psf_ccds:
+        label: The exposure count covers only CCDs with a valid PSF model
+        rationale: >-
+          exp_maps takes the CCDs whose validation_psf catalogue exp_persist
+          packed (psfex_interp writes it only when the fit succeeds) and marks
+          the sky pixels whose centres lie inside each CCD's corner polygon.
+          The CCDs of one exposure do not overlap, so summing fragments counts
+          exposures that can contribute a PSF-corrected epoch.
+        default: valid_psf_ccds
+        options:
+          valid_psf_ccds:
+            label: Count exposures whose CCD has a valid PSF model
+          all_ccds:
+            label: Count every exposure whose CCD covers the pixel
+            excluded: true
+            excluded_reason: >-
+              Counts epochs that cannot contribute a shape, overstating depth
+              exactly where the PSF fit failed.
+      defect_map_from_flags:
+        label: Instrument flags counted per sky pixel, inside the imaging area
+        rationale: >-
+          The per-CCD flag images otherwise never leave the pixel domain, so a
+          footprint built from CCD corners would include defective pixels.
+          exp_maps bins the centre of every nonzero-flag pixel inside DATASEC
+          into the sky pixel it falls in, and nflagged sums the counts over
+          exposures (~74 MegaCam pixels per sky pixel). The overscan border
+          (flag 3, 3.6% of a split) lies outside DATASEC and outside the
+          coverage polygon. On exposures 2603236/2603237 the remaining flags
+          are 0.28% of the imaging pixels; any-touch rasterization of them
+          masks 2.2% of the sky, because single-pixel columns and cosmics
+          widen to 1.6", so the map keeps the count and leaves the threshold
+          to the consumer (nflagged > 0 is the any-touch mask).
+        default: count_flagged_pixels
+        options:
+          count_flagged_pixels:
+            label: Count flagged CCD pixels per sky pixel, summed over exposures
+          any_touch_boolean:
+            label: Boolean mask of every sky pixel a flagged pixel touches
+            excluded: true
+            excluded_reason: >-
+              Widens one-pixel defects to the 1.6" sky pixel and loses how
+              much of each sky pixel is flagged; recoverable from the count as
+              nflagged > 0.
       pixel_mask_source:
         label: >-
           Of the external masks, only the instrument flag image reaches pixels

tests/unit/test_campaign_lineage.py

@@ -39,6 +39,9 @@
     "prod_exp_dir": (("2605805",), False),
     "prod_exp_manifest": (("2605805", "exp_persist"), False),
     "prod_exp_tar": (("2605805",), False),
+    "prod_exp_maps": (("2605805",), False),
+    "nexp_map": ((), True),
+    "nflagged_map": ((), True),
 }
 PRODUCT_TEMPLATES = ("PROD_TILE_DIR", "PROD_EXP_DIR")
 

tests/unit/test_exp_maps.py

@@ -0,0 +1,197 @@
+"""The per-exposure fragment and the campaign sum (``exp_maps.py``,
+``merge_exposure_maps.py``).
+
+The fixture is a synthetic store of two CCDs on a rotated TAN WCS at MegaCam's
+0.187"/pixel, so a sky pixel covers ~74 CCD pixels as on the sky. Its flag
+image carries a one-pixel bad column, a saturated blob, isolated pixels on the
+CCD edges, a lattice of hot pixels, and a flagged overscan border outside
+DATASEC, which must count for nothing.
+
+Needs healsparse, hpgeom and astropy, so it runs inside the container and
+skips outside.
+"""
+
+import importlib.util
+import json
+import sys
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+pytestmark = [pytest.mark.unions,
+              pytest.mark.decision("masking.defect_map_from_flags")]
+
+SCRIPTS = Path(__file__).resolve().parents[2] / "workflow" / "scripts"
+EXP = "2079612"
+NY, NX = 320, 240                 # DATASEC
+PAD_X, PAD_Y = 8, 6               # overscan columns either side, rows on top
+SCALE_DEG = 0.187 / 3600.0
+ROTATION_DEG = 23.0
+
+
+def _load(name):
+    sys.path.insert(0, str(SCRIPTS))
+    try:
+        spec = importlib.util.spec_from_file_location(name, SCRIPTS / f"{name}.py")
+        module = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(module)
+        return module
+    finally:
+        sys.path.remove(str(SCRIPTS))
+
+
+@pytest.fixture(scope="module")
+def maps():
+    for dep in ("healsparse", "hpgeom", "astropy"):
+        pytest.importorskip(dep)
+    return _load("exp_maps")
+
+
+def _wcs(crval):
+    from astropy.wcs import WCS
+
+    theta = np.deg2rad(ROTATION_DEG)
+    wcs = WCS(naxis=2)
+    wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
+    wcs.wcs.crval = list(crval)
+    wcs.wcs.crpix = [PAD_X + NX / 2 + 0.5, NY / 2 + 0.5]
+    wcs.wcs.cd = SCALE_DEG * np.array([[-np.cos(theta), np.sin(theta)],
+                                       [np.sin(theta), np.cos(theta)]])
+    return wcs
+
+
+def _flags():
+    """The DATASEC flags; ``_split`` adds the overscan border."""
+    flags = np.zeros((NY, NX), dtype=np.int16)
+    flags[:, 101] = 1                                       # bad column
+    yy, xx = np.mgrid[:NY, :NX]
+    flags[(yy - 210) ** 2 + (xx - 60) ** 2 <= 7 ** 2] = 2   # saturated blob
+    for row, col in [(0, 0), (NY - 1, NX - 1), (40, NX - 1), (77, 33)]:
+        flags[row, col] = 8                                 # isolated pixels
+    flags[5::11, 3::11] = 8                                 # sparse hot pixels
+    return flags
+
+
+def _split(flags):
+    """The full split: DATASEC flags inside a border flagged 3."""
+    full = np.full((NY + PAD_Y, NX + 2 * PAD_X), 3, dtype=np.int16)
+    full[:NY, PAD_X:PAD_X + NX] = flags
+    return full
+
+
+def _store(tmp_path, ccds, crvals, valid):
+    """A split dir with image (header only) and flag splits for ``ccds``, and
+    an exp_persist manifest naming the ``valid`` CCDs' validation_psf files."""
+    from astropy.io import fits
+
+    split = tmp_path / "store" / "output/run_sp_exp_Sp/split_exp_runner/output"
+    split.mkdir(parents=True)
+    for ccd, crval in zip(ccds, crvals):
+        header = _wcs(crval).to_header()
+        header["DATASEC"] = f"[{PAD_X + 1}:{PAD_X + NX},1:{NY}]"
+        fits.PrimaryHDU(header=header).writeto(split / f"image-{EXP}-{ccd}.fits")
+        fits.PrimaryHDU(data=_split(_flags())).writeto(
+            split / f"flag-{EXP}-{ccd}.fits")
+    manifest = tmp_path / "exp_persist.json"
+    manifest.write_text(json.dumps({"files": [
+        {"name": f"validation_psf-{EXP}-{c}.fits"} for c in valid]
+        + [{"name": f"psf-{EXP}-0.psf"}]}))
+    return tmp_path / "store", manifest
+
+
+def _run(maps, tmp_path, monkeypatch, **store):
+    import healsparse as hsp
+
+    exp_dir, persist = _store(tmp_path, **store)
+    fragment = tmp_path / "maps" / f"maps-{EXP}.hsp"
+    monkeypatch.setattr(sys, "argv", [
+        "exp_maps.py", "--exp-dir", str(exp_dir), "--exp", EXP,
+        "--persist-manifest", str(persist), "--fragment", str(fragment),
+        "--manifest", str(tmp_path / "exp_maps.json")])
+    maps.main()
+    return hsp.HealSparseMap.read(str(fragment)), json.loads(
+        (tmp_path / "exp_maps.json").read_text())
+
+
+def _sky_pixels(maps, crval, x, y):
+    """Sky pixels of DATASEC pixel positions (0-based within DATASEC)."""
+    import hpgeom as hpg
+
+    ra, dec = _wcs(crval).pixel_to_world_values(np.asarray(x) + PAD_X, y)
+    return hpg.angle_to_pixel(maps.NSIDE, ra, dec, nest=True)
+
+
+def test_fragment_counts_every_flagged_pixel(maps, tmp_path, monkeypatch):
+    """Each flagged DATASEC pixel is counted once, in the sky pixel holding
+    its centre; the overscan border counts for nothing."""
+    crval = (150.3, 31.7)
+    fragment, record = _run(maps, tmp_path, monkeypatch,
+                            ccds=[0], crvals=[crval], valid=[0])
+    rows, cols = np.nonzero(_flags())
+    expected = dict(zip(*np.unique(_sky_pixels(maps, crval, cols, rows),
+                                   return_counts=True)))
+    pixels = fragment.valid_pixels
+    values = fragment.get_values_pix(pixels).astype(int)
+    got = {p: v - 1 for p, v in zip(pixels, values) if v > 1}
+    # Centres on the very edge of DATASEC can fall in a sky pixel whose own
+    # centre is outside the coverage polygon; those are dropped.
+    assert set(got) <= set(expected)
+    assert all(got[p] == expected[p] for p in got)
+    assert sum(got.values()) > 0.97 * rows.size
+    assert record["ccds"]["0"]["flagged_pixels"] == rows.size
+
+
+def test_coverage_is_the_imaging_area_of_valid_psf_ccds(maps, tmp_path,
+                                                        monkeypatch):
+    """Interior pixel centres are covered, the overscan is not, and a CCD
+    without a PSF model contributes nothing."""
+    crvals = [(150.3, 31.7), (150.5, 31.7)]
+    fragment, record = _run(maps, tmp_path, monkeypatch,
+                            ccds=[0, 1], crvals=crvals, valid=[0])
+    yy, xx = np.mgrid[13:NY - 13:7, 13:NX - 13:7]   # > one sky pixel in
+    inside = _sky_pixels(maps, crvals[0], xx.ravel(), yy.ravel())
+    assert (fragment.get_values_pix(inside) > 0).all()
+    overscan = _sky_pixels(maps, crvals[0], np.full(50, -PAD_X - 4.0),
+                           np.linspace(20, NY - 20, 50))
+    assert (fragment.get_values_pix(overscan) == 0).all()
+    other = _sky_pixels(maps, crvals[1], xx.ravel(), yy.ravel())
+    assert (fragment.get_values_pix(other) == 0).all()
+    assert list(record["ccds"]) == ["0"]
+
+
+def test_merge_counts_exposures_and_flagged_pixels(maps, tmp_path,
+                                                   monkeypatch):
+    """Two fragments overlapping on half their pixels: nexp is 2 in the
+    overlap, nflagged sums the counts, and a fragment-less exposure is
+    skipped."""
+    import healsparse as hsp
+
+    merge = _load("merge_exposure_maps")
+    products = tmp_path / "products"
+    a = np.arange(1000, dtype=np.int64) + 10**9
+    for exp, pix in (("2000001", a), ("2000002", a + 500)):
+        frag = hsp.HealSparseMap.make_empty(maps.NSIDE_COVERAGE, maps.NSIDE,
+                                            np.uint8)
+        frag[pix] = np.ones(pix.size, np.uint8)
+        frag[pix[500:510]] = np.full(10, 1 + 7, np.uint8)
+        path = merge.fragment_path(products, exp)
+        path.parent.mkdir(parents=True)
+        frag.write(str(path))
+    monkeypatch.setattr(merge.build_index, "campaign_exposures",
+                        lambda *_: ["2000001", "2000002", "2000003"])
+    out = {k: tmp_path / f"{k}.hsp" for k in ("nexp", "nflagged")}
+    monkeypatch.setattr(sys, "argv", [
+        "merge_exposure_maps.py", "--products-dir", str(products),
+        "--tile-list", "-", "--index-db", "-",
+        "--nexp", str(out["nexp"]), "--nflagged", str(out["nflagged"])])
+    merge.main()
+
+    nexp = hsp.HealSparseMap.read(str(out["nexp"]))
+    nflagged = hsp.HealSparseMap.read(str(out["nflagged"]))
+    assert nexp.get_values_pix(a[:500]).tolist() == [1] * 500
+    assert nexp.get_values_pix(a[500:]).tolist() == [2] * 500
+    assert nexp.get_values_pix(a[500:] + 500).tolist() == [1] * 500
+    assert nflagged.get_values_pix(a[500:510]).tolist() == [7] * 10
+    assert nflagged.get_values_pix(a[500:510] + 500).tolist() == [7] * 10
+    assert nflagged.valid_pixels.size == 20

tests/workflow/harness.py

@@ -157,10 +157,10 @@ def final_cat(self, tile):
         return (self.products_dir / "tiles" / tile[:2] / tile
                 / f"final_cat-{tile}.fits")
 
-    def persist_manifest(self, exp):
-        """Return the expected persistent manifest for one exposure."""
+    def persist_manifest(self, exp, stage="exp_persist"):
+        """Return one exposure's manifest on the persistent root."""
         return (self.products_dir / "exp" / exp[:2] / exp
-                / "manifests" / "exp_persist.json")
+                / "manifests" / f"{stage}.json")
 
 
 @dataclass

tests/workflow/params_pin.json

@@ -5,9 +5,11 @@
     "clean_exposure": "22cb76b13a5205d20a02a9bd3b8c8bea5ea2801555e24dd2b84b11145f7e79d9",
     "clean_tile": "a5c07b0461526ed407df36a291deb866d4181524b4c8e0fbc3dd047fd9d28479",
     "exp_get_images": "a9e9eca0a99a348e43e0cd44d0c0aae0f890773d03e35007fb77f886b7d21f5f",
+    "exp_maps": "6114b2bb508217eb1a87c89a295485c4cc17f329391a27832548f525a3ef1d2e",
     "exp_persist": "302e2837542bc1102430c27c81c600b7cda32e8bddcb5fd60d33950987609fff",
     "exp_psf": "87b9fd799a54130eb3dd9c78bb5124838854c23f365df7a991a53f4c6ac4c3ff",
     "exp_split": "8fc97dd48e76cd6c7fdf9ee388ba1c1a3dd63428cbe1961641fbebe76b05786e",
+    "exposure_maps": "15bc1023a759ad61b67cc27f1740a74bb3be4d268969073820766fb64d6aace1",
     "final_cat_merge": "e7f46859c4503a2220713d7bb2507555515d0a9632d780b20f14c59e32210023",
     "prepare_all_tiles": "b8f872a22adf014e25a7fa5198f49b71a6fe9e56042ed82b682bc8763970a844",
     "star_cat_merge": "6277450958474af5270982fa35360f2f237a29f7533c526ee9265dfd5acc07a0",
@@ -24,7 +26,7 @@
     "tile_vignets": "89d95bd9107c124e86dfb029d5ba2a2fc1871311a71a8fe18466302ef2213b14"
   },
   "schema": 1,
-  "sha256": "eab6930e0fd1044b9ce7a99172dbf9cbdcdec21291f8fa640d6e60c2efd3dd03",
+  "sha256": "23b5477740ce0ffbe6e86ee3427ed1fa9d42473519a0b8137767b8fd1ad356b9",
   "unit_pre": {
     "exp_get_images": "dac6685a207dae3ea81d296ce53e9ba75d4636cf1f81e6399f2dda452680c992",
     "exp_psf": "87fc8ea1153a709ab7ba56542cfacc041f5f8b7f6cf2f38a9bc9d8004d1cecb4",

tests/workflow/test_dag.py

@@ -17,7 +17,7 @@
     "final_cat_merge",
 }
 TILE_SHAPE_STORE_RULES = ("tile_vignets", "tile_ngmix", "tile_make_cat")
-PSF_RULES = {"exp_persist", "star_cat_merge"}
+PSF_RULES = {"exp_persist", "star_cat_merge", "exp_maps", "exposure_maps"}
 CATALOGUE_RULES = {"tile_get_catalogue"}
 
 
@@ -63,6 +63,7 @@ def test_clean_exposure_waits_on_persist_iff_psf(campaign, dag):
         ]
         if campaign.psf_model != "fake":
             expected.append(campaign.persist_manifest(exp))
+            expected.append(campaign.persist_manifest(exp, "exp_maps"))
         assert Counter(map(str, job.input)) == Counter(map(str, expected)), (
             "clean-exposure-waits-on-persist-iff-psf", exp, list(job.input)
         )

universes/committed.yaml

@@ -10,6 +10,8 @@ analyses:
       pixel_mask_source: instrument_flags_only
       psf_star_mask_veto: instrument_flags_only
       sky_mask_application: catalogue_columns
+      nexp_map_valid_psf_ccds: valid_psf_ccds
+      defect_map_from_flags: count_flagged_pixels
   detection:
     decisions:
       tile_detection: unions_catalogue

workflow/CONTRACTS

@@ -9,7 +9,8 @@ prose names it.
 A campaign's name is the run config's `run:` and nothing else. The Snakefile
 binds it once, `CAMPAIGN = config["run"]`, and every product that carries a
 campaign name takes it from there: `final_cat_<run>.hdf5` and its
-`patches/<run>` group, `full_starcat_<run>.hdf5`, and the `$run` in the
+`patches/<run>` group, `full_starcat_<run>.hdf5`, `nexp_<run>.hsp`,
+`nflagged_<run>.hsp`, and the `$run` in the
 machine defaults' `products_dir`/`index_db`. No rule or script reads a
 `campaign` key or derives a name from a directory (`PRODUCTS_DIR.name`); two
 sources that can disagree would file one campaign's merge under another's name.

workflow/README.md

@@ -69,8 +69,9 @@ each overlay file to its `cfis/` original confined to input naming.
 `psf_model: fake` is the simulations' true PSF: the exposure stage runs only
 SExtractor (for the background maps the vignets read), and `tile_vignets` runs
 `fake_interp_runner`, which writes the `galaxy_psf` product from `psf_dict`.
-With no PSF model there is nothing to persist per exposure, so `exp_persist` and
-`star_cat_merge` do not run and `clean_exposure` does not wait on them.
+With no PSF model there is nothing to persist per exposure, so `exp_persist`,
+`exp_maps` and the two merges after them do not run and `clean_exposure` does
+not wait on them.
 Simulations that contain stars can run `psfex` as the data do. `mccd` is
 refused until `persist_exp.py` and `merge_star_cat.py` read MCCD products.
 
@@ -239,7 +240,7 @@ workflow/
   bin/sp                 committed launcher (module load + /project venv + launch code snapshot + run/report/container/cancel)
   rules/
     prepare.smk          tile get_images/uncompress/find_exposures
-    exposure.smk         per-exposure: get_images, split, psf, persist (no temp()); campaign star_cat_merge
+    exposure.smk         per-exposure: get_images, split, psf, persist, maps (no temp()); campaign star_cat_merge, exposure_maps
     tile.smk             per-tile: exp forest, merge_headers, detect (SExtractor, joined to the UNIONS catalogue on data), vignets, ngmix, merge, make_cat; campaign final_cat_merge
   scripts/
     build_index.py       prepare-phase run_index.sqlite builder (plain script)
@@ -252,6 +253,8 @@ workflow/
     merge_star_cat.py    ALL exposures' validation_psf, out of the tars -> full_starcat_<run>.hdf5
     merge_final_cat.py   ALL tiles' final_cat -> final_cat_<run>.hdf5 (the final_cat_merge rule)
     clean_exposure.py    ONE exposure's store + manifests + logs -> tombstone (the clean_exposure rule)
+    exp_maps.py          ONE exposure's valid-PSF CCD footprints + flagged-pixel counts -> healsparse fragment (the exp_maps rule)
+    merge_exposure_maps.py  ALL fragments -> nexp_<run>.hsp, nflagged_<run>.hsp (the exposure_maps rule)
 profiles/nibi/config.yaml  SLURM executor; apptainer SDM; per-user jobs cap; keep-going
 ```
 
@@ -302,6 +305,20 @@ profiles/nibi/config.yaml  SLURM executor; apptainer SDM; per-user jobs cap; kee
   catalogue server, staged, or rasterized, which is why the old
   `star_catalogue` / `exp_star_cat` / `exp_mask` rules and their cache root are
   gone.
+- **Two HealSparse maps come out of the exposures.** At the mask ladder's
+  resolution (nside 131072), `nexp_<run>.hsp` counts per sky pixel the
+  exposures whose CCD with a valid PSF model covers it, and
+  `nflagged_<run>.hsp` counts the flagged CCD pixels (bad columns, saturated
+  pixels, bleed trails, cosmics) of those exposures that fall in it, ~74 CCD
+  pixels filling one sky pixel. The flag image is the one mask that otherwise
+  never leaves the pixel domain, so a footprint built from CCD corners could not
+  subtract it (#878). `exp_maps` writes one fragment per exposure beside its PSF
+  tar, from the split's `DATASEC` (the overscan border is neither coverage nor
+  defect) and the CCDs `exp_persist` packed a PSF for; `exposure_maps` sums the
+  campaign's fragments. Nothing in the workflow reads either map. `nexp` is
+  what randoms need to match an `N_EPOCH` cut; `nflagged` is a diagnostic of
+  where the detector is bad, not a cut (it does not count exposures lost to
+  the defect veto).
 - **External masks are wired, on the tile side only (data runs).** `inputs.masks`
   is a third input root beside tiles and exposures, set per machine in the
   `machines:` table, exported as `$SP_INPUT_MASKS` and