Consume external healsparse masks: rasterize to pipeline flag images + per-band catalogue flags

[cailmdaley]

Closes #846. Part of the masking-unification epic #845; design reviewed in person by Martin + Cail 2026-07-16 (starting from the proposal on #839); Axel's architecture comments from the Slack thread (2026-07-16) are folded in below.

Motivation

UNIONS has converged on a single mask system: the external healsparse masks produced by the photometry-side mask effort (PhotoPipe ugriz footprint + bright stars, MaxiMask r, manual 2MASS galaxy masks — merged and converted to healsparse at nside 131072 ≈ 1.6″, conservatively fattened). These are the masks that set the catalogue flags and feed cosmology. ShapePipe should consume them rather than generate its own THELI-vintage masks (GSC-2.3 star halos/spikes rasterized by WeightWatcher, NGC/Messier region files, live Vizier queries).

What this buys:

• One mask product across photometry, shape measurement, and cosmology.
• Controlled, consistent flags upstream of PSF selection — the old GSC path misses stars (NaN-magnitude bug Extend size of bright star masks [BUG] #662, close pairs absent from GSC by design); the external masks give a complete, configurable flag source (which bits actually cut star candidates is policy, §3).
• Removes the two flakiest external dependencies from the critical path: the WeightWatcher executable and the live Vizier/GSC-2.3 query.

The interface

Everything downstream of the mask module consumes a single artifact: the per-tile / per-exposure-CCD pixel flag image pipeline_flag<num>.fits. SExtractor turns it into IMAFLAGS_ISO; setools star selection cuts IMAFLAGS_ISO == 0; vignetmaker extracts flag stamps; ngmix zeroes weights on flagged pixels. Nothing downstream cares how the flag image was made — so the swap is local: read the healsparse map, rasterize it onto the image pixel grid through the WCS, write the same flag FITS.

Design

1. New module: mask_ext_runner

A sibling of mask_runner (the old module stays untouched and config-selectable — see §5). Per input image:

1. Read the image header → astropy.wcs.WCS (same machinery as Mask._set_image_coordinates).
2. Evaluate the pixel grid → (RA, Dec) in row chunks (bounded memory; a tile is 10⁸ pixels, an exposure CCD 10⁷).
3. healsparse.HealSparseMap.get_values_pos(ra, dec) → per-pixel mask bits.
4. Map healsparse bits → ShapePipe flag values through a config-driven bit→flag mapping, producing an int16 flag image.
5. Optionally sum in the instrument flag image (the existing USE_EXT_FLAG pattern — needed for exposures, where saturation/bleeding/bad columns arrive with the data).
6. Write pipeline_flag<num>.fits with the image WCS in the header.

Config sketch:

[MASK_EXT_RUNNER]
INPUT_DIR = last:split_exp_runner        # or tile input for tile runs
NUMBERING_SCHEME = -0000000-0
MASK_PATH = /path/to/mask_r.hsp          # healsparse file; band-agnostic
BIT_FLAG_MAP = 1:1                       # healsparse bit -> output flag value(s)
USE_EXT_FLAG = True                      # exposures: sum instrument flags
PREFIX = pipeline

The same module serves tiles and exposure CCDs — only the input image (hence WCS) differs, exactly as the current mask module runs twice. On the content of exposure-level masks, Axel (Slack, 2026-07-16) answers the previously-open question: pixel masks belong primarily at the single-exposure level (where PSF and shape measurement happen), composed of Stephen's bad pixels + cosmic rays, updated star spikes/halos, and bright nearby objects; the healsparse maps can be created per exposure and combined (Mike's alternative: reproject the tile-built maps). Either resolution leaves this module unchanged — it rasterizes whatever healsparse file MASK_PATH points at through each image's WCS.

2. Strictly mask-file-agnostic

Mask paths and bit meanings live only in config. The mask products are actively evolving (DR6 three-set merge in progress; Axel also proposes moving the catalogue-level product to healsparse polygons — effectively infinite resolution, straight from the DS9 regions); we start with the 2025 all-band masks and swap files later at zero code cost.

3. PSF star selection: consumer-aware bit policy

setools cuts IMAFLAGS_ISO == 0, so any bit mapped into the exposure-level flag image rejects PSF star candidates. Per the 2026-07-21 mask-force telecon (notes), star-halo bits (1 faint, 2 bright) must not reject PSF stars — halo masks flag, selection happens at the catalogue level — while defect bits (e.g. MaxiMask streaks) may cut. The BIT_FLAG_MAP config expresses this directly: map only the bits that should reach IMAFLAGS_ISO, and carry the rest as catalogue flags. Provisional on the surviving-PSF-star overlay check in #850.

4. Per-band catalogue flags via make_cat

Shape measurement only needs the r-band mask rasterized. The other bands are catalogue metadata, and healsparse maps are query-shaped products — so make_cat gains an optional per-band healsparse lookup at each object's (RA, Dec) (the existing add_col pattern), writing e.g. MASK_u, MASK_g, … columns into the final catalogue. Off-footprint objects carry the map's sentinel (False for the boolean maps) in these columns, while the rasterized pixel mask maps off-footprint to a configurable OFF_MAP_FLAG (default 0) — a deliberate split (pixel flags feed IMAFLAGS_ISO; catalogue columns are metadata). This is the ShapePipe end of UNIONS-WL/spherex#38, and mirrors how the healsparse masks already set catalogue flags on the photometry side.

[MAKE_CAT_RUNNER]
...
MASK_EXT_PATHS = u:/path/mask_u.hsp, g:/path/mask_g.hsp, ...

5. The old path stays, and gets compared

Nothing is deleted. mask_runner (WeightWatcher/GSC/region files) remains fully functional and config-selectable. The comparison gate is #850: old-vs-new on real data (the Nibi P3 tiles/exposures) — visual overlays on real pixels, a component-provenance audit (where saturation/bleeding/bad columns/cosmic rays/halos/spikes/galaxies enter under each path), masked-area fractions, IMAFLAGS_ISO populations, PSF star selections and PSF model quality. Removal of the old path is gated on #850. Axel's further proposal — drop the tile-level mask input to SExtractor entirely (redundant if flags are applied at catalogue level) — is evaluated there too, not here. Mask-size validation via tangential shear around bright stars is #851.

Out of scope

• Deleting the internal masking code (gated on Old-vs-new mask comparison on real data (tiles + exposures) #850).
• Generating masks (MaxiMask, star-halo scaling, the DR6 merge, healsparse polygons) — that is mask-force/PhotoPipe territory; we consume the product.
• Exposure-level mask content (per-exposure generation + combine vs reprojection of tile maps, Exposure-level star/defect masking strategy #839/Masking unification: external healsparse masks throughout ShapePipe #845): either resolution keeps this PR's artifact contract unchanged.

Data

All known mask products are boolean healsparse maps (dtype bool, True = masked, valid pixels only where masked), nside 131072, coverage nside 128:

• The 2025 all-band masks (arc:home/mhudson/masks/, now mirrored to Nibi at ~/projects/def-mjhudson/cdaley/masks/) come as one file per bit, mask_ugriz_nside131072_n<bit>.hsp, produced by Mike's Healsparse_mask_pool.py from Ludo's tile-by-tile pixel FITS masks. Documented reason bits (script header): 1 faint star halo / 2 bright star halo / 4 bright stars / 8 manual mask / 1024 maximask; 512 = out-of-bounds. Empirical check (P3 region, 200k random points): OR{1,2,4,8,64,1024} reproduces mask_r_nside131072.hsp exactly — the 2025 series is a per-bit decomposition of the same masking (which also says the r masking is unchanged Oct 2024 → May 2025 there). The earlier per-band reading (16_u/32_g/64_r/128_i/256_z/2048_z2 as band masks) is contradicted: 64 is a reason bit inside Mike's default bitmask (0b10001001111), and 16/32/128/256 behave like per-band coverage flags (n128 fires on 100% of the P3 region, where i-band is absent; n256 on 71%; n2048 is empty). Definitive bit semantics need Mike — until then, per-band MASK_<band> catalogue columns (§4) should not be wired to the n16/n32/… files.
• The r-band map on candide (/n17data/UNIONS/WL/masks/mask_r_nside131072.hsp, 932 MB, RICE/bit-packed) is the same bool flavour.

For bool maps the config is BIT_FLAG_MAP = 1:<flag>; any other bit would silently select nothing (True & 64 == 0), so the module fails loudly on that misconfiguration. Integer bit-flag maps (a single map with per-band/per-reason bits) don't currently exist, but the bit→flag mapping supports them if the DR6 merge produces one.

Implementation plan

1. ✅ Deps + skeleton — healsparse, hpgeom in pyproject.toml/uv.lock (also declared on the coverage branch — no code overlap, only dep-file merge surface).
2. ✅ Rasterizer + unit tests — mask_ext_package + mask_ext_runner; 7 tests on synthetic maps (bit mapping, chunk-seam independence, RA wrap, off-map, ext-flag summing, WCS round-trip).
3. ✅ Pipeline wiring — example tile + exposure configs added (config_tile_MaExt.ini, config_exp_MaExt.ini); shapepipe_run drives mask_ext_runner end-to-end (smoke below).
4. ✅ make_cat per-band columns + 3 tests.
5. ✅ Real-data smoke, tile path — the real r-band mask (mask_r_nside131072.hsp) rasterized onto the exact CFIS.233.293 tile pixel grid (10000², TAN, 0.187″/px) via shapepipe_run: 38.0% masked, matching a direct healsparse query of the region (37.0%); morphology is unmistakably real (polygonal bright-star halos, diffraction-spike crosses, a bleed column, defect speckle). Overlay against actual image pixels (do masks land on the stars?) folds into Old-vs-new mask comparison on real data (tiles + exposures) #850.
6. ✅ Real-data smoke, exposure path (USE_EXT_FLAG summing) — the same mask rasterized through the WCS of P3 split CCD 2086317-0 (2112×4644, from the Set up ShapePipe on Nibi and run end-to-end on a first batch of P3 tiles #808 run outputs on Nibi) with the instrument flag summed in, via shapepipe_run on the develop-runtime container: external mask lands on 20.4% of the CCD as bit 512; instrument bits (1/2/3/8) preserved exactly; combined values (513/515/523) confirm clean summing; morphology shows spike crosses, a bright-star polygon, and faint-star speckle against the instrument panel's border + bad columns.

Implementation was independently reviewed against this description (fresh-context adversarial pass; tests rerun in the container): no correctness findings.

— Fable, on behalf of Cail

[cailmdaley]

Exposure-path smoke: PASS (Nibi, real data). mask_ext_runner run standalone on split CCD 2114049-28 (MegaCam 2112×4644, the densest instrument-flag CCD of the exposure: 4.9% nonzero, values {1,2,3,8,11} — bad columns/saturation), with MASK_PATH = mask_ugriz_nside131072_n2.hsp (bright-star halos, bool), BIT_FLAG_MAP = 1:16, USE_EXT_FLAG = True.

• Healsparse component vs an independent full-resolution rasterization (separate code, same WCS): 100.000000% pixel agreement, 0/9,808,128 mismatches.
• Instrument flag summed through everywhere; output − inst_flag ∈ {0, 16} — clean additive decomposition, no bit collisions.
• Fractions: healsparse 11.84%, instrument 4.94%, overlap 0.92%, total 15.86% (union checks out).
• ~13 s for the 9.8M-px CCD, single chunk, WCS carried into the output header.

With the tile smoke above, both halves of the interface (§1) are validated on real data end-to-end. Old-vs-new comparison continues in #850.

— Fable, on behalf of Cail

[cailmdaley]

Policy update from the 2026-07-21 mask-force telecon (details on #845): the rasterizer design here is unchanged, but the BIT_FLAG_MAP policy is now consumer-aware. Star-halo bits (1, 2) must not enter the flags that setools' IMAFLAGS_ISO==0 cut sees for PSF star selection — halo masks flag, they don't reject. Defect bits (e.g. MaxiMask/streaks) may still cut. The issue body's PSF-star-selection paragraph has been updated accordingly; provisional on the surviving-star overlay check (#850).

— Claude (fable) on behalf of Cail