diff --git a/astra.yaml b/astra.yaml index b58cdafed..e95afea59 100644 --- a/astra.yaml +++ b/astra.yaml @@ -231,7 +231,61 @@ analyses: catalogue's MASK__ 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_.hsp and nflagged_.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 diff --git a/tests/unit/test_campaign_lineage.py b/tests/unit/test_campaign_lineage.py index 492ec0ec8..77f937f28 100644 --- a/tests/unit/test_campaign_lineage.py +++ b/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") diff --git a/tests/unit/test_exp_maps.py b/tests/unit/test_exp_maps.py new file mode 100644 index 000000000..84501f600 --- /dev/null +++ b/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 diff --git a/tests/workflow/harness.py b/tests/workflow/harness.py index 1c8420f26..bfed3d553 100644 --- a/tests/workflow/harness.py +++ b/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 diff --git a/tests/workflow/params_pin.json b/tests/workflow/params_pin.json index 2284a7222..456530c00 100644 --- a/tests/workflow/params_pin.json +++ b/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", diff --git a/tests/workflow/test_dag.py b/tests/workflow/test_dag.py index b23493ecc..37d7af508 100644 --- a/tests/workflow/test_dag.py +++ b/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) ) diff --git a/universes/committed.yaml b/universes/committed.yaml index eac627dd1..265df3264 100644 --- a/universes/committed.yaml +++ b/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 diff --git a/workflow/CONTRACTS b/workflow/CONTRACTS index 9538e78cf..a940062d4 100644 --- a/workflow/CONTRACTS +++ b/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_.hdf5` and its -`patches/` group, `full_starcat_.hdf5`, and the `$run` in the +`patches/` group, `full_starcat_.hdf5`, `nexp_.hsp`, +`nflagged_.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. diff --git a/workflow/README.md b/workflow/README.md index 059546dab..d8aa395ed 100644 --- a/workflow/README.md +++ b/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_.hdf5 merge_final_cat.py ALL tiles' final_cat -> final_cat_.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_.hsp, nflagged_.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_.hsp` counts per sky pixel the + exposures whose CCD with a valid PSF model covers it, and + `nflagged_.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 diff --git a/workflow/Snakefile b/workflow/Snakefile index 84e816f67..5f1198265 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -465,6 +465,8 @@ FOREST_HASH = script_hash("build_forest.py") CLEAN_HASH = script_hash("clean_exposure.py") CLEAN_TILE_HASH = script_hash("clean_tile.py") PERSIST_HASH = script_hash("persist_exp.py") +EXP_MAPS_HASH = script_hash("exp_maps.py") +MERGE_MAPS_HASH = script_hash("merge_exposure_maps.py") MERGE_STAR_HASH = script_hash("merge_star_cat.py") # Three files, one trigger: the rule's script, the column extraction it calls, # and the parameter file that says which columns (path_hash argues why). @@ -844,13 +846,19 @@ def star_cat_inputs(): manifest and is in no set at all. Nothing short of rebuilding its chain from VOS recovers it; the merge reports how many exposures it found. """ - live, reclaimed = [], [] - for exp in psf_exposures(): - if not exp_store_reclaimed(exp): - live.append(prod_exp_manifest(exp, "exp_persist")) - elif Path(prod_exp_tar(exp)).exists(): - reclaimed.append(prod_exp_tar(exp)) - return live + reclaimed + return [p for exp in psf_exposures() + for p in durable_edge(exp, "exp_persist", prod_exp_tar(exp))] + + +def durable_edge(exp, stage, product): + """The DAG edge to a per-exposure product on the persistent root: the + stage's manifest while the scratch store is live, the product itself once + it is reclaimed (no rule declares it, so it is a leaf that builds nothing), + and nothing when a reclaimed store left no product. star_cat_inputs() + argues why a reclaimed exposure's manifest must never be named.""" + if not exp_store_reclaimed(exp): + return [prod_exp_manifest(exp, stage)] + return [product] if Path(product).exists() else [] @functools.lru_cache(maxsize=1) @@ -870,6 +878,49 @@ def star_cat_exposures(): or not exp_store_reclaimed(e)] +# --- the exposure maps ------------------------------------------------------ +# exp_maps rasterizes each exposure's valid-PSF CCDs and their instrument flags +# into one HealSparse fragment on the persistent root; exposure_maps sums the +# campaign's fragments into nexp_.hsp and nflagged_.hsp. Fitted-PSF +# runs only (psf_exposures()): the valid-PSF CCD set comes from exp_persist. +# The fragment follows exp_persist's tar exactly: requested for live stores, +# depended on through durable_edge() everywhere else. + +def prod_exp_maps(exp): + """The fragment exp_maps writes; not a declared output (see durable_edge).""" + return f"{prod_exp_dir(exp)}/maps/maps-{exp}.hsp" + + +def nexp_map(): + return f"{PRODUCTS_DIR}/nexp_{CAMPAIGN}.hsp" + + +def nflagged_map(): + return f"{PRODUCTS_DIR}/nflagged_{CAMPAIGN}.hsp" + + +@functools.lru_cache(maxsize=1) +def exposure_maps_inputs(): + return [p for exp in psf_exposures() + for p in durable_edge(exp, "exp_maps", prod_exp_maps(exp))] + + +def exposure_maps_exposures(): + """The exposure ids exposure_maps sums: its rerun fingerprint.""" + return [e for e in psf_exposures() + if Path(prod_exp_maps(e)).exists() or not exp_store_reclaimed(e)] + + +def exposure_maps_targets(): + """The per-exposure fragments to build, and the campaign maps when there is + at least one fragment to sum.""" + if not workflow.is_main_process: + return [] + inputs = exposure_maps_inputs() + live = [p for p in inputs if p.endswith("exp_maps.json")] + return live + ([nexp_map(), nflagged_map()] if inputs else []) + + # --- sizing the two merges (D4) --------------------------------------------- # MEASURED, not guessed, and measured as a SLOPE rather than a single number: # these are the only two rules whose one job's footprint grows with the whole @@ -1274,6 +1325,7 @@ rule all: [final_cat(t) for t in TILES_READY], persist_targets(), star_cat_targets(), + exposure_maps_targets(), final_cat_targets(), clean_targets(), clean_tile_targets(), diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk index 221d10928..12ad09a6f 100644 --- a/workflow/rules/exposure.smk +++ b/workflow/rules/exposure.smk @@ -1,6 +1,6 @@ """Exposure chain — per exposure, keyed by exp base id (dedup is structural). - exp_get_images -> exp_split -> exp_psf -> exp_persist + exp_get_images -> exp_split -> exp_psf -> exp_persist -> exp_maps Each in the exposure's own sharded work dir, chained by manifests; every config reads fixed ``$SP_RUN/output/run_sp_exp_*`` INPUT_DIRs, so nothing resolves a @@ -19,7 +19,11 @@ opt-in, see ``star_selection.setools``), and ``make_cat`` writes the per-band star catalogue, or a network fetch — hence no ``star_catalogue`` / ``exp_star_cat`` here, and no ``exp_mask``. -``exp_persist`` is the one rule here that writes to the PERSISTENT root: it +``exp_maps`` exports that flag image, with the CCD footprints, into a +HealSparse fragment so the survey footprint can subtract it; nothing in this +workflow reads it back. + +``exp_persist`` and ``exp_maps`` write to the PERSISTENT root. ``exp_persist`` packs the PSF products named by `persist_exp:` into one tar per exposure off /scratch before the purge (or clean_exposure) can take them. It is a separate rule from exp_psf precisely so that editing that list costs a re-pack and not a @@ -178,6 +182,38 @@ rule exp_persist: " {params.patterns}" +# --- the exposure's footprint and defects ----------------------------------- +# One HealSparse fragment per exposure on the persistent root: the sky pixels +# its valid-PSF CCDs cover, and how many flagged CCD pixels fall in each +# (workflow/scripts/exp_maps.py). It reads the split images' headers and flag +# splits from the scratch store, so clean_exposure waits for it exactly as for +# exp_persist; its one declared input is exp_persist's manifest, which names +# the valid-PSF CCDs. ~8 s per exposure. +# @sc [decision:masking.defect_map_from_flags] +# @sc [decision:masking.nexp_map_valid_psf_ccds] +rule exp_maps: + input: + persist = lambda wc: prod_exp_manifest(wc.exp, "exp_persist") + output: + manifest = f"{PROD_EXP_DIR}/manifests/exp_maps.json" + params: + exp_dir = lambda wc: exp_dir(wc.exp), + fragment = lambda wc: prod_exp_maps(wc.exp), + script_hash = EXP_MAPS_HASH + threads: 1 + retries: 2 + resources: + # 0.15 GB typical; 1.5 GB for a fully flagged CCD (9.4M pixels) + mem_mb = 3000, + runtime = 20 + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/exp_maps.py" + " --exp-dir '{params.exp_dir}' --exp {wildcards.exp}" + " --persist-manifest {input.persist}" + " --fragment '{params.fragment}' --manifest {output.manifest}" + + # --- reclamation (D5) ------------------------------------------------------- # The one exception to "no reclamation in this file": clean_exposure OWNS # exposure-level deletion, and it is a real job, not temp() bookkeeping, because @@ -224,11 +260,11 @@ rule clean_exposure: # changes exp_persist's params, the manifest reruns, and it sits behind # exp_psf's manifest, which went with the store, so snakemake rebuilds # the exposure from VOS. - lambda wc: ([] if not PERSISTS_PSF - else [prod_exp_manifest(wc.exp, "exp_persist")] - if not exp_store_reclaimed(wc.exp) - else [prod_exp_tar(wc.exp)] - if Path(prod_exp_tar(wc.exp)).exists() else []) + lambda wc: (durable_edge(wc.exp, "exp_persist", prod_exp_tar(wc.exp)) + if PERSISTS_PSF else []), + # The exposure maps' fragment, by the same rule. + lambda wc: (durable_edge(wc.exp, "exp_maps", prod_exp_maps(wc.exp)) + if PERSISTS_PSF else []) output: tombstone = f"{EXP_DIR}/cleaned.json" params: @@ -325,3 +361,40 @@ rule star_cat_merge: " --output {output.star_cat}" " --campaign '{params.campaign}'" " --snapshot-json '{params.snapshot}'" + + +# --- the campaign's exposure maps ------------------------------------------- +# ONE job per campaign: every fragment of the campaign's exposures, summed into +# /nexp_.hsp (exposures with a valid PSF model per sky pixel) +# and nflagged_.hsp (their flagged CCD pixels per sky pixel). Rebuilt +# whole when the exposure set or a fragment changes; inputs, fingerprint and +# the job's own rediscovery of the set follow star_cat_merge. Memory is the two +# maps, 3 MiB per nside-128 coverage pixel the campaign touches: ~13 per +# exposure, capped at the ~23k of the UNIONS footprint (~70 GB at DR6). +# @sc [decision:masking.defect_map_from_flags] +# @sc [decision:masking.nexp_map_valid_psf_ccds] +rule exposure_maps: + input: + lambda wc: exposure_maps_inputs() + output: + nexp = nexp_map(), + nflagged = nflagged_map() + params: + products_dir = str(PRODUCTS_DIR), + tile_list = str(config["tile_list"]), + index_db = str(INDEX_DB), + inputs = unit_fingerprint(exposure_maps_exposures()), + script_hash = MERGE_MAPS_HASH + threads: 1 + resources: + mem_mb = lambda wc, attempt: capped_mem(attempt * ( + 1000 + 3.2 * min(13 * len(exposure_maps_exposures()), 23_000)), + "exposure_maps"), + runtime = lambda wc, attempt: attempt * ( + 30 + len(exposure_maps_exposures()) // 30) + shell: + "set -euo pipefail\n" + f"python {SCRIPTS}/merge_exposure_maps.py" + " --products-dir '{params.products_dir}'" + " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" + " --nexp {output.nexp} --nflagged {output.nflagged}" diff --git a/workflow/scripts/exp_maps.py b/workflow/scripts/exp_maps.py new file mode 100644 index 000000000..d814856df --- /dev/null +++ b/workflow/scripts/exp_maps.py @@ -0,0 +1,148 @@ +"""Rasterize ONE exposure's footprint and instrument defects into a HealSparse +fragment: the shell of the ``exp_maps`` rule. + +The fragment is a uint8 map at the UNIONS mask ladder's resolution (nside +131072 over coverage 128, 1.6" pixels, ~74 MegaCam pixels each), so it aligns +pixel for pixel with the sky masks. A sky pixel holds + + 0 not covered by a CCD of this exposure with a PSF model + 1 + n covered, and ``n`` flagged CCD pixels have their centre in it + +``merge_exposure_maps.py`` sums the fragments of a campaign into the +exposure-count map (``nexp``) and the flagged-pixel map (``nflagged``). The +CCDs of one exposure do not overlap, so the first sum counts exposures. + +Flagged pixels are counted rather than OR-ed into a boolean mask because most +of them are thin: bad columns and cosmic rays a single CCD pixel wide. Any-touch +rasterization widens a one-pixel column to a 1.6" strip; on real exposures it +masks ~7x more sky than the flagged pixels cover. The count keeps the area +(``nflagged / 74`` is the fraction of one exposure lost there), and +``nflagged > 0`` is still the any-touch mask. + +Inputs, from the exposure's scratch store and its exp_persist manifest: + +* CCDs with a PSF model: the ``validation_psf--.fits`` members of + ``exp_persist.json``. psfex_interp writes that file only when the fit + succeeds, so these are the CCDs that can contribute a shape. +* WCS and imaging area: the header of ``image--.fits``. The split is + 2112 x 4644 pixels with a flagged overscan border (flag value 3); only + ``DATASEC`` (2048 x 4612) sees the sky. The flag split carries no WCS. +* defects: ``flag--.fits``, any nonzero bit inside ``DATASEC``. + +A sky pixel is covered when its centre lies inside the polygon through the +imaging area's outer corners. +""" + +import argparse +import filecmp +import json +import warnings +from fnmatch import fnmatch +from pathlib import Path + +import healsparse as hsp +import hpgeom as hpg +import numpy as np +from astropy.io import fits +from astropy.wcs import WCS + +import persist_exp + +NSIDE = 131072 +NSIDE_COVERAGE = 128 +SPLIT_DIR = "output/run_sp_exp_Sp/split_exp_runner/output" +PSF_PATTERN = persist_exp.resolve(persist_exp.ALWAYS) + + +def valid_ccds(manifest: Path, exp: str) -> list: + """CCD numbers with a PSF model, from the exp_persist manifest's members.""" + prefix = PSF_PATTERN.split("*")[0] + f"{exp}-" + names = [f["name"] for f in json.loads(manifest.read_text())["files"]] + return sorted(int(n[len(prefix):].removesuffix(".fits")) for n in names + if fnmatch(n, PSF_PATTERN) and n.startswith(prefix)) + + +def ccd_wcs(image: Path): + """The CCD's WCS and imaging area ``(x0, x1, y0, y1)``, 1-based inclusive, + from ``DATASEC`` (or the whole array without one).""" + header = fits.getheader(image) + if "DATASEC" in header: + x, y = header["DATASEC"].strip("[]").split(",") + bounds = [int(v) for v in x.split(":") + y.split(":")] + else: + bounds = [1, header["NAXIS1"], 1, header["NAXIS2"]] + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + return WCS(header), bounds + + +def coverage_pixels(wcs, bounds) -> np.ndarray: + """Sky pixels (NEST) whose centres lie inside the imaging area.""" + x0, x1, y0, y1 = bounds + ra, dec = wcs.all_pix2world([x0 - 0.5, x1 + 0.5, x1 + 0.5, x0 - 0.5], + [y0 - 0.5, y0 - 0.5, y1 + 0.5, y1 + 0.5], 1) + return hpg.query_polygon(NSIDE, ra, dec, nest=True) + + +def flagged_counts(wcs, flags: np.ndarray, bounds): + """Sky pixels (NEST) holding flagged imaging-area pixel centres, and how + many each holds.""" + x0, x1, y0, y1 = bounds + rows, cols = np.nonzero(flags[y0 - 1:y1, x0 - 1:x1]) + ra, dec = wcs.all_pix2world(cols + x0, rows + y0, 1) + return np.unique(hpg.angle_to_pixel(NSIDE, ra, dec, nest=True), + return_counts=True) + + +def write_stable(tmp: Path, dest: Path) -> None: + """Replace ``dest`` with ``tmp`` unless the bytes match, keeping the mtime + of an unchanged product (mtime is a snakemake rerun trigger).""" + if dest.exists() and filecmp.cmp(tmp, dest, shallow=False): + tmp.unlink() + else: + tmp.replace(dest) + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--exp-dir", required=True, type=Path) + parser.add_argument("--exp", required=True) + parser.add_argument("--persist-manifest", required=True, type=Path) + parser.add_argument("--fragment", required=True, type=Path) + parser.add_argument("--manifest", required=True, type=Path) + args = parser.parse_args() + + split = args.exp_dir / SPLIT_DIR + fragment = hsp.HealSparseMap.make_empty(NSIDE_COVERAGE, NSIDE, np.uint8) + ccds = {} + for ccd in valid_ccds(args.persist_manifest, args.exp): + wcs, bounds = ccd_wcs(split / f"image-{args.exp}-{ccd}.fits") + covered = coverage_pixels(wcs, bounds) + pixels, counts = flagged_counts( + wcs, fits.getdata(split / f"flag-{args.exp}-{ccd}.fits"), bounds) + keep = np.isin(pixels, covered, assume_unique=True) + fragment[covered] = np.ones(covered.size, np.uint8) + fragment[pixels[keep]] = (1 + np.minimum(counts[keep], 254)).astype(np.uint8) + ccds[ccd] = {"covered": int(covered.size), + "flagged_pixels": int(counts.sum())} + + args.fragment.parent.mkdir(parents=True, exist_ok=True) + tmp = args.fragment.with_name(f"{args.fragment.stem}.tmp.hsp") + fragment.write(str(tmp), clobber=True) + write_stable(tmp, args.fragment) + + body = {"stage": "exp_maps", "level": "exp", "unit": args.exp, + "status": "complete", "fragment": str(args.fragment), + "nside": NSIDE, "ccds": ccds} + args.manifest.parent.mkdir(parents=True, exist_ok=True) + tmp = args.manifest.with_name(args.manifest.name + ".tmp") + tmp.write_text(json.dumps(body, indent=1, sort_keys=True) + "\n") + write_stable(tmp, args.manifest) + covered = sum(c["covered"] for c in ccds.values()) + flagged = sum(c["flagged_pixels"] for c in ccds.values()) + print(f"[exp_maps] {args.exp}: {len(ccds)} CCD(s) with a PSF model, " + f"{covered} sky pixel(s), {flagged} flagged CCD pixel(s)") + + +if __name__ == "__main__": + main() diff --git a/workflow/scripts/merge_exposure_maps.py b/workflow/scripts/merge_exposure_maps.py new file mode 100644 index 000000000..d5a79caee --- /dev/null +++ b/workflow/scripts/merge_exposure_maps.py @@ -0,0 +1,84 @@ +"""Sum a campaign's per-exposure fragments into its two HealSparse maps: the +shell of the ``exposure_maps`` rule. + +* ``nexp_.hsp`` (uint8): per sky pixel, the number of exposures whose CCD + with a valid PSF model covers it, the count behind sp_validation's + ``npoint`` cut. +* ``nflagged_.hsp`` (uint16): per sky pixel, the flagged CCD pixels of + those exposures whose centres fall in it. A sky pixel holds ~74 MegaCam + pixels, so ``nflagged / 74`` is the number of exposures' worth of area lost + there, and ``nflagged > 0`` the any-touch defect mask. + +Both are 0 where nothing is counted. The campaign's exposures come from the tile list +and the index (``build_index``), as for the star-catalogue merge; an exposure +without a fragment (its store reclaimed before ``exp_maps`` existed) is left +out and counted. The maps are rebuilt from every fragment on each run. +""" + +import argparse +from pathlib import Path + +import healsparse as hsp +import numpy as np + +import build_index +from exp_maps import NSIDE, NSIDE_COVERAGE + + +def fragment_path(products_dir: Path, exp: str) -> Path: + """Where ``exp_maps`` writes an exposure's fragment (Snakefile: + ``prod_exp_maps``).""" + return products_dir / "exp" / exp[:2] / exp / "maps" / f"maps-{exp}.hsp" + + +def add(acc, pixels, values): + """Add ``values`` to ``acc`` at ``pixels``.""" + acc[pixels] = acc.get_values_pix(pixels) + values.astype(acc.dtype) + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--products-dir", required=True, type=Path) + parser.add_argument("--tile-list", required=True, type=Path) + parser.add_argument("--index-db", required=True, type=Path) + parser.add_argument("--nexp", required=True, type=Path) + parser.add_argument("--nflagged", required=True, type=Path) + args = parser.parse_args() + + exposures = build_index.campaign_exposures(args.tile_list, args.index_db) + nexp = hsp.HealSparseMap.make_empty(NSIDE_COVERAGE, NSIDE, np.uint8) + nflagged = hsp.HealSparseMap.make_empty(NSIDE_COVERAGE, NSIDE, np.uint16) + missing = [] + for i, exp in enumerate(exposures): + path = fragment_path(args.products_dir, exp) + if not path.exists(): + missing.append(exp) + continue + fragment = hsp.HealSparseMap.read(str(path)) + pixels = fragment.valid_pixels + values = fragment.get_values_pix(pixels) + add(nexp, pixels, np.ones_like(values)) + flagged = values > 1 + add(nflagged, pixels[flagged], values[flagged] - 1) + if i % 500 == 0: + print(f"[exposure_maps] {i + 1}/{len(exposures)} {exp}", flush=True) + + used = len(exposures) - len(missing) + if not used: + raise SystemExit("merge_exposure_maps: no exposure of the campaign has " + "a fragment") + for acc, out in ((nexp, args.nexp), (nflagged, args.nflagged)): + out.parent.mkdir(parents=True, exist_ok=True) + tmp = out.with_name(f"{out.stem}.tmp.hsp") + acc.write(str(tmp), clobber=True) + tmp.replace(out) + print(f"[exposure_maps] {used} exposure(s) -> {args.nexp.name}, " + f"{args.nflagged.name}") + if missing: + print(f"[exposure_maps] WARNING: {len(missing)} exposure(s) in the " + f"campaign have no fragment (reclaimed before exp_maps ran) and " + f"are not counted, e.g. {missing[:5]}") + + +if __name__ == "__main__": + main()