diff --git a/deployment/aws/build_function.sh b/deployment/aws/build_function.sh index ee37b8228..5cab0e31c 100755 --- a/deployment/aws/build_function.sh +++ b/deployment/aws/build_function.sh @@ -126,10 +126,12 @@ ITEM_COUNT=$(ls -1 "$BUILD_DIR" | wc -l) echo "" echo "Function code: ${UNZIPPED_SIZE} (${UNZIPPED_BYTES} bytes)" -# Function code budget: 30MB leaves room for the ~220MB layer -FUNCTION_BUDGET=$((30 * 1024 * 1024)) +# Function code budget: 32MB (espg ruling 2026-08-24, PR #511 question 1) — +# an early-warning tripwire under AWS's 50MB direct-upload zip limit, leaving +# room for the ~220MB layer; mirrored in tests/test_lambda_build.py. +FUNCTION_BUDGET=$((32 * 1024 * 1024)) if [ "$UNZIPPED_BYTES" -gt "$FUNCTION_BUDGET" ]; then - echo "WARNING: Function code exceeds 30MB budget!" + echo "WARNING: Function code exceeds 32MB budget!" echo " Top directories by size:" du -sh "$BUILD_DIR"/*/ 2>/dev/null | sort -rh | head -10 fi diff --git a/docs/api/catalog.md b/docs/api/catalog.md index d5208acf5..a4f7ecce6 100644 --- a/docs/api/catalog.md +++ b/docs/api/catalog.md @@ -152,3 +152,82 @@ cannot see. ::: zagg.catalog.polygon_to_bbox ::: zagg.catalog.load_antarctic_basins + +## Closest-observation pairing (issue #509) + +One raster store (e.g. Sentinel-2 L2A) can serve several point-cloud +reference stores (ATL03 + GEDI): for every reference *epoch* a shard's +stores actually observed, ingest the single **nearest** acquisition from the +raster catalog. The pairing is a property of the ingest *query*, not the +store schema — the raster store stays a plain raster store, which granules +were ingested *is* the pairing, and coincidence at read time is toc +intersection. + +Epochs are **store-derived**, never catalog-derived: each reference store's +`coverage.toc` sibling (spec §10.5) records per-shard word-set covers of the +data that actually landed, quantized at temporal order 18 (2^45 ns ≈ 9.77 h +buckets). The builder expands each cover word into its constituent buckets +and takes one epoch per bucket midpoint — good to ±4.9 h against Sentinel-2's +~4.3-day revisit. Granule catalogs would inherit the CMR-hull +over-assignment (~70 assigned granules vs 49 contributing pass-days on a +measured Californian shard); covers reflect contribution, not assignment. + +```python +from zagg.catalog.closest_obs import closest_obs_shardmap +from zagg.catalog.sources import Catalog +from zagg.grids import HealpixGrid +import numpy as np + +grid = HealpixGrid(9, 13) # parent_order must equal the covers' shard order +s2 = Catalog.from_geoparquet("catalog_s2_ca.parquet") + +# Size the run first — the dry-run builds nothing and prices the fan-out: +est = closest_obs_shardmap( + s2, + ["s3://bucket/atl03_store", "s3://bucket/gedi_store"], + grid=grid, + aoi="california.geojson", + max_time_offset=np.timedelta64(3, "D"), + max_granules_per_shard=200, # the same gate the build below applies + estimate=True, +) +# violations is [] unless the gate is passed here too -- estimate returns +# before the build's raise, so this is the safe way to size it. +est["histogram"], est["max_cost_usd"], est["violations"] + +# Then build the map; dispatch consumes it like any other ShardMap: +sm = closest_obs_shardmap( + s2, + ["s3://bucket/atl03_store", "s3://bucket/gedi_store"], + grid=grid, + aoi="california.geojson", + max_time_offset=np.timedelta64(3, "D"), + max_granules_per_shard=200, +) +sm.to_json("s2_closest_obs.json") +``` + +Everything refuses or records **loudly**, never silently: a reference store +with no readable `coverage.toc` raises (sweep the store first); an epoch +whose nearest acquisition lies beyond `max_time_offset` selects nothing and +is recorded per-epoch in `metadata["closest_obs"]["dropped"]` with its +near-miss offset; a shard past `max_granules_per_shard` raises naming the +worst shards (`estimate=True` reports the violations instead, so the gate +can be sized first); a cover block coarsened below the §10.5 pin is warned +about and reported in `coarsened_orders` — and under a `max_time_offset`, +epochs whose coarse-bucket half-span exceeds the stated offset cannot be +paired to that precision, so they drop into the ledger as their own category +(`epochs_dropped_low_resolution`, rows naming the block's effective order; +espg tolerance ruling 2026-08-24). Selected granule entries carry +`paired_epochs` / `epoch_offsets_ns` provenance so the paired product is +reconstructable from the manifest alone. Epochs are bucket midpoints — +size `max_time_offset` with `ReferenceEpochs.tolerance()`'s half-bucket +slack in mind. + +::: zagg.catalog.closest_obs.reference_epochs + +::: zagg.catalog.closest_obs.ReferenceEpochs + +::: zagg.catalog.closest_obs.nearest_acquisitions + +::: zagg.catalog.closest_obs.closest_obs_shardmap diff --git a/src/zagg/catalog/closest_obs.py b/src/zagg/catalog/closest_obs.py new file mode 100644 index 000000000..32948bbc9 --- /dev/null +++ b/src/zagg/catalog/closest_obs.py @@ -0,0 +1,887 @@ +"""Closest-observation ingest builder: cover-driven epochs -> paired shard map. + +Issue #509 — the consumer of the #489/#507 ``coverage.toc`` surface. One +Sentinel-2 (or any raster) store serves several point-cloud reference stores +(ATL03 + GEDI): for every reference *epoch* a shard's covers claim, the +builder selects the single **nearest** acquisition from the raster catalog — +closest observation, not a two-sided bracket; multiple passes bracket +naturally (espg design ruling, 2026-08-23/24). + +Epochs are **store-derived**: each reference store's ``coverage.toc`` sibling +(spec §10.5) carries per-shard word-set covers quantized at order 18 +(2^45 ns ≈ 9.77 h buckets). A cover *word* is a maximal RUN of those buckets +(``toc_normalize`` coalesces ranges that merely abut), so the epochs are the +midpoints of a word's **constituent buckets**, one per bucket — each within +±4.9 h of every instant its bucket covers, against Sentinel-2's ~4.3-day +revisit. Granule catalogs are *not* an +epoch source — the leaf sub-maps record the dispatched assignment verbatim +and inherit the CMR-hull over-assignment (~70 assigned vs 49 contributing +pass-days on shard ``3231422244``-class cases); covers reflect only data +that landed. + +The pairing is a property of the ingest *query*, not the store schema: the +raster store stays a plain raster store, which granules were ingested **is** +the pairing, and coincidence at read time is toc intersection. + +Word semantics (bit layout, decode, midpoints) are mortie's; the §10.5 cover +accessors are :mod:`zagg.coverage_toc`'s. This module owns only the join: +covers -> epochs, epochs -> nearest acquisitions, and the resulting +:class:`~zagg.catalog.shardmap.ShardMap`. +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass, field + +import numpy as np + +from zagg.coverage_toc import TEMPORAL_COVER_ORDER + +logger = logging.getLogger(__name__) + + +def _word_midpoints(words: np.ndarray, order: int = TEMPORAL_COVER_ORDER) -> np.ndarray: + """Cover words -> UTC ``datetime64[ns]`` midpoints of their *buckets*. + + A cover word is not one bucket. :func:`zagg.coverage_toc.quantize_words` + widens each instant to an aligned order-``order`` bucket and then + canonicalizes with ``toc_normalize``, which "coalesces ranges that merely + abut" (§10.5) — so a word is a maximal RUN of contiguous buckets and its + envelope midpoint would name one epoch per *campaign*, not per pass. This + expands every word back into its constituent buckets and emits one + midpoint each, which is what restores the ±4.9 h bound: a bucket spans + ``2**(63 - order)`` ns, so its midpoint is within half a bucket of every + instant the bucket covers. + + ``toc2time`` decodes a word's conservative envelope ``(start, end)`` on + the internal-ns scale — ``end`` exclusive for a range, ``end == start`` + for an exact timestamp — so the last covered instant is + ``max(end, start + 1) - 1``, the same uniform rule ``quantize_words`` + applies, and the covered buckets are ``start >> k`` through ``last >> k`` + inclusive with ``k = 63 - order``. Bucket midpoints use that same rule + within the bucket: ``(b << k) + 2**(k - 1) - 1``. + + One pass that straddles a bucket edge (its word's envelope reaches into + the neighbouring bucket) yields **two** epochs rather than one. That is + benign over-selection, not error: both epochs sit within half a bucket of + the pass, both pick the same nearest acquisition, and the builder dedupes + granule ids per shard. + + ``order`` is the block's *effective* temporal order (§10.5 lets a block + coarsen below the object's pin), defaulting to + :data:`zagg.coverage_toc.TEMPORAL_COVER_ORDER`. + """ + import mortie + + words = np.asarray(words, dtype=np.uint64) + if words.size == 0: + return np.empty(0, dtype="datetime64[ns]") + k = int(63 - int(order)) + start, end = mortie.toc2time(words) + start = np.atleast_1d(np.asarray(start, dtype=np.uint64)) + end = np.atleast_1d(np.asarray(end, dtype=np.uint64)) + last = np.maximum(end, start + np.uint64(1)) - np.uint64(1) + # Bucket index run [b0, b1] per word; expanded with repeat/arange + # arithmetic rather than a Python loop over words. + b0 = (start >> np.uint64(k)).astype(np.int64) + b1 = (last >> np.uint64(k)).astype(np.int64) + counts = b1 - b0 + 1 + offsets = np.cumsum(counts) - counts + within = np.arange(int(counts.sum()), dtype=np.int64) - np.repeat(offsets, counts) + buckets = np.unique(np.repeat(b0, counts) + within).astype(np.uint64) + half = np.uint64((1 << k) // 2 - 1) + mid = np.minimum((buckets << np.uint64(k)) + half, np.uint64(mortie.TOC_MAX_NS)) + return np.asarray(mortie.to_datetime64(mid), dtype="datetime64[ns]") + + +def _aoi_shard_set(aoi, order: int) -> set[int] | None: + """Resolve an ``aoi`` argument to the set of shard keys it covers. + + ``None`` passes through (no restriction). Accepted forms, matching the + catalog layer's existing vocabulary: + + - ``mortie.Moc`` — cast to the flat cell list at ``order``; + - ``str`` — a GeoJSON path (:func:`zagg.catalog.load_polygon`); + - ``[(lats, lons), ...]`` ring parts — the ``coverage``/``region`` form. + """ + if aoi is None: + return None + import mortie + + if isinstance(aoi, mortie.Moc): + return {int(c) for c in aoi.to_order(order)} + if isinstance(aoi, str): + from zagg.catalog import load_polygon + + aoi = load_polygon(aoi) + from zagg.grids.aoi import healpix_aoi_moc + + moc = healpix_aoi_moc(aoi, order) + return {int(c) for c in mortie.moc_to_order(moc, order)} + + +def _shard_order(cover: dict, root: str) -> int: + """A cover object's declared shard (dispatch) ``order``, or a loud refusal. + + §10.5 makes ``order`` mandatory; a body missing it (or carrying a + non-integer) is debris this builder cannot key shards from, and the bare + ``int(...)`` it used to get raised an opaque ``TypeError`` naming nothing. + """ + value = cover.get("order") + try: + return int(cover["order"]) + except (KeyError, TypeError, ValueError) as e: + raise ValueError( + f"reference_epochs: store {root!r} declares a non-integer cover shard order " + f"{value!r} — §10.5 requires it, and shard keys are parsed against it" + ) from e + + +@dataclass +class ReferenceEpochs: + """Per-shard reference epochs decoded from one or more store covers. + + Attributes + ---------- + order : int + The covers' shard (dispatch) order — every contributing store must + agree on it, and the raster grid the epochs pair against must match. + epochs : dict of int -> np.ndarray + Shard key (packed morton word, the canonical in-memory form D1 + renders as a decimal string externally) -> sorted unique UTC + ``datetime64[ns]`` epoch midpoints, the union across the contributing + stores' covers. Shards whose cover block decodes to an empty word set + are omitted (they claim nothing). + stores : list of str + The store roots that contributed, in the order given (provenance). + orders : dict of int -> int + Shard key -> the **coarsest** temporal order that contributed to that + shard's epochs. Normally the pinned + :data:`~zagg.coverage_toc.TEMPORAL_COVER_ORDER`, but §10.5 lets a + block coarsen below the pin to fit the cover cap, and a coarser order + means a wider midpoint bound (``2**(62 - order)`` ns). Phase 2's + ``max_time_offset`` reasoning reads this rather than assuming the pin. + """ + + order: int + epochs: dict[int, np.ndarray] + stores: list[str] = field(default_factory=list) + orders: dict[int, int] = field(default_factory=dict) + epoch_orders: dict[int, np.ndarray] = field(default_factory=dict) + """Shard key -> int64 array row-aligned with ``epochs``: each epoch's own + effective temporal order (espg tolerance ruling, 2026-08-24). :attr:`orders` + is the per-shard COARSEST — right for a headline warning, too blunt for a + precision gate: one shard can mix a pinned store's epochs with a coarsened + store's, and only the coarse ones fail a stated ``max_time_offset``. A + midpoint claimed at two orders keeps the finest.""" + + @property + def total(self) -> int: + """Total epoch count across shards (post-union, post-dedupe).""" + return sum(e.size for e in self.epochs.values()) + + def tolerance(self, shard: int) -> np.timedelta64: + """Half a bucket at ``shard``'s effective order — its epoch bound.""" + order = self.orders.get(shard, TEMPORAL_COVER_ORDER) + return np.timedelta64(2 ** (62 - int(order)), "ns") + + +def reference_epochs(reference_stores, *, aoi=None, **store_kwargs) -> ReferenceEpochs: + """Per-shard epochs from the reference stores' ``coverage.toc`` covers. + + For each store root: fetch the §10.5 sibling + (:func:`zagg.coverage_toc.read_cover`), strict-decode its per-shard word + sets (:func:`zagg.coverage_toc.cover_words`), and expand each word into + its constituent buckets' midpoints (:func:`_word_midpoints`). Per shard + the result is the **union** across stores, deduplicated (espg ruling: one + raster store serves both sensors, epochs are the union across the + reference stores). + + The union is canonical at the **bucket** level, not the word level. Words + are post-``toc_normalize`` runs whose extent depends on that store's + *other* data, so two stores sharing one pass routinely emit + overlapping-but-**unequal** range words — a raw ``np.unique`` over words + would keep both and represent the shared pass twice, at two displaced + midpoints. Expanding to buckets first removes that degree of freedom: + the bucket grid is fixed by the order alone, so a shared pass contributes + the same bucket midpoint from every store and dedupes exactly. + + §10.5 lets a block coarsen **below** the object's pinned temporal order to + fit the cover cap, and its landed order is recorded in the block itself. + :func:`zagg.coverage_toc.cover_words` returns only the words, so this + reads each block's ``temporal_order`` straight from the grammar, expands + that block's words at **its** order, logs a warning whenever a + block sits below the pin (the read half of §10.5's "widening only, + *loudly recorded*"), and reports the coarsest contributing order per + shard on :attr:`ReferenceEpochs.orders`. + + A store that carries **no readable cover refuses loudly** — this builder + is cover-driven by design (store-derived epochs, never the granule + catalogs), so "no cover yet" means "sweep the store first", not "fall + back silently". Likewise a shard-order mismatch between stores: D1 ids + at two orders are not comparable. + + Parameters + ---------- + reference_stores : str or sequence of str + Store roots whose covers drive the epochs (e.g. ATL03 + GEDI). + aoi : optional + Restrict the shard set: a ``mortie.Moc``, a GeoJSON path, or + ``[(lats, lons), ...]`` ring parts (see :func:`_aoi_shard_set`). + **store_kwargs + Forwarded to the object-store open (region, credentials, ...). + + Returns + ------- + ReferenceEpochs + """ + from zagg.coverage_toc import COVER_NAME, cover_words, load_cover, read_cover + from zagg.grids.morton import morton_word + + if isinstance(reference_stores, str): + reference_stores = [reference_stores] + reference_stores = list(reference_stores) + if not reference_stores: + raise ValueError("reference_epochs: at least one reference store root is required") + + order: int | None = None + mids_by_shard: dict[int, list[tuple[np.ndarray, int]]] = {} + orders: dict[int, int] = {} + for root in reference_stores: + obj = read_cover(root, **store_kwargs) + cover = load_cover(obj) + if cover is None: + detail = ( + "no coverage.toc object" + if obj is None + else f"unreadable {COVER_NAME} (spec {obj.get('spec') if isinstance(obj, dict) else obj!r})" + ) + raise ValueError( + f"reference_epochs: store {root!r} has {detail} — epochs are cover-driven " + f"(spec §10.5, issue #509); run the rollup sweep that materializes the " + f"cover before pairing against this store" + ) + store_order = _shard_order(cover, root) + if order is None: + order = store_order + elif store_order != order: + raise ValueError( + f"reference_epochs: store {root!r} covers shard order {store_order}, " + f"previous stores cover order {order} — D1 ids at two orders are not " + f"comparable (spec §10.5)" + ) + # cover_words strict-decodes (and validates the object's pin); the + # per-block order it drops is read back off the same grammar. + decoded = cover_words(obj) or {} + pinned = int(cover.get("temporal_order", TEMPORAL_COVER_ORDER)) + blocks = cover.get("shards") or {} + for decimal, words in decoded.items(): + if not len(words): + continue + block = blocks.get(decimal) or {} + raw = block.get("temporal_order", pinned) + # Boundary REFUSAL, not a clip. §10.5's order check in + # ``coverage_toc._decode_cover_block`` is one-sided (ceiling only), + # so a corrupt block declaring a negative order decodes fine — and + # then ``62 - order`` overflows the int64 shift to a huge NEGATIVE + # half-span that silently PASSES the caller's precision bar, on + # midpoints ``_word_midpoints`` computed from a shift it cannot + # express either. Refuse at this read boundary; the durable + # one-line lower bound belongs beside that ceiling check. + if isinstance(raw, bool) or not isinstance(raw, int) or raw < 0: + raise ValueError( + f"reference_epochs: store {root!r} shard {decimal} cover block declares " + f"temporal_order {raw!r} — a block only ever coarsens BELOW the object's " + f"pin and never below 0 (spec §10.5); this block is corrupt, and decoding " + f"it would yield epochs no precision bar can hold" + ) + effective = int(raw) + if effective < TEMPORAL_COVER_ORDER: + logger.warning( + f"reference_epochs: store {root!r} shard {decimal} cover sits at " + f"temporal order {effective}, below the pinned {TEMPORAL_COVER_ORDER} " + f"(§10.5 cap coarsening) — its buckets span 2^{63 - effective} ns, so " + f"these epochs are good to ±2^{62 - effective} ns, not ±4.9 h" + ) + # Cover blocks are keyed by the D1 decimal id (the external + # string form, sign included); shard maps key on the packed + # morton word — parse at the boundary (issue #199). + shard = morton_word(decimal) + mids_by_shard.setdefault(shard, []).append( + (_word_midpoints(words, effective), effective) + ) + orders[shard] = min(orders.get(shard, effective), effective) + + assert order is not None # non-empty store list, every cover carried an order + keep = _aoi_shard_set(aoi, order) + epochs: dict[int, np.ndarray] = {} + epoch_orders: dict[int, np.ndarray] = {} + for shard in sorted(mids_by_shard): + if keep is not None and shard not in keep: + continue + parts = mids_by_shard[shard] + mids = np.concatenate([m for m, _ in parts]) + ords = np.concatenate([np.full(m.size, o, dtype=np.int64) for m, o in parts]) + # Dedupe exact midpoints keeping each epoch's FINEST claiming order + # (the precision gate reads per-epoch orders): sort by (mid, -order) + # so the first row of every equal-midpoint run is the finest. + ix = np.lexsort((-ords, mids.astype("int64"))) + mids, ords = mids[ix], ords[ix] + first = np.ones(mids.size, dtype=bool) + first[1:] = mids[1:] != mids[:-1] + mids, ords = mids[first], ords[first] + if mids.size: + epochs[shard] = mids + epoch_orders[shard] = ords + return ReferenceEpochs( + order, + epochs, + reference_stores, + {k: v for k, v in orders.items() if k in epochs}, + epoch_orders, + ) + + +def _cap_ns(max_time_offset) -> int | None: + """``max_time_offset`` as validated non-negative nanoseconds, or ``None``. + + One conversion for the selection gate (:func:`nearest_acquisitions`) and + the builder's cover-resolution gate, so both refuse the same inputs with + the same message: NaT, a duration that does not round-trip through + ``timedelta64[ns]`` (~292-year span), and a negative cap. + """ + if max_time_offset is None: + return None + offset = np.timedelta64(max_time_offset) + if np.isnat(offset): + raise ValueError(f"max_time_offset must be a real duration (got {max_time_offset!r})") + as_ns = offset.astype("timedelta64[ns]") + if as_ns.astype(offset.dtype) != offset: + raise ValueError( + "max_time_offset does not convert exactly to nanoseconds " + f"(got {max_time_offset!r}; timedelta64[ns] spans ~292 years)" + ) + cap = int(as_ns.astype("int64")) + if cap < 0: + raise ValueError(f"max_time_offset must be non-negative (got {max_time_offset!r})") + return cap + + +def nearest_acquisitions(epochs, times, *, max_time_offset=None): + """Nearest acquisition per epoch — the vectorized closest-1 core. + + The phase-2 selection of the closest-observation join (issue #509): for + each reference epoch, the single nearest acquisition — closest-1, never a + two-sided bracket (espg ruling); several epochs bracketing one + acquisition each select it, and the builder dedupes granules downstream. + + Parameters + ---------- + epochs : array-like of datetime64[ns] + One shard's reference epochs (:func:`reference_epochs`). Any order. + ``NaT`` refuses — see the note below. + times : array-like of datetime64[ns] + The shard's acquisition times, in catalog record order. Any order — + the returned selection indexes THIS array's positions. ``NaT`` + refuses — see the note below. + max_time_offset : np.timedelta64 or int, optional + An epoch whose nearest acquisition lies further than this selects + nothing. Exactly-at selects; one ns past does not. A bare int means + nanoseconds. ``None`` (default) always selects the nearest, however + far. Negative refuses, and so does a duration that will not convert + to nanoseconds (``np.timedelta64(1000, "Y")`` overflows int64 ns). Callers gating against the epochs' own + precision should widen by :meth:`ReferenceEpochs.tolerance` — a + cover epoch is a bucket midpoint, good to half a bucket, not exact. + + Returns + ------- + selection : np.ndarray of int64 + Per epoch, the index into ``times`` of the selected acquisition, or + ``-1`` where the epoch selects nothing (no acquisitions at all, or + nearest beyond ``max_time_offset``). + offsets : np.ndarray of timedelta64[ns] + Per epoch, the SIGNED offset ``times[nearest] - epoch`` of the + nearest acquisition — positive when the acquisition follows the + epoch. Reported for every epoch, dropped ones included (the loud + record a drop rides — the builder's report needs the near-miss + distance, not just the fact of the drop). ``NaT`` when ``times`` is + empty, and for the offset that will not fit ``timedelta64[ns]`` — + selection is exact over the whole ``datetime64[ns]`` span (~584 + years), but a *difference* past ~292 years is unrepresentable, so it + saturates to ``NaT`` rather than reporting a wrapped duration. Such + an epoch selects nothing under any cap (no finite tolerance reaches + it) and still selects its nearest with ``max_time_offset=None``. + + Raises + ------ + ValueError + If either input carries ``NaT``. A missing instant has no nearest + anything: ``NaT`` sorts last but casts to ``iinfo(int64).min``, so it + would leave ``ts`` unsorted (silently mis-pairing its neighbours) and + wrap the offset subtraction. :func:`reference_epochs` never emits + ``NaT`` and the catalog's time parser refuses a missing acquisition + time, so ``NaT`` here is caller debris — refused loudly, this + module's posture. + + Notes + ----- + A tie — an epoch exactly equidistant between two acquisitions — selects + the EARLIER acquisition, deterministically. Equal acquisition times (one + Sentinel-2 datatake stamps many granules with the same instant) are + broken by catalog record order: the FIRST record of the equal-time run, + whichever flank the epoch approaches it from. + """ + epochs = np.asarray(epochs, dtype="datetime64[ns]") + times = np.asarray(times, dtype="datetime64[ns]") + for name, arr in (("epochs", epochs), ("times", times)): + if np.isnat(arr).any(): + raise ValueError(f"{name} carries NaT ({int(np.isnat(arr).sum())} of {arr.size})") + cap = _cap_ns(max_time_offset) + selection = np.full(epochs.shape, -1, dtype=np.int64) + offsets = np.full(epochs.shape, np.timedelta64("NaT"), dtype="timedelta64[ns]") + if epochs.size == 0 or times.size == 0: + return selection, offsets + order = np.argsort(times, kind="stable") + ts = times[order].astype("int64") + e = epochs.astype("int64") + pos = np.searchsorted(ts, e) # left insertion point + # Flank distances as UNSIGNED magnitudes. Both are non-negative by + # construction (``e >= ts[pos - 1]`` and ``ts[pos] >= e``), and mod-2^64 + # subtraction of the int64 bit patterns is exact across the whole + # datetime64[ns] span — an int64 ``e - ts`` wraps silently for a pair + # more than ~292 years apart, which then passes the cap gate. + tsu, eu = ts.view(np.uint64), e.view(np.uint64) + far = np.uint64(2**64 - 1) # a missing flank, wider than any real span + left = np.where(pos > 0, eu - tsu[np.maximum(pos - 1, 0)], far) + right = np.where(pos < ts.size, tsu[np.minimum(pos, ts.size - 1)] - eu, far) + # Strict ``<`` keeps a tie on the LEFT (earlier) neighbor; an epoch equal + # to an acquisition has ``right == 0`` and selects it exactly. + take_right = right < left + nearest = np.where(take_right, np.minimum(pos, ts.size - 1), np.maximum(pos - 1, 0)) + # Equal acquisition times form one run in ``ts``; the left flank lands on + # its END and the right flank on its START, so snap to the run start — + # with a stable ``argsort`` that is the run's first catalog record, from + # either side. + nearest = np.searchsorted(ts, ts[nearest], side="left") + selection = order[nearest].astype(np.int64) + magnitude = np.where(take_right, right, left) + # A magnitude past int64 ns is a real distance the report cannot carry: + # saturate it to NaT rather than wrap it into a plausible-looking day. + fits = magnitude <= np.uint64(np.iinfo(np.int64).max) + signed = np.minimum(magnitude, np.uint64(np.iinfo(np.int64).max)).astype(np.int64) + offsets = np.where(take_right, signed, -signed).astype("timedelta64[ns]") + offsets = np.where(fits, offsets, np.timedelta64("NaT")) + if cap is not None: + selection = np.where(magnitude <= np.uint64(cap), selection, np.int64(-1)) + return selection, offsets + + +def _acquisition_times(entries: list[dict], shard: str) -> np.ndarray: + """Granule entries -> UTC ``datetime64[ns]`` acquisition instants. + + Reads the entry's ``datetime`` (the raster-source acquisition instant, + #218) falling back to ``time_start`` (the STAC acquisition-range start, + #246). A granule carrying neither refuses loudly — a catalog without + acquisition times cannot be temporally paired, and skipping the granule + would silently thin the product. + """ + from datetime import datetime, timezone + + out = np.empty(len(entries), dtype="datetime64[ns]") + for i, entry in enumerate(entries): + iso = entry.get("datetime") or entry.get("time_start") + if iso is None: + raise ValueError( + f"closest_obs_shardmap: granule {entry.get('id')!r} in shard {shard} " + f"carries no acquisition time (neither 'datetime' nor 'time_start') — " + f"the catalog cannot be temporally paired (issue #509)" + ) + dt = datetime.fromisoformat(str(iso).replace("Z", "+00:00")) + if dt.tzinfo is not None: + dt = dt.astimezone(timezone.utc).replace(tzinfo=None) + out[i] = np.datetime64(dt, "ns") + return out + + +def closest_obs_shardmap( + s2_catalog, + reference_stores, + *, + grid, + aoi=None, + max_time_offset=None, + max_granules_per_shard=None, + estimate=False, + backend="auto", + bytes_per_granule=None, + **store_kwargs, +): + """Closest-observation ingest map: covers -> epochs -> nearest granules. + + The issue-#509 builder. Epochs come from the reference stores' + ``coverage.toc`` covers (:func:`reference_epochs`); the raster catalog is + spatially assigned to shards through the existing machinery + (:meth:`~zagg.catalog.shardmap.ShardMap.build` — the stored-index / + batch-cover fast paths included); per shard, each epoch then selects its + single nearest acquisition (:func:`nearest_acquisitions`) and the + selected granules, deduplicated, become the shard's ingest list. The + result is a standard :class:`~zagg.catalog.shardmap.ShardMap` — JSON + round-trip, ``total_pairs`` bookkeeping — so dispatch consumes it + unchanged; the pairing is a property of this query, never of the raster + store's schema (espg ruling). + + Parameters + ---------- + s2_catalog : Catalog or str + The raster acquisition catalog (stac-geoparquet), or a path to one + (``Catalog.from_geoparquet``). Every record must carry an acquisition + time (``datetime``, or STAC ``start_datetime``). + reference_stores : str or sequence of str + Store roots whose covers drive the epochs (e.g. ATL03 + GEDI). + grid : HealpixGrid + The raster store's output grid. Its ``parent_order`` must equal the + covers' shard order — the map's shard keys and the covers' D1 ids + live on the same grid, and the emitted ``grid_signature`` is what the + ingest run validates against. + aoi : optional + Restrict the shard set (``mortie.Moc``, GeoJSON path, or ring + parts); intersected with the store coverage (:func:`reference_epochs`). + Ring parts / a GeoJSON path additionally ride as the spatial build's + ``region=``, scoping the raster intersection; a ``mortie.Moc`` has no + parts form, so it restricts the epoch side only and the spatial build + still covers the whole catalog bbox. + max_time_offset : np.timedelta64 or int (ns), optional + An epoch whose nearest acquisition lies beyond this selects nothing — + recorded per epoch in ``metadata["closest_obs"]["dropped"]`` and + warned about, never silent. It is also a **precision bar** on the + epochs themselves (espg tolerance ruling, 2026-08-24): an epoch from + a cover block coarsened far enough that its bucket half-span exceeds + this offset cannot be paired to the stated precision, and is dropped + into the ledger as its own category (rows carrying + ``temporal_order``/``cover_half_span_ns``; counted in + ``epochs_dropped_low_resolution``). Half-span exactly at the offset + stays pairable — the same side the selection gate's exactly-at rule + pins. ``None`` always selects the nearest, at whatever resolution the + covers offer (a single warning names the effective resolution when + any block sits below the §10.5 pin). + max_granules_per_shard : int, optional + Cost gate: a shard exceeding this REFUSES loudly (``ValueError`` + naming the worst shards) — never truncates. ``estimate=True`` reports + violations instead of raising, so the gate can be sized first. + estimate : bool + Dry-run: return the per-shard histogram + cost estimate dict (see + Notes) WITHOUT building the map. The espg-operated ingest runs are + cost-gated through this. + backend : {"auto", "spherely", "mortie"} + Forwarded to :meth:`ShardMap.build`. + bytes_per_granule : int, optional + Per-granule ingest volume for the ``estimate`` byte figure; without + it ``est_bytes`` is ``None`` (the catalog does not carry sizes). + **store_kwargs + Forwarded to the reference stores' object-store opens. + + Returns + ------- + ShardMap or dict + The ingest map — or, with ``estimate=True``, a dict: + ``{"shards", "granules", "pairs", "epochs_total", "epochs_paired", + "epochs_dropped", "epochs_dropped_low_resolution", "per_shard" (decimal -> granule count), + "histogram" (granule count -> shard count), "est_bytes", + "max_cost_usd", "violations"}``. ``max_cost_usd`` is the + :func:`zagg.dispatch.max_cost_usd` ceiling at the production worker + size (one invoke per shard, 900 s timeout) — a bound, not a forecast. + + Notes + ----- + A selected entry whose record carries no ``datetime`` (STAC's null-datetime + + ``start_datetime`` form) is emitted with ``datetime`` backfilled from + ``time_start`` — the same instant the pairing used, and the key raster + dispatch requires. Selected granule entries also gain two provenance keys + so the eventual paired product is reconstructable from the manifest alone: ``paired_epochs`` + (ISO instants of every epoch that selected the granule) and + ``epoch_offsets_ns`` (row-aligned SIGNED ``acquisition - epoch`` ns). + Do NOT derive a paired map with :meth:`ShardMap.reproject`: its + ``_granule_entry`` passthrough does not know these two keys, so a + reprojected map (the ``noop`` same-order branch included) drops the + provenance while ``metadata["closest_obs"]`` rides through describing the + SOURCE map — after a coarsen its shard ids name a shard set that no longer + exists. Rebuild at the target grid instead of reprojecting. + ``metadata["closest_obs"]`` records the query: the reference stores, the + epoch totals, every dropped epoch with its near-miss offset (``None`` + where there was no acquisition to measure against, so + ``epochs_total == epochs_paired + epochs_dropped`` holds), shards + whose epochs found no acquisition at all, and any cover blocks coarsened + below the §10.5 pin. Epochs are bucket midpoints, good to half a bucket + (:meth:`ReferenceEpochs.tolerance`) — size ``max_time_offset`` with that + slack in mind. + + The strict-AOI per-shard mask (``output.aoi_mask``, issue #101) is the + SPATIAL map's payload and is not carried onto the emitted map, so the + ``aoi_mask`` metadata claim is dropped with it: an ``output.aoi_mask`` + ingest run must compute its mask at run time rather than read one off + this manifest. + """ + from zagg.catalog.shardmap import ShardMap + from zagg.grids.morton import morton_decimal + + if isinstance(s2_catalog, str): + from zagg.catalog.sources import Catalog + + s2_catalog = Catalog.from_geoparquet(s2_catalog) + + parent_order = getattr(grid, "parent_order", None) + if parent_order is None: + raise ValueError( + "closest_obs_shardmap: grid must be a HEALPix grid (parent_order) — the " + "covers' D1 shard ids and the map's shard keys live on the same grid" + ) + # Resolve the aoi once. Ring parts (and a GeoJSON path, which loads to + # parts) also scope the SPATIAL build via ``region=``, so an AOI run stops + # intersecting the whole catalog bbox and discarding the outside shards + # afterwards; a ``mortie.Moc`` has no ring-parts form and stays epoch-side. + import mortie + + if isinstance(aoi, str): + from zagg.catalog import load_polygon + + aoi = load_polygon(aoi) + region = None if aoi is None or isinstance(aoi, mortie.Moc) else aoi + ref = reference_epochs(reference_stores, aoi=aoi, **store_kwargs) + if int(parent_order) != int(ref.order): + raise ValueError( + f"closest_obs_shardmap: grid.parent_order {int(parent_order)} != the covers' " + f"shard order {ref.order} — the emitted map would key shards on a different " + f"grid than the epochs (spec §10.5)" + ) + + cap_ns = _cap_ns(max_time_offset) + spatial = ShardMap.build(s2_catalog, grid, region=region, backend=backend) + spatial_idx = {int(k): i for i, k in enumerate(spatial.shard_keys)} + aoi_keys = _aoi_shard_set(aoi, ref.order) + + shard_keys: list[int] = [] + granules: list[list[dict]] = [] + dropped: list[dict] = [] + no_acquisitions: list[str] = [] + epochs_paired = 0 + low_resolution = 0 + for shard, epoch_arr in sorted(ref.epochs.items()): + decimal = morton_decimal(shard) + # Cover-resolution gate (espg tolerance ruling, 2026-08-24, thread + # r3845481805): an epoch whose bucket HALF-SPAN exceeds the caller's + # stated ``max_time_offset`` cannot be paired to that precision — the + # true pass sits anywhere inside the bucket, so any pairing within the + # cap would be arbitrary. Dropped loudly, per epoch, as its own ledger + # category (``temporal_order``/``cover_half_span_ns`` rows), BEFORE the + # nearest selection. Strictly greater: half-span exactly at the cap + # stays pairable, the same side the selection gate's exactly-at rule + # pins. With no cap the caller declared no precision bar — the build + # warns once (below) and proceeds; widening is lawful (§10.5). + # + # It runs AHEAD of the spatial lookup on purpose: unresolvability is a + # property of the epoch and its own cover block, so which category an + # epoch lands in must not flip on whether the raster catalog happens to + # reach the shard. Gating after the lookup made an unreached shard's + # coarse epochs ledger as no-acquisition rows and read + # ``epochs_dropped_low_resolution == 0``. + eo = ref.epoch_orders.get(shard) + if cap_ns is not None and eo is not None: + half_span = np.int64(1) << (np.int64(62) - eo) + unresolvable = half_span > cap_ns + if unresolvable.any(): + dropped.extend( + { + "shard": decimal, + "epoch": np.datetime_as_string(t), + "temporal_order": int(o), + "cover_half_span_ns": int(h), + } + for t, o, h in zip( + epoch_arr[unresolvable], eo[unresolvable], half_span[unresolvable] + ) + ) + low_resolution += int(unresolvable.sum()) + epoch_arr = epoch_arr[~unresolvable] + i = spatial_idx.get(shard) + if i is None: + no_acquisitions.append(decimal) + # The shard is named even when the gate took every epoch: this row + # is about the SHARD, not its epochs. Ledger the SURVIVING epochs + # too — a shard the catalog never reaches is the largest drop class + # in practice, and leaving it out of ``dropped`` made the numbers an + # operator reconciles read "nothing dropped". ``nearest_offset_ns`` + # is None -- no acquisition to measure against, the meaning the key + # already carries. + dropped.extend( + {"shard": decimal, "epoch": np.datetime_as_string(t), "nearest_offset_ns": None} + for t in epoch_arr + ) + continue + if epoch_arr.size == 0: + continue + entries = spatial.granules[i] + times = _acquisition_times(entries, decimal) + sel, off = nearest_acquisitions(epoch_arr, times, max_time_offset=max_time_offset) + off_ns = off.astype("int64") + chosen: dict[int, dict] = {} + for j in range(epoch_arr.size): + iso = np.datetime_as_string(epoch_arr[j]) + if sel[j] < 0: + dropped.append( + { + "shard": decimal, + "epoch": iso, + "nearest_offset_ns": None if np.isnat(off[j]) else int(off_ns[j]), + } + ) + continue + src = entries[sel[j]] + entry = chosen.setdefault( + int(sel[j]), + { + **src, + # STAC allows ``datetime: null`` beside start/end_datetime, + # but raster dispatch keys off ``datetime`` + # (``runner._raster_windowed_units``). Emit the instant the + # pairing actually used so the map dispatches as built. + "datetime": src.get("datetime") or src.get("time_start"), + "paired_epochs": [], + "epoch_offsets_ns": [], + }, + ) + entry["paired_epochs"].append(iso) + entry["epoch_offsets_ns"].append(None if np.isnat(off[j]) else int(off_ns[j])) + epochs_paired += 1 + if chosen: + shard_keys.append(shard) + granules.append([chosen[j] for j in sorted(chosen)]) + + coarse = {morton_decimal(k): o for k, o in ref.orders.items() if o < TEMPORAL_COVER_ORDER} + if dropped: + logger.warning( + f"closest_obs_shardmap: {len(dropped)} epoch(s) selected nothing " + f"(max_time_offset={max_time_offset!r}, no acquisitions in the shard, or a " + f"cover block too coarse for the offset); " + f"e.g. {dropped[:3]} — every drop is " + f"recorded in metadata['closest_obs']['dropped']" + ) + if low_resolution: + # The arm that DISCARDS epochs must be at least as loud as the no-cap + # arm below, which discards nothing and still gets a purpose-built + # line. The summary above cannot carry this: those rows never reached + # the selection, so its distance/catalog-gap causes are both wrong for + # them, and ``reference_epochs``' coarsening warning says only that the + # epochs are coarse — never that they were consequently dropped. + worst = min(coarse.values()) if coarse else None + logger.warning( + f"closest_obs_shardmap: {low_resolution} epoch(s) dropped as UNRESOLVABLE at " + f"max_time_offset={max_time_offset!r} — a coarsened cover, NOT distance to an " + f"acquisition: their block's bucket half-span exceeds the offset (coarsest " + f"temporal order {worst} across {len(coarse)} shard(s), e.g. " + f"{sorted(coarse)[:5]}) — every row is in metadata['closest_obs']['dropped'] " + f"with temporal_order/cover_half_span_ns, counted in " + f"metadata['closest_obs']['epochs_dropped_low_resolution']" + ) + if no_acquisitions: + logger.warning( + f"closest_obs_shardmap: {len(no_acquisitions)} shard(s) carry reference epochs " + f"but NO spatially-assigned acquisitions at all (catalog gap?): " + f"{no_acquisitions[:5]}" + ) + + counts = {morton_decimal(k): len(g) for k, g in zip(shard_keys, granules)} + violations = ( + sorted( + ((d, n) for d, n in counts.items() if n > max_granules_per_shard), + key=lambda kv: -kv[1], + ) + if max_granules_per_shard is not None + else [] + ) + + if coarse and cap_ns is None: + worst = min(coarse.values()) + logger.warning( + f"closest_obs_shardmap: cover blocks in {len(coarse)} shard(s) sit below the " + f"§10.5 pin (coarsest temporal order {worst}); effective epoch resolution is " + f"±2^{62 - worst} ns (~{2 ** (62 - worst) / 3.6e12:.1f} h) — no max_time_offset " + f"was set, so pairing proceeds at that resolution" + ) + closest_meta = { + "reference_stores": list(ref.stores), + "shard_order": int(ref.order), + "max_time_offset_ns": cap_ns, + "epochs_total": int(ref.total), + "epochs_paired": int(epochs_paired), + "epochs_dropped": len(dropped), + # The cover-resolution category's own count (espg tolerance ruling): + # these rows are inside ``dropped`` — the ledger invariant stays + # ``epochs_total == epochs_paired + epochs_dropped`` — but an operator + # sizing max_time_offset needs them told apart from near-misses. + "epochs_dropped_low_resolution": int(low_resolution), + "dropped": dropped, + "shards_without_acquisitions": no_acquisitions, + # Coverage disagreement, not self-inflicted clipping: ``ref.epochs`` is + # AOI-filtered, so an AOI-excluded shard is not a shard the reference + # stores never observed. + "spatial_shards_without_epochs": sum( + 1 + for k in spatial.shard_keys + if k not in ref.epochs and (aoi_keys is None or k in aoi_keys) + ), + "coarsened_orders": coarse, + } + + if estimate: + from zagg.dispatch import LAMBDA_MEMORY_GB, max_cost_usd + + histogram: dict[int, int] = {} + for n in counts.values(): + histogram[n] = histogram.get(n, 0) + 1 + distinct = len({g["id"] for shard in granules for g in shard}) + pairs = sum(len(g) for g in granules) + return { + **closest_meta, + "shards": len(shard_keys), + "granules": distinct, + "pairs": pairs, + "per_shard": counts, + "histogram": dict(sorted(histogram.items())), + "est_bytes": None if bytes_per_granule is None else pairs * int(bytes_per_granule), + "max_cost_usd": round( + max_cost_usd(len(shard_keys), LAMBDA_MEMORY_GB, timeout_s=900.0), 2 + ), + "violations": violations, + } + + if violations: + worst = ", ".join(f"{d}={n}" for d, n in violations[:5]) + raise ValueError( + f"closest_obs_shardmap: {len(violations)} shard(s) exceed " + f"max_granules_per_shard={max_granules_per_shard} (worst: {worst}) — refusing " + f"loudly rather than truncating; raise the gate, tighten max_time_offset/aoi, " + f"or size the run with estimate=True first (issue #509)" + ) + + meta = { + **spatial.metadata, + "total_shards": len(shard_keys), + "total_pairs": sum(len(g) for g in granules), + "granules_assigned": len({g["id"] for shard in granules for g in shard}), + "closest_obs": closest_meta, + } + # The spatial build's per-shard strict-AOI mask is NOT carried onto this + # derived map (``aoi_mask=None`` below), so the metadata must stop + # advertising one -- the same guard ``reproject`` applies to a derived map + # (``shardmap.py:1738``). + meta.pop("aoi_mask", None) + return ShardMap(spatial.grid_signature, shard_keys, granules, meta, None) + + +__all__ = [ + "ReferenceEpochs", + "closest_obs_shardmap", + "nearest_acquisitions", + "reference_epochs", +] diff --git a/tests/test_closest_obs.py b/tests/test_closest_obs.py new file mode 100644 index 000000000..1fba41b1a --- /dev/null +++ b/tests/test_closest_obs.py @@ -0,0 +1,1204 @@ +"""Tests for the closest-observation ingest builder (issue #509). + +Phase 1 — cover-driven epoch extraction: per-shard epochs from one or more +store roots' ``coverage.toc`` word-set covers (spec §10.5), union across +stores, AOI intersect, loud refusal when a store carries no readable cover. +Synthetic covers ride the same :func:`zagg.coverage_toc.build_cover_section` +producer the sweep uses; the committed golden ``coverage.toc`` fixture pins +the extraction against the frozen §10.5 grammar bytes. +""" + +import json +import logging +from pathlib import Path + +import numpy as np +import pytest +from mortie import from_datetime64, time2toc, to_datetime64 + +from zagg.catalog.closest_obs import ( + ReferenceEpochs, + _word_midpoints, + closest_obs_shardmap, + nearest_acquisitions, + reference_epochs, +) +from zagg.coverage_toc import ( + COVER_CAP, + COVER_NAME, + COVER_SPEC, + TEMPORAL_COVER_ORDER, + _encode_cover_block, + build_cover_section, + cover_words, + quantize_words, + read_cover, + write_cover, +) +from zagg.grids.morton import morton_word + +SPEC_DATA = Path(__file__).parent / "data" / "spec" +DAY_NS = 86_400 * 10**9 +#: An arbitrary but realistic base instant on the §8 internal-ns scale +#: (mirrors ``test_coverage_toc``'s convention). +BASE_NS = 5_344_000_000_000_000_000 +#: One order-18 cover bucket (2^45 ns) and half of it — the epoch-midpoint +#: error bound is half a BUCKET, not half a word (a word is a run of buckets). +BUCKET_NS = 2 ** (63 - TEMPORAL_COVER_ORDER) +HALF_BUCKET = np.timedelta64(BUCKET_NS // 2, "ns") + +#: The golden fixture's one shard, and the cell centre / far-away rings the +#: AOI tests use (order 4; centre from ``mortie.mort2geo``). +SHARD = "11213" +SHARD_KEY = morton_word(SHARD) +#: The shard next door (centre lat 14.54, lon 59.06) — the cross-shard join +#: tests need a second shard the single-shard fixture never exercised. +SHARD_B = "11212" +SHARD_B_KEY = morton_word(SHARD_B) + + +def _ring(lat, lon, half=2.0): + lats = np.array([lat - half, lat - half, lat + half, lat + half, lat - half]) + lons = np.array([lon - half, lon + half, lon + half, lon - half, lon - half]) + return [(lats, lons)] + + +AOI_AT_SHARD = _ring(14.54, 53.44) +AOI_ELSEWHERE = _ring(-40.0, -120.0) + + +def _write_store(root, shards: dict[str, np.ndarray], order: int = 4) -> str: + """A store root carrying a §10.5 cover claiming ``shards``' instants.""" + root = Path(root) + root.mkdir(parents=True, exist_ok=True) + contributions = { + decimal: [(None, None, None, quantize_words(time2toc(np.asarray(inst, dtype=np.uint64))))] + for decimal, inst in shards.items() + } + write_cover(str(root), build_cover_section(contributions, ["h"], order)) + return str(root) + + +def _instants(*day_offsets) -> np.ndarray: + return np.array([BASE_NS + d * DAY_NS for d in day_offsets], dtype=np.uint64) + + +def _utc(instants) -> np.ndarray: + return np.sort(np.asarray(to_datetime64(np.asarray(instants, np.uint64)), "datetime64[ns]")) + + +def _nearest_gap(mids: np.ndarray, true: np.ndarray) -> np.timedelta64: + """Largest distance from a true instant to its nearest epoch.""" + return np.abs(true[:, None] - mids[None, :]).min(axis=1).max() + + +class TestWordMidpoints: + def test_empty_words_decode_to_no_epochs(self): + out = _word_midpoints(np.empty(0, dtype=np.uint64)) + assert out.dtype == np.dtype("datetime64[ns]") and out.size == 0 + + def test_a_quantized_word_midpoint_stays_within_half_a_bucket(self): + inst = _instants(0, 5, 11) + words = quantize_words(time2toc(inst)) + mids = np.sort(_word_midpoints(words)) + true = np.sort(np.asarray(to_datetime64(inst), dtype="datetime64[ns]")) + assert mids.size == true.size + assert (np.abs(mids - true) <= HALF_BUCKET).all() + + def test_an_exact_timestamp_word_decodes_to_its_bucket_midpoint(self): + """An unquantized word still resolves to the bucket it falls in.""" + inst = _instants(3) + mids = _word_midpoints(np.asarray(time2toc(inst), dtype=np.uint64)) + assert mids.size == 1 + assert abs(mids[0] - _utc(inst)[0]) <= HALF_BUCKET + + def test_two_passes_one_bucket_apart_yield_two_epochs(self): + """``toc_normalize`` coalesces the abutting buckets into ONE word.""" + inst = np.array([BASE_NS, BASE_NS + BUCKET_NS], dtype=np.uint64) + words = quantize_words(time2toc(inst)) + assert len(words) == 1 # the coalescing that motivates the expansion + mids = np.sort(_word_midpoints(words)) + assert mids.size == 2 + assert _nearest_gap(mids, _utc(inst)) <= HALF_BUCKET + + def test_a_contiguous_campaign_yields_one_epoch_per_covered_bucket(self): + """A 10-day, 6-hourly campaign is one word — and 25 covered buckets.""" + inst = np.array([BASE_NS + i * 6 * 3600 * 10**9 for i in range(40)], dtype=np.uint64) + words = quantize_words(time2toc(inst)) + assert len(words) == 1 + covered = {int(t) >> (63 - TEMPORAL_COVER_ORDER) for t in inst} + mids = np.sort(_word_midpoints(words)) + assert mids.size == len(covered) == 25 + assert _nearest_gap(mids, _utc(inst)) <= HALF_BUCKET + + def test_a_pass_straddling_a_bucket_edge_yields_two_epochs(self): + """Benign over-selection: both epochs pick the same nearest granule.""" + edge = ((BASE_NS >> (63 - TEMPORAL_COVER_ORDER)) + 1) << (63 - TEMPORAL_COVER_ORDER) + inst = np.array([edge - 60 * 10**9, edge + 60 * 10**9], dtype=np.uint64) + words = quantize_words(time2toc(inst)) + assert len(words) == 1 + mids = np.sort(_word_midpoints(words)) + assert mids.size == 2 + assert _nearest_gap(mids, _utc(inst)) <= HALF_BUCKET + + def test_a_coarser_order_widens_the_buckets(self): + """The block's effective order drives the expansion, not the pin.""" + inst = np.array([BASE_NS, BASE_NS + BUCKET_NS], dtype=np.uint64) + coarse = quantize_words(time2toc(inst), TEMPORAL_COVER_ORDER - 2) + mids = _word_midpoints(coarse, TEMPORAL_COVER_ORDER - 2) + assert mids.size == 1 # both passes now share one 39 h bucket + assert _nearest_gap(mids, _utc(inst)) <= 4 * HALF_BUCKET + + +class TestReferenceEpochs: + def test_one_store_one_shard(self, tmp_path): + root = _write_store(tmp_path, {SHARD: _instants(0, 5, 11)}) + out = reference_epochs(root) + assert isinstance(out, ReferenceEpochs) + assert out.order == 4 + assert list(out.epochs) == [SHARD_KEY] + assert out.epochs[SHARD_KEY].size == 3 == out.total + assert out.stores == [root] + # Sorted unique datetime64[ns], the contract phase 2 searchsorted rides. + e = out.epochs[SHARD_KEY] + assert e.dtype == np.dtype("datetime64[ns]") + assert (np.diff(e.astype("int64")) > 0).all() + + def test_union_across_stores_is_deduplicated(self, tmp_path): + """Two stores quantized on the same grid: shared passes count once.""" + a = _write_store(tmp_path / "a", {SHARD: _instants(0, 5)}) + b = _write_store(tmp_path / "b", {SHARD: _instants(5, 11)}) + out = reference_epochs([a, b]) + assert out.epochs[SHARD_KEY].size == 3 + # Parity: the union equals each store's own epochs united. + ea = reference_epochs(a).epochs[SHARD_KEY] + eb = reference_epochs(b).epochs[SHARD_KEY] + assert np.array_equal(out.epochs[SHARD_KEY], np.union1d(ea, eb)) + + def test_partially_overlapping_covers_union_at_the_bucket_grid(self, tmp_path): + """Unequal words for a shared pass must not double it (bucket union). + + Three passes one bucket apart; store A saw passes 1+2 and store B + saw 2+3, so each store's cover carries a *different* two-bucket range + word for the shared pass. Word-level ``np.unique`` keeps both and + yields two displaced epochs for three passes; bucket-level union + yields exactly one epoch per covered bucket. + """ + passes = np.array([BASE_NS, BASE_NS + BUCKET_NS, BASE_NS + 2 * BUCKET_NS], np.uint64) + a = _write_store(tmp_path / "a", {SHARD: passes[:2]}) + b = _write_store(tmp_path / "b", {SHARD: passes[1:]}) + # Each store really does emit one coalesced (and unequal) word. + wa = cover_words(read_cover(a))[SHARD] + wb = cover_words(read_cover(b))[SHARD] + assert len(wa) == len(wb) == 1 and wa[0] != wb[0] + out = reference_epochs([a, b]) + e = out.epochs[SHARD_KEY] + assert e.size == 3 + assert _nearest_gap(e, _utc(passes)) <= HALF_BUCKET + + def test_a_shard_only_one_store_covers_still_contributes(self, tmp_path): + a = _write_store(tmp_path / "a", {SHARD: _instants(0)}) + b = _write_store(tmp_path / "b", {SHARD: _instants(40), "11212": _instants(7)}) + out = reference_epochs([a, b]) + assert set(out.epochs) == {SHARD_KEY, morton_word("11212")} + assert out.epochs[SHARD_KEY].size == 2 + + def test_a_store_without_a_cover_refuses_loudly(self, tmp_path): + (tmp_path / "empty").mkdir() + with pytest.raises(ValueError, match="cover-driven"): + reference_epochs(str(tmp_path / "empty")) + + def test_an_unknown_revision_cover_refuses_loudly(self, tmp_path): + root = tmp_path / "future" + root.mkdir() + (root / COVER_NAME).write_text( + json.dumps({"spec": "zagg-coverage-toc-cover/9", "shards": {}}) + ) + with pytest.raises(ValueError, match="unreadable"): + reference_epochs(str(root)) + + def test_a_shard_order_mismatch_between_stores_refuses(self, tmp_path): + a = _write_store(tmp_path / "a", {SHARD: _instants(0)}, order=4) + b = _write_store(tmp_path / "b", {SHARD: _instants(5)}, order=5) + with pytest.raises(ValueError, match="not.*comparable|shard order"): + reference_epochs([a, b]) + + def test_a_cover_without_a_shard_order_refuses_by_name(self, tmp_path): + """A missing ``order`` used to surface as an opaque ``TypeError``.""" + root = tmp_path / "orderless" + root.mkdir() + (root / COVER_NAME).write_text( + json.dumps({"spec": "zagg-coverage-toc-cover/1", "shards": {}}) + ) + with pytest.raises(ValueError, match="non-integer cover shard order"): + reference_epochs(str(root)) + + def test_no_stores_refuses(self): + with pytest.raises(ValueError, match="at least one"): + reference_epochs([]) + + def test_aoi_parts_restrict_the_shard_set(self, tmp_path): + root = _write_store(tmp_path, {SHARD: _instants(0, 5)}) + assert set(reference_epochs(root, aoi=AOI_AT_SHARD).epochs) == {SHARD_KEY} + assert reference_epochs(root, aoi=AOI_ELSEWHERE).epochs == {} + + def test_aoi_moc_restricts_the_shard_set(self, tmp_path): + from mortie import Moc + + root = _write_store(tmp_path, {SHARD: _instants(0, 5)}) + keep = Moc.from_polygon(*AOI_AT_SHARD[0]) + drop = Moc.from_polygon(*AOI_ELSEWHERE[0]) + assert set(reference_epochs(root, aoi=keep).epochs) == {SHARD_KEY} + assert reference_epochs(root, aoi=drop).epochs == {} + + def test_aoi_geojson_path_restricts_the_shard_set(self, tmp_path): + root = _write_store(tmp_path, {SHARD: _instants(0)}) + lats, lons = AOI_AT_SHARD[0] + geojson = { + "type": "Polygon", + "coordinates": [[[float(x), float(y)] for x, y in zip(lons, lats)]], + } + path = tmp_path / "aoi.geojson" + path.write_text(json.dumps(geojson)) + assert set(reference_epochs(root, aoi=str(path)).epochs) == {SHARD_KEY} + + +class TestCoarsenedBlock: + """§10.5 lets a block coarsen below the pin — the read half must be loud.""" + + #: Enough single-bucket claims (2 buckets apart, so no gap coalesces) to + #: blow the 512-word cap and force ``_cap_cover`` down an order. + COARSE_INSTANTS = np.array( + [BASE_NS + i * 2 * BUCKET_NS for i in range(COVER_CAP + 88)], dtype=np.uint64 + ) + + def _root(self, tmp_path): + return _write_store(tmp_path, {SHARD: self.COARSE_INSTANTS}) + + def test_the_block_really_coarsened(self, tmp_path): + root = self._root(tmp_path) + block = read_cover(root)["shards"][SHARD] + assert block["temporal_order"] < TEMPORAL_COVER_ORDER + + def test_a_coarsened_block_warns_and_reports_its_order(self, tmp_path, caplog): + root = self._root(tmp_path) + landed = read_cover(root)["shards"][SHARD]["temporal_order"] + with caplog.at_level("WARNING", logger="zagg.catalog.closest_obs"): + out = reference_epochs(root) + assert f"temporal order {landed}" in caplog.text + assert "below the pinned" in caplog.text + assert out.orders == {SHARD_KEY: landed} + assert out.tolerance(SHARD_KEY) == np.timedelta64(2 ** (62 - landed), "ns") + + def test_the_epochs_are_the_coarse_buckets_midpoints(self, tmp_path): + root = self._root(tmp_path) + landed = read_cover(root)["shards"][SHARD]["temporal_order"] + k = 63 - landed + internal = np.asarray(from_datetime64(reference_epochs(root).epochs[SHARD_KEY]), np.uint64) + # Every epoch sits at an aligned coarse-bucket midpoint, and at the + # coarse grid's buckets — not the pinned order's. + assert set(int(t) % 2**k for t in internal) == {2 ** (k - 1) - 1} + covered = {int(t) >> k for t in self.COARSE_INSTANTS} + assert set(int(t) >> k for t in internal) == covered + assert _nearest_gap( + reference_epochs(root).epochs[SHARD_KEY], _utc(self.COARSE_INSTANTS) + ) <= np.timedelta64(2 ** (k - 1), "ns") + + def test_an_uncoarsened_block_neither_warns_nor_hides_its_order(self, tmp_path, caplog): + root = _write_store(tmp_path, {SHARD: _instants(0, 5)}) + with caplog.at_level("WARNING", logger="zagg.catalog.closest_obs"): + out = reference_epochs(root) + assert "below the pinned" not in caplog.text + assert out.orders == {SHARD_KEY: TEMPORAL_COVER_ORDER} + assert out.tolerance(SHARD_KEY) == HALF_BUCKET + + +class TestGoldenFixture: + """The committed §7 ``temporal/`` fixture pins the frozen grammar bytes.""" + + #: The two epochs the committed fixture decodes to, pinned as literals + #: rather than recomputed with the function under test. Verified by hand: + #: each is the midpoint of an order-18 bucket (internal ns congruent to + #: 2^44 - 1 mod 2^45, buckets 151903 and 151915 — twelve apart, the + #: fixture's five-day gap), and each sits within half a bucket of the + #: generator's two campaign clusters (``TEMPORAL_BASE`` and +5 days). + GOLDEN = np.array( + ["2019-05-14T03:14:07.595891711", "2019-05-19T00:31:00.060957695"], + dtype="datetime64[ns]", + ) + + def test_the_golden_cover_decodes_to_the_pinned_epochs(self): + out = reference_epochs(str(SPEC_DATA / "temporal")) + assert out.order == 4 + assert list(out.epochs) == [SHARD_KEY] + assert np.array_equal(out.epochs[SHARD_KEY], self.GOLDEN) + assert out.orders == {SHARD_KEY: TEMPORAL_COVER_ORDER} + + def test_the_golden_epochs_are_order_18_bucket_midpoints(self): + """The arithmetic the ±4.9 h claim rests on, checked independently.""" + k = 63 - TEMPORAL_COVER_ORDER + internal = np.asarray(from_datetime64(self.GOLDEN), dtype=np.uint64) + assert [int(t) % 2**k for t in internal] == [2 ** (k - 1) - 1] * 2 + assert [int(t) >> k for t in internal] == [151903, 151915] + + def test_every_fixture_observation_has_an_epoch_within_half_a_bucket(self): + """The property the ruling actually claims, against the generator's + own recorded instants (``temporal.expected.json``), not against + anything :mod:`zagg.catalog.closest_obs` computed.""" + expected = json.loads((SPEC_DATA / "temporal.expected.json").read_text()) + true = np.array( + sorted({int(ns) for c in expected["cells"] for ns in c["obs_span_ns"]}), + dtype="datetime64[ns]", + ) + epochs = reference_epochs(str(SPEC_DATA / "temporal")).epochs[SHARD_KEY] + assert _nearest_gap(epochs, true) <= HALF_BUCKET + + def test_the_golden_epochs_sit_inside_the_fixture_campaign(self): + """Midpoints land inside the leaf's synthetic campaign window.""" + out = reference_epochs(str(SPEC_DATA / "temporal")) + lo = np.datetime64("2019-01-01") + hi = np.datetime64("2020-01-01") + e = out.epochs[SHARD_KEY] + assert ((e > lo) & (e < hi)).all() + + +class TestNearestAcquisitions: + """Phase 2 — the vectorized closest-1 selection core.""" + + def _dt(self, *hours): + return np.array( + [np.datetime64("2025-06-01T00:00") + np.timedelta64(h, "h") for h in hours] + ).astype("datetime64[ns]") + + def test_each_epoch_selects_its_nearest_acquisition(self): + times = self._dt(0, 10, 24) + epochs = self._dt(1, 9, 23) + sel, off = nearest_acquisitions(epochs, times) + assert sel.tolist() == [0, 1, 2] + assert off.astype("timedelta64[h]").astype(int).tolist() == [-1, 1, 1] + + def test_offsets_are_signed_acquisition_minus_epoch(self): + times = self._dt(12) + epochs = self._dt(10, 14) + sel, off = nearest_acquisitions(epochs, times) + assert sel.tolist() == [0, 0] + assert off[0] == np.timedelta64(2, "h") and off[1] == -np.timedelta64(2, "h") + + def test_an_exact_match_has_zero_offset(self): + times = self._dt(5, 7) + sel, off = nearest_acquisitions(self._dt(7), times) + assert sel.tolist() == [1] and off[0] == np.timedelta64(0, "ns") + + def test_a_tie_selects_the_earlier_acquisition(self): + times = self._dt(0, 10) + sel, off = nearest_acquisitions(self._dt(5), times) + assert sel.tolist() == [0] + assert off[0] == -np.timedelta64(5, "h") + + def test_selection_indices_refer_to_the_input_order(self): + times = self._dt(24, 0, 10) # unsorted catalog order + sel, _ = nearest_acquisitions(self._dt(23, 1), times) + assert sel.tolist() == [0, 1] + + def test_max_time_offset_boundary_exactly_at_selects(self): + times = self._dt(0) + sel, off = nearest_acquisitions(self._dt(6), times, max_time_offset=np.timedelta64(6, "h")) + assert sel.tolist() == [0] + + def test_max_time_offset_one_ns_past_drops_but_still_reports(self): + times = self._dt(0) + cap = np.timedelta64(6, "h") - np.timedelta64(1, "ns") + sel, off = nearest_acquisitions(self._dt(6), times, max_time_offset=cap) + assert sel.tolist() == [-1] + # The nearest offset is still reported — the loud record the drop rides. + assert off[0] == -np.timedelta64(6, "h") + + def test_no_acquisitions_selects_nothing_and_reports_nat(self): + sel, off = nearest_acquisitions(self._dt(1, 2), np.array([], dtype="datetime64[ns]")) + assert sel.tolist() == [-1, -1] + assert np.isnat(off).all() + + def test_no_epochs_is_empty(self): + sel, off = nearest_acquisitions(np.array([], dtype="datetime64[ns]"), self._dt(0)) + assert sel.size == 0 and off.size == 0 + + def test_a_negative_max_time_offset_refuses(self): + with pytest.raises(ValueError, match="non-negative"): + nearest_acquisitions(self._dt(1), self._dt(0), max_time_offset=-np.timedelta64(1, "h")) + + def test_duplicate_acquisition_times_select_the_first_record(self): + # One datatake, three granules at the same instant, plus a far one. + times = self._dt(0, 0, 0, 96) + for epoch in (self._dt(-1), self._dt(0), self._dt(1)): + sel, _ = nearest_acquisitions(epoch, times) + assert sel.tolist() == [0] + + def _check_oracle(self, epochs, times, cap): + """Brute force: nearest, ties to the earlier, then the lowest index.""" + sel, off = nearest_acquisitions(epochs, times, max_time_offset=cap) + t = times.astype("int64") + for k, e in enumerate(epochs.astype("int64")): + d = np.abs(t - e) + best = np.flatnonzero(d == d.min()) + # tie -> earlier acquisition; equal times -> lowest catalog index + want = best[np.argmin(t[best])] if best.size > 1 else best[0] + assert off[k] == np.timedelta64(int(t[want] - e), "ns") + if d.min() <= cap.astype("timedelta64[ns]").astype("int64"): + assert sel[k] == want + else: + assert sel[k] == -1 + + def test_an_unconvertible_max_time_offset_refuses_by_name(self): + with pytest.raises(ValueError, match="does not convert exactly to nanoseconds"): + nearest_acquisitions( + self._dt(1), self._dt(0), max_time_offset=np.timedelta64(1000, "Y") + ) + + def test_a_nat_max_time_offset_refuses(self): + with pytest.raises(ValueError, match="must be a real duration"): + nearest_acquisitions(self._dt(1), self._dt(0), max_time_offset=np.timedelta64("NaT")) + + def test_a_gap_past_int64_nanoseconds_never_pairs_under_a_cap(self): + # 584 years apart — a real distance that timedelta64[ns] cannot carry. + times = np.array(["1677-09-22T00:12:44"], dtype="datetime64[ns]") + epochs = np.array(["2262-04-11T23:47:16"], dtype="datetime64[ns]") + for cap in (np.timedelta64(1, "D"), np.timedelta64(120, "D")): + sel, off = nearest_acquisitions(epochs, times, max_time_offset=cap) + assert sel.tolist() == [-1] + assert np.isnat(off).all() + # Uncapped, the nearest is still the nearest; the offset saturates to NaT. + sel, off = nearest_acquisitions(epochs, times) + assert sel.tolist() == [0] + assert np.isnat(off).all() + + def test_a_two_century_gap_still_reports_an_exact_offset(self): + times = np.array(["1800-01-01", "2150-01-01"], dtype="datetime64[ns]") + epochs = np.array(["2000-01-01"], dtype="datetime64[ns]") + sel, off = nearest_acquisitions(epochs, times) + assert sel.tolist() == [1] + assert off[0] == times[1] - epochs[0] + + def test_a_nat_acquisition_time_refuses(self): + times = np.array(["2025-06-01T00", "NaT", "2025-06-01T10"], dtype="datetime64[ns]") + with pytest.raises(ValueError, match="times carries NaT"): + nearest_acquisitions(self._dt(1, 9), times) + + def test_a_nat_epoch_refuses(self): + epochs = np.array(["2025-06-01T01", "NaT"], dtype="datetime64[ns]") + with pytest.raises(ValueError, match="epochs carries NaT"): + nearest_acquisitions(epochs, self._dt(0, 10)) + + def test_matches_a_brute_force_oracle(self): + rng = np.random.default_rng(7) + base = np.datetime64("2025-01-01").astype("datetime64[ns]").astype("int64") + times = (base + rng.integers(0, 400 * 86_400, 60) * 10**9).astype("datetime64[ns]") + epochs = (base + rng.integers(0, 400 * 86_400, 45) * 10**9).astype("datetime64[ns]") + self._check_oracle(epochs, times, np.timedelta64(2, "D")) + + def test_matches_the_oracle_with_every_acquisition_time_duplicated(self): + rng = np.random.default_rng(7) + base = np.datetime64("2025-01-01").astype("datetime64[ns]").astype("int64") + raw = base + rng.integers(0, 400 * 86_400, 20) * 10**9 + # Every instant carried by two catalog records, interleaved out of order. + times = np.concatenate([raw, raw]).astype("datetime64[ns]") + epochs = (base + rng.integers(0, 400 * 86_400, 30) * 10**9).astype("datetime64[ns]") + self._check_oracle(epochs, times, np.timedelta64(2, "D")) + + +# ── phase 3: the builder ───────────────────────────────────────────────────── +# +# Geometry: everything happens in shard 11213 (order 4; centre lat 14.54, +# lon 53.44 from mortie.mort2geo) — the same shard the golden fixture covers — +# with S2-like STAC items (multi-band assets, NO canonical data asset, a +# per-item datetime) footprinted around the centre so the spatial build +# assigns them there. + + +def _s2_item(gid, iso, lat=14.54, lon=53.44, half=0.4): + ring = [ + [lon - half, lat - half], + [lon + half, lat - half], + [lon + half, lat + half], + [lon - half, lat + half], + [lon - half, lat - half], + ] + return { + "type": "Feature", + "stac_version": "1.0.0", + "id": gid, + "geometry": {"type": "Polygon", "coordinates": [ring]}, + "bbox": [lon - half, lat - half, lon + half, lat + half], + "properties": {"datetime": iso}, + "collection": "sentinel-2-l2a", + "stac_extensions": [], + "links": [], + "assets": { + "red": {"href": f"https://h/{gid}/B04.tif", "roles": ["data"]}, + "nir": {"href": f"https://h/{gid}/B08.tif", "roles": ["data"]}, + }, + } + + +def _s2_catalog(items, bbox=(52.0, 13.0, 55.0, 16.0)): + import pyarrow as pa + import stac_geoparquet.arrow as sga + + from zagg.catalog.sources import Catalog + + return Catalog( + pa.table(sga.parse_stac_items_to_arrow(items)), + {"collection": "sentinel-2-l2a", "bbox": list(bbox)}, + ) + + +def _grid(parent_order=4): + from zagg.grids import HealpixGrid + + return HealpixGrid(parent_order, 6) + + +def _epoch_iso(day): + """An ISO instant near BASE_NS + day*DAY (the synthetic covers' passes).""" + ns = np.datetime64(to_datetime64(np.array([BASE_NS + day * DAY_NS], dtype=np.uint64))[0], "ns") + return np.datetime_as_string(ns.astype("datetime64[s]")) + "Z" + + +class TestClosestObsShardmap: + def _setup(self, tmp_path, days=(0, 5, 11), s2_days=(0.1, 5.2, 20.0)): + store = _write_store(tmp_path / "ref", {SHARD: _instants(*days)}) + items = [_s2_item(f"S2_{i}", _epoch_iso(d)) for i, d in enumerate(s2_days)] + return store, _s2_catalog(items) + + def test_builds_a_standard_shardmap(self, tmp_path): + store, cat = self._setup(tmp_path) + sm = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + assert sm.shard_keys == [SHARD_KEY] + ids = {g["id"] for g in sm.granules[0]} + # Epochs at days 0/5/11 pick S2_0 (0.1), S2_1 (5.2), S2_2 (20 vs 5.2: + # day 11 is 5.8 days from S2_1 and 9 days from S2_2 -> S2_1), deduped. + assert ids == {"S2_0", "S2_1"} + assert sm.metadata["total_pairs"] == 2 + assert sm.metadata["granules_assigned"] == 2 + assert sm.metadata["closest_obs"]["epochs_paired"] == 3 + assert sm.metadata["closest_obs"]["epochs_dropped"] == 0 + + def test_provenance_rows_reconstruct_the_pairing(self, tmp_path): + store, cat = self._setup(tmp_path) + sm = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + by_id = {g["id"]: g for g in sm.granules[0]} + # S2_1 was selected by two epochs (days 5 and 11): dedupe keeps ONE + # entry carrying BOTH provenance rows, row-aligned. + assert len(by_id["S2_1"]["paired_epochs"]) == 2 + assert len(by_id["S2_1"]["epoch_offsets_ns"]) == 2 + # Signed offsets: acquisition - epoch. Day-5 epoch precedes the day-5.2 + # acquisition -> positive; day-11 epoch follows it -> negative. + offs = sorted(by_id["S2_1"]["epoch_offsets_ns"]) + assert offs[0] < 0 < offs[1] + # And each is within half a bucket + the true separation. + assert len(by_id["S2_0"]["paired_epochs"]) == 1 + + def test_the_map_json_round_trips_with_provenance(self, tmp_path): + store, cat = self._setup(tmp_path) + sm = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + path = tmp_path / "map.json" + sm.to_json(str(path)) + from zagg.catalog.shardmap import ShardMap + + back = ShardMap.from_json(str(path)) + assert back.shard_keys == sm.shard_keys + assert back.granules == sm.granules + assert back.metadata["closest_obs"] == sm.metadata["closest_obs"] + + def test_max_time_offset_drops_are_recorded_loudly(self, tmp_path, caplog): + # Day-11 epoch's nearest acquisition is 5.8 days away: beyond a 2-day + # cap it selects nothing, and the drop carries the near-miss offset. + store, cat = self._setup(tmp_path) + with caplog.at_level(logging.WARNING, logger="zagg.catalog.closest_obs"): + sm = closest_obs_shardmap( + cat, + store, + grid=_grid(), + backend="mortie", + max_time_offset=np.timedelta64(2, "D"), + ) + rec = sm.metadata["closest_obs"] + assert rec["epochs_dropped"] == 1 and len(rec["dropped"]) == 1 + assert rec["epochs_total"] == rec["epochs_paired"] + rec["epochs_dropped"] + assert rec["dropped"][0]["shard"] == SHARD + assert rec["dropped"][0]["nearest_offset_ns"] is not None + assert any("selected nothing" in m for m in caplog.messages) + ids = {g["id"] for g in sm.granules[0]} + assert ids == {"S2_0", "S2_1"} + + def test_filtered_map_is_a_subset_of_the_spatial_map(self, tmp_path): + from zagg.catalog.shardmap import ShardMap + + store, cat = self._setup(tmp_path) + spatial = ShardMap.build(cat, _grid(), backend="mortie") + sm = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + spatial_pairs = { + (k, g["id"]) for k, gr in zip(spatial.shard_keys, spatial.granules) for g in gr + } + paired = {(k, g["id"]) for k, gr in zip(sm.shard_keys, sm.granules) for g in gr} + assert paired <= spatial_pairs + + def test_grid_order_mismatch_refuses(self, tmp_path): + store, cat = self._setup(tmp_path) + with pytest.raises(ValueError, match="parent_order"): + closest_obs_shardmap(cat, store, grid=_grid(parent_order=5), backend="mortie") + + def test_max_granules_per_shard_refuses_loudly(self, tmp_path): + store, cat = self._setup(tmp_path) + with pytest.raises(ValueError, match="max_granules_per_shard"): + closest_obs_shardmap( + cat, store, grid=_grid(), backend="mortie", max_granules_per_shard=1 + ) + + def test_estimate_reports_without_building(self, tmp_path): + store, cat = self._setup(tmp_path) + est = closest_obs_shardmap( + cat, + store, + grid=_grid(), + backend="mortie", + estimate=True, + max_granules_per_shard=1, + bytes_per_granule=10**6, + ) + assert isinstance(est, dict) + assert est["shards"] == 1 and est["granules"] == 2 and est["pairs"] == 2 + assert est["per_shard"] == {SHARD: 2} + assert est["histogram"] == {2: 1} + assert est["est_bytes"] == 2 * 10**6 + assert est["max_cost_usd"] > 0 + # The cost gate REPORTS violations in a dry run instead of raising. + assert est["violations"] == [(SHARD, 2)] + + def test_a_granule_without_acquisition_time_refuses(self, tmp_path): + store, _ = self._setup(tmp_path) + item = _s2_item("S2_bare", _epoch_iso(0)) + del item["properties"]["datetime"] # -> null datetime column + # stac-geoparquet requires a datetime key; set None explicitly instead + item["properties"]["datetime"] = None + cat = _s2_catalog([item]) + with pytest.raises(ValueError, match="acquisition time"): + closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + + def test_the_spatial_aoi_mask_is_neither_carried_nor_claimed(self, tmp_path): + # The strict-AOI payload (#101) belongs to the spatial map; the derived + # map carries none, so its metadata must not advertise one either. + from zagg.catalog.shardmap import ShardMap + from zagg.config import default_config + from zagg.grids import HealpixGrid + + cfg = default_config("atl06") + cfg.output = {**cfg.output, "aoi_mask": True} + grid = HealpixGrid(4, 6, config=cfg) + store, cat = self._setup(tmp_path) + spatial = ShardMap.build(cat, grid, backend="mortie") + assert spatial.aoi_mask is not None and spatial.metadata["aoi_mask"] is True + sm = closest_obs_shardmap(cat, store, grid=grid, backend="mortie") + assert sm.aoi_mask is None + assert "aoi_mask" not in sm.metadata + + def test_a_shard_the_catalog_never_reaches_ledgers_its_epochs(self, tmp_path, caplog): + # Two-shard cover, one-shard catalog: the unreached shard's epochs are + # dropped ROWS, not silence — epochs_total reconciles by construction. + store = _write_store( + tmp_path / "ref", {SHARD: _instants(0, 5, 11), SHARD_B: _instants(0, 5, 11, 17)} + ) + items = [_s2_item(f"S2_{i}", _epoch_iso(d)) for i, d in enumerate((0.1, 5.2, 20.0))] + with caplog.at_level(logging.WARNING, logger="zagg.catalog.closest_obs"): + sm = closest_obs_shardmap(_s2_catalog(items), store, grid=_grid(), backend="mortie") + rec = sm.metadata["closest_obs"] + assert sm.shard_keys == [SHARD_KEY] + assert rec["shards_without_acquisitions"] == [SHARD_B] + assert rec["epochs_paired"] == 3 + assert rec["epochs_total"] == rec["epochs_paired"] + rec["epochs_dropped"] + rows = [d for d in rec["dropped"] if d["shard"] == SHARD_B] + assert len(rows) == rec["epochs_dropped"] > 0 + assert all(d["nearest_offset_ns"] is None for d in rows) + assert any("NO spatially-assigned acquisitions" in m for m in caplog.messages) + + def test_a_start_datetime_only_record_still_emits_a_datetime(self, tmp_path): + # STAC's null-datetime + start/end_datetime form pairs on ``time_start``; + # the emitted entry must carry the ``datetime`` raster dispatch reads. + store, _ = self._setup(tmp_path) + item = _s2_item("S2_range", _epoch_iso(0)) + item["properties"]["datetime"] = None + item["properties"]["start_datetime"] = _epoch_iso(0) + item["properties"]["end_datetime"] = _epoch_iso(1) + cat = _s2_catalog([item]) + assert "datetime" not in cat.granule_records()[0] # the fixture really lacks it + sm = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + entry = sm.granules[0][0] + assert entry["datetime"] == entry["time_start"] + + def _two_shard(self, tmp_path): + """Cover + catalog spanning 11213 and its lon neighbour 11212.""" + store = _write_store( + tmp_path / "ref", {SHARD: _instants(0, 5, 11), SHARD_B: _instants(0, 5, 11)} + ) + items = [_s2_item(f"A{i}", _epoch_iso(d)) for i, d in enumerate((0.1, 5.2))] + items += [_s2_item(f"B{i}", _epoch_iso(d), lon=59.06) for i, d in enumerate((0.1, 5.2))] + return store, _s2_catalog(items, bbox=(51.0, 12.0, 61.0, 18.0)) + + def _spy_build(self, monkeypatch): + """Record the kwargs the builder hands ``ShardMap.build``.""" + from zagg.catalog.shardmap import ShardMap + + calls = [] + real = ShardMap.build.__func__ + + def spy(cls, catalog, grid, **kw): + calls.append(kw) + return real(cls, catalog, grid, **kw) + + monkeypatch.setattr(ShardMap, "build", classmethod(spy)) + return calls + + def _pairs(self, sm): + return {(k, g["id"]) for k, gr in zip(sm.shard_keys, sm.granules) for g in gr} + + def test_an_aoi_excluded_shard_is_not_coverage_disagreement(self, tmp_path, monkeypatch): + # A Moc aoi has no ring-parts form, so the spatial build stays unscoped + # and still assigns 11212 — whose epochs the aoi clipped away. That is + # self-inflicted, not "the reference stores never observed this ground". + from mortie import Moc + + calls = self._spy_build(monkeypatch) + store, cat = self._two_shard(tmp_path) + sm = closest_obs_shardmap( + cat, store, grid=_grid(), backend="mortie", aoi=Moc.from_polygon(*AOI_AT_SHARD[0]) + ) + assert calls[0]["region"] is None + assert sm.shard_keys == [SHARD_KEY] + assert sm.metadata["closest_obs"]["spatial_shards_without_epochs"] == 0 + + def test_ring_parts_aoi_scopes_the_spatial_build(self, tmp_path, monkeypatch): + calls = self._spy_build(monkeypatch) + store, cat = self._two_shard(tmp_path) + full = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + scoped = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie", aoi=AOI_AT_SHARD) + assert calls[0]["region"] is None # aoi=None is unchanged + assert calls[1]["region"] is AOI_AT_SHARD + # Scoping the intersection changes the cost, never the answer. + assert SHARD_B_KEY in full.shard_keys + assert self._pairs(scoped) == {p for p in self._pairs(full) if p[0] == SHARD_KEY} + + def test_a_geojson_aoi_scopes_the_spatial_build_as_parts(self, tmp_path, monkeypatch): + calls = self._spy_build(monkeypatch) + store, cat = self._two_shard(tmp_path) + lats, lons = AOI_AT_SHARD[0] + path = tmp_path / "aoi.geojson" + path.write_text( + json.dumps( + { + "type": "Polygon", + "coordinates": [[[float(x), float(y)] for x, y in zip(lons, lats)]], + } + ) + ) + sm = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie", aoi=str(path)) + # Resolved to ring parts before it reaches the build, never the path. + assert isinstance(calls[0]["region"], list) + assert sm.shard_keys == [SHARD_KEY] + + def test_reprojecting_a_paired_map_drops_the_provenance(self, tmp_path): + # Pins the documented trap: ShardMap._granule_entry does not know the + # two provenance keys, so even the same-order noop branch strips them + # while metadata["closest_obs"] rides through describing the SOURCE + # map. Carrying them through reproject is a shardmap.py change left + # standing for review — until then, rebuild, never reproject. + store, cat = self._setup(tmp_path) + sm = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + assert "paired_epochs" in sm.granules[0][0] + noop = sm.reproject(_grid()) + assert "paired_epochs" not in noop.granules[0][0] + assert "epoch_offsets_ns" not in noop.granules[0][0] + assert noop.metadata["closest_obs"] == sm.metadata["closest_obs"] + + +# ── phase 4: two-store scenarios ───────────────────────────────────────────── +# +# An ATL03-like sparse cover beside a GEDI-like denser cover with an interior +# gap, both over shard 11213; the GEDI-like store also claims 11212, where the +# S2 catalog has no acquisitions at all. S2 revisit ~4.3 days across the span. + + +class TestTwoStoreScenarios: + # Sparse, nothing near the middle. Day 13 is what makes A's selection + # NOT a subset of B's (it reaches S2_3, which no B epoch does), so the + # union-parity test below can actually fail on a dropped store; it is + # still clear of the gap set (its S2_3 sits at day 12.9, and the first + # gap acquisition at 17.2 d is 4.2 d away, past the 3 d cap). + A_DAYS = (0, 13, 55) + B_DAYS = (0, 2, 4, 6, 8, 10, 50, 52, 54, 56, 58, 60) # dense, gap 10..50 + + def _stores(self, tmp_path): + a = _write_store(tmp_path / "atl03", {SHARD: _instants(*self.A_DAYS)}) + b = _write_store( + tmp_path / "gedi", + {SHARD: _instants(*self.B_DAYS), "11212": _instants(1, 3)}, + ) + return a, b + + def _catalog(self, days=None): + days = days if days is not None else [d / 10 for d in range(0, 600, 43)] + return _s2_catalog([_s2_item(f"S2_{i}", _epoch_iso(d)) for i, d in enumerate(days)]) + + def test_an_acquisition_in_the_cover_gap_is_never_selected(self, tmp_path): + """Days 10..50 are a gap in BOTH stores: no epoch reaches into it.""" + a, b = self._stores(tmp_path) + cat = self._catalog() + sm = closest_obs_shardmap( + cat, + [a, b], + grid=_grid(), + backend="mortie", + max_time_offset=np.timedelta64(3, "D"), + ) + ids = {g["id"] for g in sm.granules[sm.shard_keys.index(SHARD_KEY)]} + # Acquisitions land every 4.3 d; those in (13, 47) days sit >3 d from + # every epoch (epochs live at passes 0..10 and 50..60) so none of them + # may appear -- the cover gap prunes them even though they are + # spatially assigned. + gap = {f"S2_{i}" for i, d in enumerate(d / 10 for d in range(0, 600, 43)) if 13 < d < 47} + assert gap and not (ids & gap) + # And the near-gap acquisitions ARE selected -- pinned by name, since + # a bare `assert ids` passes on any non-empty selection and so cannot + # fail when the shoulders stop reaching. S2_2 (day 8.6) is reached by + # the epochs at days 8.098/10.134, S2_12 (day 51.6) by 50.042/52.078. + assert {"S2_2", "S2_12"} <= ids # the gap's shoulder epochs reach across + + def test_union_parity_across_stores(self, tmp_path): + """map(A ∪ B) selects exactly union(map(A), map(B)) per shard. + + Closest-1 is per-epoch independent and epochs are the union across + stores, so the selection commutes with the union. + """ + a, b = self._stores(tmp_path) + cat = self._catalog() + kw = dict(grid=_grid(), backend="mortie", max_time_offset=np.timedelta64(3, "D")) + + def _pairs(sm): + return {(k, g["id"]) for k, gr in zip(sm.shard_keys, sm.granules) for g in gr} + + both = _pairs(closest_obs_shardmap(cat, [a, b], **kw)) + only_a = _pairs(closest_obs_shardmap(cat, a, **kw)) + only_b = _pairs(closest_obs_shardmap(cat, b, **kw)) + assert both == only_a | only_b + # Both stores must CONTRIBUTE, or parity is satisfied by a builder that + # keeps only one of them: pin each side as a proper subset of the union + # so neither store's epochs can be silently dropped. + assert only_a < both and only_b < both + + def test_the_filtered_map_is_a_subset_of_the_spatial_map(self, tmp_path): + from zagg.catalog.shardmap import ShardMap + + a, b = self._stores(tmp_path) + cat = self._catalog() + spatial = ShardMap.build(cat, _grid(), backend="mortie") + sm = closest_obs_shardmap(cat, [a, b], grid=_grid(), backend="mortie") + spatial_pairs = { + (k, g["id"]) for k, gr in zip(spatial.shard_keys, spatial.granules) for g in gr + } + got = {(k, g["id"]) for k, gr in zip(sm.shard_keys, sm.granules) for g in gr} + assert got and got <= spatial_pairs + + def test_a_covered_shard_with_no_acquisitions_is_recorded(self, tmp_path): + a, b = self._stores(tmp_path) + cat = self._catalog() # items only around shard 11213's centre + sm = closest_obs_shardmap(cat, [a, b], grid=_grid(), backend="mortie") + rec = sm.metadata["closest_obs"] + assert "11212" in rec["shards_without_acquisitions"] + # Membership alone would pass while the epochs went unledgered; the + # shard's two cover epochs (days 1 and 3) must each show up in + # ``dropped`` with no offset -- there is no acquisition to measure. + rows = [d for d in rec["dropped"] if d["shard"] == "11212"] + assert len(rows) == 2 + assert all(d["nearest_offset_ns"] is None for d in rows) + got = np.sort(np.array([d["epoch"] for d in rows], dtype="datetime64[ns]")) + assert _nearest_gap(got, _utc(_instants(1, 3))) <= HALF_BUCKET + assert morton_word("11212") not in sm.shard_keys + + def test_offset_boundary_exactly_at_selects_one_ns_past_drops(self, tmp_path): + """The cap boundary at BUILDER level, against a cover-derived epoch.""" + store = _write_store(tmp_path / "ref", {SHARD: _instants(0)}) + epochs = reference_epochs(store).epochs[SHARD_KEY] + assert epochs.size == 1 + acq_iso = _epoch_iso(2.0) + cat = _s2_catalog([_s2_item("S2_only", acq_iso)]) + acq = np.datetime64(acq_iso.rstrip("Z")).astype("datetime64[ns]") + exact = acq - epochs[0] # signed timedelta64[ns], acquisition - epoch + assert exact > np.timedelta64(0, "ns") + kw = dict(grid=_grid(), backend="mortie") + at = closest_obs_shardmap(cat, store, max_time_offset=exact, **kw) + assert [g["id"] for g in at.granules[0]] == ["S2_only"] + assert at.granules[0][0]["epoch_offsets_ns"] == [int(exact.astype("int64"))] + past = closest_obs_shardmap( + cat, store, max_time_offset=exact - np.timedelta64(1, "ns"), **kw + ) + assert past.shard_keys == [] + rec = past.metadata["closest_obs"] + assert rec["epochs_dropped"] == 1 + assert rec["dropped"][0]["nearest_offset_ns"] == int(exact.astype("int64")) + + +# ── espg tolerance ruling (2026-08-24, thread r3845481805) ─────────────────── + + +def _write_store_at_order(root, decimal, instants, block_order, order=4): + """A store whose cover block sits EXPLICITLY at ``block_order`` < the pin.""" + root = Path(root) + root.mkdir(parents=True, exist_ok=True) + words = quantize_words(time2toc(np.asarray(instants, dtype=np.uint64)), block_order) + section = { + "spec": COVER_SPEC, + "source": "test", + "order": order, + "temporal_order": TEMPORAL_COVER_ORDER, + "cap": COVER_CAP, + "fields": ["h"], + "element": {"dtype": "uint64", "shape": [-1]}, + "encoding": "base64", + "shards": {decimal: _encode_cover_block(np.asarray(words, np.uint64), block_order)}, + } + (root / COVER_NAME).write_text(json.dumps(section)) + return str(root) + + +class TestCoarsenedCoverTolerance: + """max_time_offset is a precision bar on the epochs, not just the offsets. + + The ruling of record: with a cap set, an epoch whose bucket half-span + exceeds it cannot be paired to the stated precision and drops loudly as + its own ledger category; with no cap, one warning names the effective + resolution and pairing proceeds (widening is lawful, §10.5). Flat-warn + risks silently arbitrary pairings from a cap-degraded store; flat-refuse + fails whole builds over blocks that may not even intersect the AOI. + """ + + #: An order-12 bucket's half-span: 2^50 ns ≈ 13.03 days. + ORDER12_HALF = np.timedelta64(2**50, "ns") + + def _coarse_store(self, tmp_path, days=(40.0,), name="coarse"): + return _write_store_at_order(tmp_path / name, SHARD, _instants(*days), 12) + + def test_low_resolution_epochs_drop_under_a_cap(self, tmp_path): + store = self._coarse_store(tmp_path) + cat = _s2_catalog([_s2_item("S2_0", _epoch_iso(40.0))]) + sm = closest_obs_shardmap( + cat, store, grid=_grid(), backend="mortie", max_time_offset=np.timedelta64(3, "D") + ) + rec = sm.metadata["closest_obs"] + # Every epoch in the order-12 block is unresolvable at ±3 d. + assert sm.shard_keys == [] + assert rec["epochs_total"] > 0 + assert rec["epochs_dropped_low_resolution"] == rec["epochs_dropped"] == rec["epochs_total"] + assert rec["epochs_total"] == rec["epochs_paired"] + rec["epochs_dropped"] + row = rec["dropped"][0] + assert row["temporal_order"] == 12 + assert row["cover_half_span_ns"] == 2**50 + # Its own category: no near-miss offset — the selection never ran. + assert "nearest_offset_ns" not in row + + def test_the_estimate_reports_the_category(self, tmp_path): + store = self._coarse_store(tmp_path) + cat = _s2_catalog([_s2_item("S2_0", _epoch_iso(40.0))]) + est = closest_obs_shardmap( + cat, + store, + grid=_grid(), + backend="mortie", + estimate=True, + max_time_offset=np.timedelta64(3, "D"), + ) + assert est["epochs_dropped_low_resolution"] == est["epochs_total"] > 0 + assert est["dropped"][0]["temporal_order"] == 12 + + def test_no_cap_warns_once_and_pairs_everything(self, tmp_path, caplog): + store = self._coarse_store(tmp_path) + cat = _s2_catalog([_s2_item("S2_0", _epoch_iso(40.0))]) + with caplog.at_level(logging.WARNING, logger="zagg.catalog.closest_obs"): + sm = closest_obs_shardmap(cat, store, grid=_grid(), backend="mortie") + build_warnings = [m for m in caplog.messages if "no max_time_offset" in m] + assert len(build_warnings) == 1 + assert "temporal order 12" in build_warnings[0] + rec = sm.metadata["closest_obs"] + assert rec["epochs_dropped_low_resolution"] == 0 + assert rec["epochs_paired"] == rec["epochs_total"] > 0 + assert {g["id"] for g in sm.granules[0]} == {"S2_0"} + + def test_the_cap_arm_warns_that_a_coarsened_cover_dropped_them(self, tmp_path, caplog): + """Symmetric to the no-cap arm: the arm that DISCARDS is at least as loud. + + The no-cap arm names the effective resolution and proceeds; the + with-cap arm throws epochs away, so it gets its own line naming the + count and the cause. The generic summary must stop mis-attributing + these rows to distance or a catalog gap — neither ran for them. + """ + store = self._coarse_store(tmp_path) + cat = _s2_catalog([_s2_item("S2_0", _epoch_iso(40.0))]) + with caplog.at_level(logging.WARNING, logger="zagg.catalog.closest_obs"): + sm = closest_obs_shardmap( + cat, store, grid=_grid(), backend="mortie", max_time_offset=np.timedelta64(3, "D") + ) + n = sm.metadata["closest_obs"]["epochs_dropped_low_resolution"] + assert n > 0 + drops = [m for m in caplog.messages if "UNRESOLVABLE" in m] + assert len(drops) == 1 + assert f"{n} epoch(s) dropped as UNRESOLVABLE" in drops[0] + assert "a coarsened cover, NOT distance" in drops[0] + assert "temporal order 12" in drops[0] + # ...and the generic summary now names the third cause. + summary = [m for m in caplog.messages if "selected nothing" in m] + assert len(summary) == 1 + assert "cover block too coarse for the offset" in summary[0] + # The no-cap line is the OTHER arm's; it must not fire under a cap. + assert not [m for m in caplog.messages if "no max_time_offset" in m] + + def test_half_span_exactly_at_the_cap_still_pairs(self, tmp_path): + """The pinned boundary side: half-span == max_time_offset is pairable. + + Strictly-greater drops, matching the selection gate where an offset + exactly at the cap SELECTS — both boundaries sit on the permissive + side. (The acquisition inside the epoch's own bucket is within + half-span of its midpoint, so the selection also passes.) + """ + store = self._coarse_store(tmp_path) + cat = _s2_catalog([_s2_item("S2_0", _epoch_iso(40.0))]) + sm = closest_obs_shardmap( + cat, store, grid=_grid(), backend="mortie", max_time_offset=self.ORDER12_HALF + ) + rec = sm.metadata["closest_obs"] + assert rec["epochs_dropped_low_resolution"] == 0 + assert rec["epochs_paired"] == rec["epochs_total"] > 0 + one_ns_under = self.ORDER12_HALF - np.timedelta64(1, "ns") + dropped = closest_obs_shardmap( + cat, store, grid=_grid(), backend="mortie", max_time_offset=one_ns_under + ) + assert ( + dropped.metadata["closest_obs"]["epochs_dropped_low_resolution"] == rec["epochs_total"] + ) + + def test_mixed_orders_drop_only_the_coarse_epochs(self, tmp_path): + """One shard, a pinned store beside a coarsened one: per-epoch gating.""" + fine = _write_store(tmp_path / "fine", {SHARD: _instants(0, 5)}) + coarse = self._coarse_store(tmp_path) + ref = reference_epochs([fine, coarse]) + assert set(ref.epoch_orders[SHARD_KEY].tolist()) == {TEMPORAL_COVER_ORDER, 12} + n_coarse = int((ref.epoch_orders[SHARD_KEY] == 12).sum()) + assert n_coarse > 0 + cat = _s2_catalog( + [_s2_item(f"S2_{i}", _epoch_iso(d)) for i, d in enumerate((0.1, 5.2, 40.0))] + ) + sm = closest_obs_shardmap( + cat, + [fine, coarse], + grid=_grid(), + backend="mortie", + max_time_offset=np.timedelta64(3, "D"), + ) + rec = sm.metadata["closest_obs"] + # Only the coarse-block epochs fail the precision bar; the pinned + # store's epochs pair through the same shard. + assert rec["epochs_dropped_low_resolution"] == n_coarse + assert rec["epochs_total"] == rec["epochs_paired"] + rec["epochs_dropped"] + ids = {g["id"] for g in sm.granules[0]} + assert ids == {"S2_0", "S2_1"} + + def test_an_unreached_shard_still_counts_its_resolution_drops(self, tmp_path): + """The category is a property of the epoch, not of the catalog's reach. + + Unresolvability follows from the epoch's own cover block, so it must + not flip to a no-acquisition row just because the raster catalog + never reaches the shard — an operator sizing ``max_time_offset`` off + ``epochs_dropped_low_resolution == 0`` would fix the catalog gap and + watch the same epochs reappear as resolution drops. + """ + store = self._coarse_store(tmp_path) + cat = _s2_catalog( + [_s2_item("S2_far", _epoch_iso(40.0), lat=-40.0, lon=-120.0)], + bbox=(-121.0, -41.0, -119.0, -39.0), + ) + sm = closest_obs_shardmap( + cat, store, grid=_grid(), backend="mortie", max_time_offset=np.timedelta64(3, "D") + ) + rec = sm.metadata["closest_obs"] + assert rec["epochs_dropped_low_resolution"] == rec["epochs_total"] > 0 + assert all("temporal_order" in d for d in rec["dropped"]) + # The shard is still named: that row is about the SHARD, not its epochs. + assert rec["shards_without_acquisitions"] == [SHARD] + assert rec["epochs_total"] == rec["epochs_paired"] + rec["epochs_dropped"] + + def test_an_unreached_shard_splits_coarse_from_surviving_epochs(self, tmp_path): + """Mixed orders + no acquisitions: each epoch lands in exactly one class.""" + fine = _write_store(tmp_path / "fine", {SHARD: _instants(0, 5)}) + coarse = self._coarse_store(tmp_path) + ref = reference_epochs([fine, coarse]) + n_coarse = int((ref.epoch_orders[SHARD_KEY] == 12).sum()) + n_fine = int(ref.epoch_orders[SHARD_KEY].size - n_coarse) + assert n_coarse > 0 and n_fine > 0 + cat = _s2_catalog( + [_s2_item("S2_far", _epoch_iso(0.1), lat=-40.0, lon=-120.0)], + bbox=(-121.0, -41.0, -119.0, -39.0), + ) + sm = closest_obs_shardmap( + cat, + [fine, coarse], + grid=_grid(), + backend="mortie", + max_time_offset=np.timedelta64(3, "D"), + ) + rec = sm.metadata["closest_obs"] + assert rec["epochs_dropped_low_resolution"] == n_coarse + # The survivors of the precision bar become no-acquisition rows. + gap_rows = [d for d in rec["dropped"] if "nearest_offset_ns" in d] + assert len(gap_rows) == n_fine + assert all(d["nearest_offset_ns"] is None for d in gap_rows) + assert rec["shards_without_acquisitions"] == [SHARD] + assert rec["epochs_total"] == rec["epochs_paired"] + rec["epochs_dropped"] + + def test_a_negative_block_order_refuses_by_name(self, tmp_path): + """Fail-CLOSED on a corrupt block order, rather than through the bar. + + ``coverage_toc._decode_cover_block``'s §10.5 order check is one-sided + (ceiling only), so ``temporal_order: -1`` decodes; the half-span shift + then overflows int64 to a huge NEGATIVE, which passes any cap, on a + midpoint ``_word_midpoints`` invents (2142-04-11 for this block). The + durable lower bound belongs beside that ceiling check; this is the + read-boundary refusal on the builder's side. + """ + root = Path(tmp_path / "corrupt") + root.mkdir(parents=True, exist_ok=True) + words = quantize_words(time2toc(_instants(40.0)), 12) + block = _encode_cover_block(np.asarray(words, np.uint64), 12) + block["temporal_order"] = -1 + (root / COVER_NAME).write_text( + json.dumps( + { + "spec": COVER_SPEC, + "source": "test", + "order": 4, + "temporal_order": TEMPORAL_COVER_ORDER, + "cap": COVER_CAP, + "fields": ["h"], + "element": {"dtype": "uint64", "shape": [-1]}, + "encoding": "base64", + "shards": {SHARD: block}, + } + ) + ) + for call in ( + lambda: reference_epochs(str(root)), + lambda: closest_obs_shardmap( + _s2_catalog([_s2_item("S2_0", _epoch_iso(40.0))]), + str(root), + grid=_grid(), + backend="mortie", + max_time_offset=np.timedelta64(3, "D"), + ), + ): + with pytest.raises(ValueError, match="temporal_order -1"): + call() diff --git a/tests/test_lambda_build.py b/tests/test_lambda_build.py index 69d4518e1..4078d5473 100644 --- a/tests/test_lambda_build.py +++ b/tests/test_lambda_build.py @@ -20,7 +20,11 @@ LAMBDA_UNZIPPED_LIMIT = 250 * 1024 * 1024 # 250MB combined (layer + function) # Budget allocation — layer gets most of the space, function code should be small -FUNCTION_SIZE_BUDGET = 30 * 1024 * 1024 # 30MB for function code +# 32MB function-code budget (espg ruling 2026-08-24, PR #511 question 1): AWS's +# hard limit for direct-upload zips is 50MB, so this is an early-warning +# tripwire, not the platform cap — 30MB left ~19KB of headroom on main and any +# source addition tripped it. Mirrored in deployment/aws/build_function.sh. +FUNCTION_SIZE_BUDGET = 32 * 1024 * 1024 class TestLambdaImports: