Exposure maps: npoint (pointing squares) and ntaken beside nexp

astra.yaml

@@ -235,28 +235,52 @@ analyses:
         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.
+          npoint_<run>.hsp, ntaken_<run>.hsp, nexp_<run>.hsp and
+          nflagged_<run>.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: >-
           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.
+          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

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)

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",

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."""

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:

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_<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).
+- **Four HealSparse maps come out of the exposures.** At nside 131072,
+  `npoint_<run>.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_<run>.hsp` counts exposures covering each pixel with any CCD's DATASEC, ears included.
+  `nexp_<run>.hsp` counts only CCDs with a validated PSF; `ntaken - nexp` measures CCD-level PSF losses.
+  These three counts are uint8.
+  `nflagged_<run>.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