diff --git a/astra.yaml b/astra.yaml index 111d4e538..36a914c40 100644 --- a/astra.yaml +++ b/astra.yaml @@ -235,14 +235,34 @@ analyses: 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. + npoint_.hsp, ntaken_.hsp, nexp_.hsp and + nflagged_.hsp at the mask ladder's resolution: filled main-array + pointings, all CCD DATASEC exposures, validated-PSF CCD exposures, + and flagged CCD pixel centres of validated CCDs per sky pixel. + Products; nothing in the workflow reads them. inputs: [exposure_flags] - decisions: [nexp_map_valid_psf_ccds, defect_map_from_flags] + decisions: [pointing_square, nexp_map_valid_psf_ccds, defect_map_from_flags] decisions: + pointing_square: + label: Filled main-array pointing coverage, excluding ears + rationale: >- + npoint is deliberately smooth pointing coverage, the v2 counterpart + of v1.X's npoint3 >= 3 cut. CCD tiling, gaps, ears and dithers produce + patterned imaging-area counts instead. The convex hull of DATASEC + outer corners (half a pixel beyond pixel centres) of present CCDs + 0–35 is computed in a tangent plane about the exposure centre, mapped + back to RA/Dec, and rasterized by sky-pixel centre inclusion. CCDs + 36–39 are ears outside the main 4x9 array and are excluded. Every + fragment contributes regardless of PSF success; gaps are filled. + default: main_array_hull + options: + main_array_hull: + label: Convex hull of present CCDs 0–35, gaps filled and ears excluded + ccd_tiling: + label: Individual CCD DATASEC footprints including gaps and ears + excluded: true + excluded_reason: >- + Measures taken imaging area (ntaken), not smooth pointing coverage. nexp_map_valid_psf_ccds: label: The exposure count covers only CCDs with a valid PSF model rationale: >- @@ -250,13 +270,17 @@ analyses: 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. + exposures that can contribute a PSF-corrected epoch. ntaken counts + all present CCD DATASEC footprints, ears included, without the PSF + requirement; ntaken - nexp measures epochs lost to CCD-level PSF + failure at each pixel. Both use centre-in-polygon inclusion and count + each exposure once even at overlapping footprint edges. 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 + label: Use all CCD coverage for nexp rather than separately in ntaken excluded: true excluded_reason: >- Counts epochs that cannot contribute a shape, overstating depth diff --git a/tests/unit/test_exp_maps.py b/tests/unit/test_exp_maps.py index 84501f600..1176cae0d 100644 --- a/tests/unit/test_exp_maps.py +++ b/tests/unit/test_exp_maps.py @@ -1,7 +1,7 @@ """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 +The fixture is a synthetic CCD grid 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 @@ -20,7 +20,9 @@ import pytest pytestmark = [pytest.mark.unions, - pytest.mark.decision("masking.defect_map_from_flags")] + pytest.mark.decision("masking.defect_map_from_flags"), + pytest.mark.decision("masking.nexp_map_valid_psf_ccds"), + pytest.mark.decision("masking.pointing_square")] SCRIPTS = Path(__file__).resolve().parents[2] / "workflow" / "scripts" EXP = "2079612" @@ -132,8 +134,8 @@ def test_fragment_counts_every_flagged_pixel(maps, tmp_path, monkeypatch): 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} + values = (fragment.get_values_pix(pixels) >> maps.FLAG_SHIFT).astype(int) + got = {p: v for p, v in zip(pixels, values) if v > 0} # 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) @@ -144,20 +146,22 @@ def test_fragment_counts_every_flagged_pixel(maps, tmp_path, monkeypatch): 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.""" + """Validated CCD coverage excludes overscan; unvalidated CCDs are taken.""" 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() + assert (fragment.get_values_pix(inside) & 4 != 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() + assert (fragment.get_values_pix(overscan) & 2 == 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"] + values = fragment.get_values_pix(other) + assert (values & 2 != 0).all() + assert (values & 4 == 0).all() + assert (values >> maps.FLAG_SHIFT == 0).all() + assert list(record["ccds"]) == ["0", "1"] def test_merge_counts_exposures_and_flagged_pixels(maps, tmp_path, @@ -172,19 +176,21 @@ def test_merge_counts_exposures_and_flagged_pixels(maps, tmp_path, 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) + maps.FRAGMENT_DTYPE) + values = np.full(pix.size, 7, maps.FRAGMENT_DTYPE) + values[500:510] |= 7 << maps.FLAG_SHIFT + frag[pix] = values 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")} + out = {k: tmp_path / f"{k}.hsp" for k in + ("npoint", "ntaken", "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"])]) + *[arg for k, path in out.items() for arg in (f"--{k}", str(path))]]) merge.main() nexp = hsp.HealSparseMap.read(str(out["nexp"])) @@ -195,3 +201,117 @@ def test_merge_counts_exposures_and_flagged_pixels(maps, tmp_path, 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 + for name in ("npoint", "ntaken"): + count = hsp.HealSparseMap.read(str(out[name])) + assert np.array_equal(count.get_values_pix(a), nexp.get_values_pix(a)) + + +def test_pointing_fills_gaps_excludes_ears_and_matches_mef(maps, tmp_path, + monkeypatch): + """A main-array grid has gaps, an unvalidated CCD, and an external ear.""" + import healsparse as hsp + import hpgeom as hpg + from astropy.io import fits + + # Arrange centres in the rotated WCS's own pixel plane; 80-pixel gaps. + base = _wcs((359.99, 31.7)) + offsets = [(0, 0), (NX + 80, 0), (0, NY + 80), + (NX + 80, NY + 80), (2 * (NX + 80), 0)] + crvals = [base.all_pix2world(base.wcs.crpix[0] + x, + base.wcs.crpix[1] + y, 1) for x, y in offsets] + ccds = [0, 1, 9, 10, 36] + fragment, _ = _run(maps, tmp_path, monkeypatch, + ccds=ccds, crvals=crvals, valid=[0, 9, 10, 36]) + ra, dec = base.all_pix2world(base.wcs.crpix[0] + (NX + 80) / 2, + base.wcs.crpix[1], 1) + gap = hpg.angle_to_pixel(maps.NSIDE, ra, dec) + assert fragment[gap] & maps.COVER_MASK == 1 + ear = _sky_pixels(maps, crvals[-1], NX // 2, NY // 2) + assert fragment[ear] & maps.COVER_MASK == 6 + unvalidated = _sky_pixels(maps, crvals[1], NX // 2, NY // 2) + assert fragment[unvalidated] & maps.COVER_MASK == 3 + assert fragment[unvalidated] >> maps.FLAG_SHIFT == 0 + + split = tmp_path / "store" / maps.SPLIT_DIR + images, flags = [fits.PrimaryHDU()], [fits.PrimaryHDU()] + for ccd in range(37): + if ccd in ccds: + header = fits.getheader(split / f"image-{EXP}-{ccd}.fits") + images.append(fits.ImageHDU(data=np.zeros(_split(_flags()).shape), + header=header)) + # Delivered floating flags follow split_exp's int16 conversion. + flags.append(fits.ImageHDU( + data=_split(_flags()).astype(np.float32) + 0.25)) + else: + # Raw MEFs index every CCD; put absent grid CCDs at an existing + # unvalidated footprint so coverage and flags remain unchanged. + header = fits.getheader(split / f"image-{EXP}-1.fits") + images.append(fits.ImageHDU(data=np.zeros(_split(_flags()).shape), + header=header)) + flags.append(fits.ImageHDU( + data=_split(_flags()).astype(np.float32) + 0.25)) + image_path, flag_path = tmp_path / "image.fits", tmp_path / "flag.fits" + fits.HDUList(images).writeto(image_path) + fits.HDUList(flags).writeto(flag_path) + raw_path = tmp_path / "raw.hsp" + monkeypatch.setattr(sys, "argv", [ + "exp_maps.py", "--image-mef", str(image_path), + "--flag-mef", str(flag_path), + "--exp", EXP, "--persist-manifest", str(tmp_path / "exp_persist.json"), + "--fragment", str(raw_path), "--manifest", str(tmp_path / "raw.json")]) + maps.main() + raw = hsp.HealSparseMap.read(str(raw_path)) + assert np.array_equal(raw.valid_pixels, fragment.valid_pixels) + assert np.array_equal(raw[raw.valid_pixels], + fragment[fragment.valid_pixels]) + + # Pointing and taken coverage do not depend on any PSF succeeding. + manifest = tmp_path / "exp_persist.json" + manifest.write_text(json.dumps({"files": []})) + maps.main() + empty_psf = hsp.HealSparseMap.read(str(raw_path)) + values = empty_psf[raw.valid_pixels] + assert np.array_equal(values & 3, raw[raw.valid_pixels] & 3) + assert (values & 4 == 0).all() + assert (values >> maps.FLAG_SHIFT == 0).all() + + +def test_fragment_clips_flag_counts_and_preserves_mtime(maps, tmp_path, + monkeypatch): + """Flag counts saturate at 255; identical writes preserve mtimes.""" + crval = (150.3, 31.7) + pixel = _sky_pixels(maps, crval, NX // 2, NY // 2) + monkeypatch.setattr(maps, "flagged_counts", lambda *_: ( + np.atleast_1d(pixel), np.array([400]))) + fragment, _ = _run(maps, tmp_path, monkeypatch, + ccds=[0], crvals=[crval], valid=[0]) + assert fragment[pixel] >> maps.FLAG_SHIFT == 255 + assert fragment[pixel] & maps.COVER_MASK == 7 + paths = [tmp_path / "maps" / f"maps-{EXP}.hsp", tmp_path / "exp_maps.json"] + mtimes = [p.stat().st_mtime_ns for p in paths] + maps.main() + assert [p.stat().st_mtime_ns for p in paths] == mtimes + + +def test_merge_refuses_plain_uint8(maps, tmp_path, monkeypatch): + """Legacy fragments require an explicit exposure-specific raw backfill.""" + import healsparse as hsp + + merge = _load("merge_exposure_maps") + path = merge.fragment_path(tmp_path, EXP) + path.parent.mkdir(parents=True) + old = hsp.HealSparseMap.make_empty( + maps.NSIDE_COVERAGE, maps.NSIDE, np.uint8) + old[np.array([10**9])] = np.array([1], dtype=np.uint8) + old.write(str(path)) + monkeypatch.setattr(merge.build_index, "campaign_exposures", + lambda *_: [EXP]) + monkeypatch.setattr(sys, "argv", [ + "merge_exposure_maps.py", "--products-dir", str(tmp_path), + "--tile-list", "-", "--index-db", "-", + *[arg for name in ("npoint", "ntaken", "nexp", "nflagged") + for arg in (f"--{name}", str(tmp_path / f"{name}.hsp"))]]) + with pytest.raises(SystemExit, match=f"Exposure {EXP}.*") as error: + merge.main() + assert "--image-mef" in str(error.value) + assert "--persist-manifest" in str(error.value) diff --git a/tests/workflow/params_pin.json b/tests/workflow/params_pin.json index 37bde1770..df726b2f7 100644 --- a/tests/workflow/params_pin.json +++ b/tests/workflow/params_pin.json @@ -9,7 +9,7 @@ "exp_persist": "302e2837542bc1102430c27c81c600b7cda32e8bddcb5fd60d33950987609fff", "exp_psf": "87b9fd799a54130eb3dd9c78bb5124838854c23f365df7a991a53f4c6ac4c3ff", "exp_split": "8fc97dd48e76cd6c7fdf9ee388ba1c1a3dd63428cbe1961641fbebe76b05786e", - "exposure_maps": "15bc1023a759ad61b67cc27f1740a74bb3be4d268969073820766fb64d6aace1", + "exposure_maps": "da8d3af52e651793a1a158590b802f37cf03bf754e35afce4aab1030082737bf", "final_cat_merge": "e7f46859c4503a2220713d7bb2507555515d0a9632d780b20f14c59e32210023", "prepare_all_tiles": "b8f872a22adf014e25a7fa5198f49b71a6fe9e56042ed82b682bc8763970a844", "star_cat_merge": "6277450958474af5270982fa35360f2f237a29f7533c526ee9265dfd5acc07a0", @@ -26,7 +26,7 @@ "tile_vignets": "1d434b4a6a4bcdeb8f1e5e9a5e4f228b67ade400245b41b033a2c883fb6a3267" }, "schema": 1, - "sha256": "389e2d5d5be672d4a6477a882e56eafb8a445e6280e95362d3c951bb126938f6", + "sha256": "a0286385d66743ac26e7d8e07081c70a5639fc2af5a21bb2a5d12907861bc1d9", "unit_pre": { "exp_get_images": "dac6685a207dae3ea81d296ce53e9ba75d4636cf1f81e6399f2dda452680c992", "exp_psf": "87fc8ea1153a709ab7ba56542cfacc041f5f8b7f6cf2f38a9bc9d8004d1cecb4", diff --git a/tests/workflow/test_dag.py b/tests/workflow/test_dag.py index 80ee3680f..0084aee0e 100644 --- a/tests/workflow/test_dag.py +++ b/tests/workflow/test_dag.py @@ -123,6 +123,21 @@ def test_products_use_products_dir_and_run_name(campaign, dag): assert Path(dag.namespace["INDEX_DB"]) == campaign.index_db +def test_exposure_maps_are_campaign_targets(campaign, dag): + """Fitted-PSF campaigns publish all four maps on the persistent root.""" + jobs = dag.jobs_for("exposure_maps") + if campaign.psf_model == "fake": + assert not jobs + return + assert len(jobs) == 1 + expected = {campaign.products_dir / f"{name}_{campaign.name}.hsp" + for name in ("npoint", "ntaken", "nexp", "nflagged")} + assert set(map(Path, jobs[0].output)) == expected + assert expected <= set(map(Path, dag.namespace["exposure_maps_targets"]())) + for name in ("npoint", "ntaken", "nexp", "nflagged"): + assert f"--{name} " in jobs[0].shellcmd + + def test_tile_store_is_unique_per_campaign(campaign, tmp_path, resolve_dag): """Concurrent campaigns over the same tiles get distinct node-local stores; within a campaign, every tile_shape member of a tile names the same one.""" diff --git a/universes/committed.yaml b/universes/committed.yaml index 265df3264..6d41807fb 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -10,6 +10,7 @@ analyses: pixel_mask_source: instrument_flags_only psf_star_mask_veto: instrument_flags_only sky_mask_application: catalogue_columns + pointing_square: main_array_hull nexp_map_valid_psf_ccds: valid_psf_ccds defect_map_from_flags: count_flagged_pixels detection: diff --git a/workflow/README.md b/workflow/README.md index a415e1b55..4384afba2 100644 --- a/workflow/README.md +++ b/workflow/README.md @@ -305,20 +305,28 @@ 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). +- **Four HealSparse maps come out of the exposures.** At nside 131072, + `npoint_.hsp` counts filled pointing hulls of present main-array CCDs 0–35. + The hull fills CCD gaps and excludes ear CCDs 36–39, regardless of PSF success. + `ntaken_.hsp` counts exposures covering each pixel with any CCD's DATASEC, ears included. + `nexp_.hsp` counts only CCDs with a validated PSF; `ntaken - nexp` measures CCD-level PSF losses. + These three counts are uint8. + `nflagged_.hsp` (uint16) counts flagged CCD pixel centres inside validated CCD DATASEC coverage, with ~74 CCD pixels per sky pixel. + The overscan border contributes neither coverage nor flags. + `exp_maps` writes one packed uint16 fragment beside each PSF tar: the low three bits store 1 = pointing, 2 = taken, 4 = validated, and bits 3–10 store flag counts clipped at 255. + Scalar fragments use compressed image storage (~2.6 MB per exposure); record-array table storage costs ~31 MB. + `exposure_maps` sums the fragments; nothing in the workflow reads the four maps. + `npoint` follows the v1.X npoint3 coverage-cut lineage, `nexp` supports an `N_EPOCH` cut, and `nflagged` diagnoses detector defects rather than exposures lost to the defect veto. + A plain uint8 fragment is rejected with an exposure-specific backfill command. + For reclaimed scratch stores, raw-MEF mode rewrites both the fragment and its completion manifest (the fragment is the durable DAG edge): + + ```bash + python workflow/scripts/exp_maps.py --exp "$exp" \ + --image-mef "$raw/${exp}p.fits.fz" --flag-mef "$raw/${exp}p.flag.fits.fz" \ + --persist-manifest "$products/exp/${exp:0:2}/$exp/manifests/exp_persist.json" \ + --fragment "$products/exp/${exp:0:2}/$exp/maps/maps-$exp.hsp" \ + --manifest "$products/exp/${exp:0:2}/$exp/manifests/exp_maps.json" + ``` - **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 20780e1be..04a24e0a0 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -879,9 +879,9 @@ def star_cat_exposures(): # --- 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 +# exp_maps rasterizes pointing hulls, all CCD DATASEC footprints, and validated +# CCD coverage and flags into one persistent fragment; exposure_maps sums these +# into npoint, ntaken, nexp and nflagged campaign maps. 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. @@ -891,6 +891,14 @@ def prod_exp_maps(exp): return f"{prod_exp_dir(exp)}/maps/maps-{exp}.hsp" +def npoint_map(): + return f"{PRODUCTS_DIR}/npoint_{CAMPAIGN}.hsp" + + +def ntaken_map(): + return f"{PRODUCTS_DIR}/ntaken_{CAMPAIGN}.hsp" + + def nexp_map(): return f"{PRODUCTS_DIR}/nexp_{CAMPAIGN}.hsp" @@ -918,7 +926,8 @@ def exposure_maps_targets(): 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 []) + return live + ([npoint_map(), ntaken_map(), nexp_map(), nflagged_map()] + if inputs else []) # --- sizing the two merges (D4) --------------------------------------------- diff --git a/workflow/rules/exposure.smk b/workflow/rules/exposure.smk index 12ad09a6f..6a401691c 100644 --- a/workflow/rules/exposure.smk +++ b/workflow/rules/exposure.smk @@ -19,9 +19,10 @@ 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_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_maps`` exports the pointing hull, all CCD DATASEC footprints, +validated-PSF coverage and instrument flag counts into a HealSparse fragment. +``exposure_maps`` sums these into four campaign maps; nothing in this +workflow reads them 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 @@ -184,13 +185,14 @@ rule exp_persist: # --- 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 +# covered by its pointing hull, all CCDs, and validated CCDs, with flag counts # (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. +# the valid-PSF CCDs. # @sc [decision:masking.defect_map_from_flags] # @sc [decision:masking.nexp_map_valid_psf_ccds] +# @sc [decision:masking.pointing_square] rule exp_maps: input: persist = lambda wc: prod_exp_manifest(wc.exp, "exp_persist") @@ -365,18 +367,20 @@ rule star_cat_merge: # --- 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 +# /{npoint,ntaken,nexp,nflagged}_.hsp: pointing hulls, +# all CCD DATASEC exposures, validated exposures, and their flag counts. 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). +# the job's own rediscovery of the set follow star_cat_merge. The four maps +# require 5 MiB per nside-128 coverage pixel: ~13 per exposure, capped at +# the ~23k of the UNIONS footprint. # @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: + npoint = npoint_map(), + ntaken = ntaken_map(), nexp = nexp_map(), nflagged = nflagged_map() params: @@ -388,7 +392,7 @@ rule exposure_maps: threads: 1 resources: mem_mb = lambda wc, attempt: capped_mem(attempt * ( - 1000 + 3.2 * min(13 * len(exposure_maps_exposures()), 23_000)), + 1000 + 5.3 * min(13 * len(exposure_maps_exposures()), 23_000)), "exposure_maps"), runtime = lambda wc, attempt: attempt * ( 30 + len(exposure_maps_exposures()) // 30) @@ -397,4 +401,5 @@ rule exposure_maps: f"python {SCRIPTS}/merge_exposure_maps.py" " --products-dir '{params.products_dir}'" " --tile-list '{params.tile_list}' --index-db '{params.index_db}'" + " --npoint {output.npoint} --ntaken {output.ntaken}" " --nexp {output.nexp} --nflagged {output.nflagged}" diff --git a/workflow/scripts/exp_maps.py b/workflow/scripts/exp_maps.py index d814856df..7cb4e9601 100644 --- a/workflow/scripts/exp_maps.py +++ b/workflow/scripts/exp_maps.py @@ -1,71 +1,53 @@ -"""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. +"""Rasterize one exposure into a packed uint16 HealSparse fragment. + +Bits 0–2 store coverage: 1 = pointing hull, 2 = CCD DATASEC, +4 = validated-PSF CCD DATASEC. Bits 3–10 store nonzero flagged CCD pixel +centre counts inside validated CCDs' coverage, clipped at 255. +A scalar map uses HealSparse's compressed image storage; record-array maps +use uncompressed tables and cost about 13 times as much disk space. +The maps use nside 131072 over coverage 128, aligned with the sky masks. + +Headers and flag arrays come from either the split scratch store or delivered +image/flag MEFs (CCD + 1 is the HDU index). The exp_persist manifest names the +validated CCDs. Both the fragment and exp_maps manifest are written, including +in raw-MEF backfills; durable_edge reads the fragment for reclaimed stores. """ import argparse import filecmp import json import warnings +from contextlib import ExitStack from fnmatch import fnmatch from pathlib import Path import healsparse as hsp import hpgeom as hpg import numpy as np +import persist_exp from astropy.io import fits from astropy.wcs import WCS - -import persist_exp +from scipy.spatial import ConvexHull NSIDE = 131072 NSIDE_COVERAGE = 128 +FRAGMENT_DTYPE = np.dtype(np.uint16) +COVER_MASK = 7 +FLAG_SHIFT = 3 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.""" + """Return CCD numbers with a PSF model from exp_persist 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) +def ccd_wcs(header): + """Return WCS and DATASEC bounds (1-based inclusive) from a CCD header.""" if "DATASEC" in header: x, y = header["DATASEC"].strip("[]").split(",") bounds = [int(v) for v in x.split(":") + y.split(":")] @@ -76,17 +58,45 @@ def ccd_wcs(image: Path): return WCS(header), bounds +def outer_corners(wcs, bounds): + """RA/Dec of DATASEC outer corners, half a pixel beyond its centres.""" + x0, x1, y0, y1 = bounds + return 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) + + 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) + ra, dec = outer_corners(wcs, bounds) + return hpg.query_polygon(NSIDE, ra, dec, nest=True) + + +def pointing_pixels(corners): + """Rasterize the main-array hull in a tangent plane about its centre. + + @sc [decision:masking.pointing_square] + CCDs 0–35 supply DATASEC outer corners; ears 36–39 are excluded and + gaps are filled. Pixel centres inside the convex hull count once. + """ + if not corners: + return np.empty(0, dtype=np.int64) + ra, dec = np.concatenate(corners, axis=1) + r, d = np.deg2rad([ra, dec]) + xyz = np.array([np.cos(d) * np.cos(r), np.cos(d) * np.sin(r), np.sin(d)]) + centre = xyz.mean(axis=1) + centre /= np.linalg.norm(centre) + tangent = WCS(naxis=2) + tangent.wcs.ctype = ["RA---TAN", "DEC--TAN"] + tangent.wcs.crval = [np.rad2deg(np.arctan2(centre[1], centre[0])) % 360, + np.rad2deg(np.arcsin(centre[2]))] + xy = np.column_stack(tangent.all_world2pix(ra, dec, 0)) + hull = xy[ConvexHull(xy).vertices] + ra, dec = tangent.all_pix2world(hull[:, 0], hull[:, 1], 0) 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.""" + """Sky pixels holding flagged DATASEC pixel centres, and their counts.""" 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) @@ -94,9 +104,78 @@ def flagged_counts(wcs, flags: np.ndarray, bounds): return_counts=True) +def rasterize(ccds, validated): + """Build a fragment from an iterable of (CCD number, header, flag array). + + @sc [decision:masking.nexp_map_valid_psf_ccds] + All CCD DATASEC footprints count as taken; only validated CCDs count as + successful. A footprint counts an exposure once even at overlapping edges. + @sc [decision:masking.defect_map_from_flags] + Flags count only inside validated CCD DATASEC coverage, not overscan. + """ + fragment = hsp.HealSparseMap.make_empty( + NSIDE_COVERAGE, NSIDE, FRAGMENT_DTYPE) + records, corners = {}, [] + present = set() + for ccd, header, flags in ccds: + present.add(ccd) + wcs, bounds = ccd_wcs(header) + if ccd < 36: + corners.append(outer_corners(wcs, bounds)) + covered = coverage_pixels(wcs, bounds) + values = fragment.get_values_pix(covered) + values |= 2 | (4 if ccd in validated else 0) + fragment[covered] = values + record = {"covered": int(covered.size), "validated": ccd in validated} + if ccd in validated: + pixels, counts = flagged_counts(wcs, flags, bounds) + keep = np.isin(pixels, covered, assume_unique=True) + values = fragment.get_values_pix(pixels[keep]) + flagged = np.minimum((values >> FLAG_SHIFT) + counts[keep], 255) + values = ((values & COVER_MASK) + | (flagged.astype(np.uint16) << FLAG_SHIFT)) + fragment[pixels[keep]] = values + record["flagged_pixels"] = int(counts.sum()) + records[ccd] = record + if set(validated) - present: + raise ValueError( + f"Validated CCDs absent from source: {set(validated) - present}") + if not present: + raise ValueError("No CCD images found") + pixels = pointing_pixels(corners) + values = fragment.get_values_pix(pixels) + values |= 1 + fragment[pixels] = values + return fragment, records + + +def source_ccds(args): + """Yield CCD headers and flags from splits or paired delivered MEFs.""" + if args.exp_dir is not None: + split = args.exp_dir / SPLIT_DIR + images = sorted(split.glob(f"image-{args.exp}-*.fits")) + for image in images: + ccd = int(image.stem.rsplit("-", 1)[1]) + yield ccd, fits.getheader(image), fits.getdata( + split / f"flag-{args.exp}-{ccd}.fits") + else: + with ExitStack() as stack: + images = stack.enter_context(fits.open(args.image_mef)) + flags = stack.enter_context(fits.open(args.flag_mef)) + if len(images) != len(flags): + raise ValueError( + "Image and flag MEFs have different HDU counts") + for ccd in range(len(images) - 1): + # Match split_exp's flag conversion, including delivered + # floating-point arrays and int16 overflow semantics. + yield (ccd, images[ccd + 1].header, + flags[ccd + 1].data.astype(np.int16)) + # Release decompressed flag arrays between CCDs. + del flags[ccd + 1].data + + 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).""" + """Replace dest unless bytes match, preserving unchanged product mtimes.""" if dest.exists() and filecmp.cmp(tmp, dest, shallow=False): tmp.unlink() else: @@ -104,44 +183,35 @@ def write_stable(tmp: Path, dest: Path) -> None: def main() -> None: + """Write a fragment and completion manifest from the selected source.""" parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument("--exp-dir", required=True, type=Path) + source = parser.add_mutually_exclusive_group(required=True) + source.add_argument("--exp-dir", type=Path) + source.add_argument("--image-mef", type=Path) + parser.add_argument("--flag-mef", 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())} - + if (args.image_mef is not None) != (args.flag_mef is not None): + parser.error("--image-mef and --flag-mef must be supplied together") + fragment, ccds = rasterize(source_ccds(args), + valid_ccds(args.persist_manifest, args.exp)) 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} + "format": "packed-cover-nflagged-v1", "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)") + print(f"[exp_maps] {args.exp}: {len(ccds)} CCD(s), " + f"{sum(c['validated'] for c in ccds.values())} validated, " + f"{fragment.valid_pixels.size} sky pixel(s)") if __name__ == "__main__": diff --git a/workflow/scripts/merge_exposure_maps.py b/workflow/scripts/merge_exposure_maps.py index d5a79caee..d106d8fc5 100644 --- a/workflow/scripts/merge_exposure_maps.py +++ b/workflow/scripts/merge_exposure_maps.py @@ -1,51 +1,61 @@ -"""Sum a campaign's per-exposure fragments into its two HealSparse maps: the -shell of the ``exposure_maps`` rule. +"""Sum per-exposure fragments into four campaign HealSparse maps. + +* ``npoint_.hsp`` (uint8): filled main-array pointing hulls. +* ``ntaken_.hsp`` (uint8): CCD DATASEC coverage, including ears. * ``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. + with a valid PSF model covers it. * ``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. +All 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. +Missing fragments are omitted and reported; incompatible formats fail with a +backfill command. The maps are rebuilt from every fragment on each run. """ import argparse from pathlib import Path +import build_index import healsparse as hsp import numpy as np - -import build_index -from exp_maps import NSIDE, NSIDE_COVERAGE +from exp_maps import FLAG_SHIFT, FRAGMENT_DTYPE, 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 the fragment path matching Snakefile's 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) + """Add ``values`` to ``acc`` at ``pixels``; a sum past the map's dtype + is an error, since a wrapped count reads as zero, i.e. uncovered.""" + total = acc.get_values_pix(pixels).astype(np.int64) + values + if total.size and total.max() > np.iinfo(acc.dtype).max: + raise SystemExit(f"merge_exposure_maps: a pixel count exceeds " + f"{np.iinfo(acc.dtype).max} ({acc.dtype}); widen " + f"the map dtype") + acc[pixels] = total.astype(acc.dtype) def main() -> None: + """Sum compatible fragments into the four campaign maps.""" 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("--npoint", required=True, type=Path) + parser.add_argument("--ntaken", 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) + npoint = hsp.HealSparseMap.make_empty(NSIDE_COVERAGE, NSIDE, np.uint8) + ntaken = hsp.HealSparseMap.make_empty(NSIDE_COVERAGE, NSIDE, np.uint8) nexp = hsp.HealSparseMap.make_empty(NSIDE_COVERAGE, NSIDE, np.uint8) nflagged = hsp.HealSparseMap.make_empty(NSIDE_COVERAGE, NSIDE, np.uint16) missing = [] @@ -55,24 +65,41 @@ def main() -> None: missing.append(exp) continue fragment = hsp.HealSparseMap.read(str(path)) + if fragment.dtype != FRAGMENT_DTYPE: + persist = path.parent.parent / "manifests" / "exp_persist.json" + manifest = path.parent.parent / "manifests" / "exp_maps.json" + raise SystemExit( + f"Exposure {exp}: incompatible fragment {path}. " + "Backfill with:\n" + f"python workflow/scripts/exp_maps.py --exp {exp} " + f"--image-mef /{exp}p.fits.fz " + f"--flag-mef /{exp}p.flag.fits.fz " + f"--persist-manifest {persist} --fragment {path} " + f"--manifest {manifest}") 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) + for acc, bit in ((npoint, 1), (ntaken, 2), (nexp, 4)): + selected = pixels[(values & bit) != 0] + add(acc, selected, np.ones(selected.size, dtype=np.uint8)) + counts = values >> FLAG_SHIFT + flagged = counts > 0 + add(nflagged, pixels[flagged], counts[flagged]) if i % 500 == 0: - print(f"[exposure_maps] {i + 1}/{len(exposures)} {exp}", flush=True) + 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)): + raise SystemExit("merge_exposure_maps: no exposure of the campaign " + "has a fragment") + for acc, out in ((npoint, args.npoint), (ntaken, args.ntaken), + (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}, " + print(f"[exposure_maps] {used} exposure(s) -> {args.npoint.name}, " + f"{args.ntaken.name}, {args.nexp.name}, " f"{args.nflagged.name}") if missing: print(f"[exposure_maps] WARNING: {len(missing)} exposure(s) in the "