Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
408 changes: 408 additions & 0 deletions benchmarks/region_timing.py

Large diffs are not rendered by default.

31 changes: 24 additions & 7 deletions src/zagg/config.py
Original file line number Diff line number Diff line change
Expand Up @@ -455,14 +455,31 @@ def _normalize_filter(f: dict) -> dict:


def _validate_filters(data_source: dict) -> None:
"""Validate the ``filters`` list (and that a flat ``quality_filter`` is sane).

Raises ``ValueError`` on: unknown op, missing ``dataset``, ``in``/``not_in``
without a list ``values``, scalar ops without ``value``, non-int ``column``,
a non-base-level ``expression`` filter, or wrong ``value`` type. ``column`` is
required for the N-D flag case but cannot be checked against array rank here
(no data); rank checks happen at read time.
"""Validate the ``filters`` list and the flat ``quality_filter`` sugar.

For the structured ``filters:`` list, raises ``ValueError`` on: unknown op,
missing ``dataset``, ``in``/``not_in`` without a list ``values``, scalar ops
without ``value``, non-int ``column``, a non-base-level ``expression`` filter,
or wrong ``value`` type. ``column`` is required for the N-D flag case but
cannot be checked against array rank here (no data); rank checks happen at
read time.

For the flat ``quality_filter`` sugar, only ``dataset`` and ``value`` are
honored (``filters_from_data_source`` synthesizes ``op: eq, column: null``).
Reject any extra keys at load time so a user-typoed ``op: gt`` or stray
``column:`` is not silently dropped on the floor — the structured ``filters:``
list is the right form when those knobs are wanted.
"""
qf = data_source.get("quality_filter")
if qf is not None:
allowed = {"dataset", "value"}
unknown = set(qf) - allowed
if unknown:
raise ValueError(
f"data_source.quality_filter only honors {sorted(allowed)} "
f"(got extra keys {sorted(unknown)}); use the structured "
"'filters:' list to set 'op', 'column', 'keep', etc."
)
filters = data_source.get("filters")
if filters is None:
return
Expand Down
95 changes: 95 additions & 0 deletions src/zagg/configs/atl03.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,95 @@
# ATL03 photon-height aggregation template.
#
# Scalar-only stats (non-scalar outputs are deferred to #29): per-cell density
# (photon count), min/max elevation, mean and median height, and variance.
#
# Confidence filter drops only TEP (transmit-echo path) photons and keeps
# everything else: signal_conf_ph != -2 (column 0, the land surface type), per
# the ATL03 v3 data dictionary
# (https://nsidc.org/sites/default/files/icesat2_atl03_data_dict_v003_0.pdf).
#
# signal_conf_ph is 2-D (n_photons x 5 surface types) and TEP is flagged -2
# uniformly across every surface type, so the land column (0) is operationally
# equivalent to any other column for the TEP drop. The structured filter list
# replaces the older flat ``quality_filter`` form to access the ``column`` and
# ``op: ne`` knobs (issue #43, Phase A). Grid is rectilinear (HEALPix order-19,
# the ~10 m match, waits on mortie #35).
#
# Multi-level form (issue #43 Phase B + C) is declared so the read path can
# bound base-rate IO via the coarse ``segments`` rep-point coords + the
# ``ph_index_beg``/``segment_ph_cnt`` link. Without this, ``lat_ph``/``lon_ph``
# would be full-read (~245 MB per beam) before any spatial filter, which is the
# dominant driver of the Lambda OOMs on ATL03 region runs (#43 motivation).
data_source:
reader: h5coro
driver: s3
groups: [gt1l, gt1r, gt2l, gt2r, gt3l, gt3r]
coordinates:
latitude: "/{group}/heights/lat_ph"
longitude: "/{group}/heights/lon_ph"
variables:
h_ph: "/{group}/heights/h_ph"
filters:
- dataset: "/{group}/heights/signal_conf_ph"
column: 0 # land surface type; TEP is uniform across columns
op: ne
value: -2
base_level: photons
# ``coordinates`` uses bare names joined to the level's ``path``; link
# ``index_beg`` / ``count`` use absolute path templates so they remain
# readable when the link arrays live outside the level's path (the read
# path resolves both forms -- see ``_level_coord_paths`` in processing.py).
# ``variables`` at the level scope is documentation only; the read path
# consults the top-level ``data_source.variables`` + filter datasets.
levels:
photons:
path: "/{group}/heights"
coordinates: {latitude: lat_ph, longitude: lon_ph}
link: null
segments:
path: "/{group}/geolocation"
coordinates: {latitude: reference_photon_lat, longitude: reference_photon_lon}
link:
to: photons
index_beg: "/{group}/geolocation/ph_index_beg"
count: "/{group}/geolocation/segment_ph_cnt"

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

Asymmetric path conventions: coordinates use bare names (joined to path), link fields use absolute /{group}/... templates.

The level's coordinates here are bare names — reference_photon_lat / reference_photon_lon — which _level_coord_paths (src/zagg/processing.py:837) joins to the level's path template (/{group}/geolocation) to produce /{group}/geolocation/reference_photon_lat. That join behavior is explicit:

lat_path = lat_name if lat_name.startswith("/") else f"{base}/{lat_name}"

But the link fields one line below — index_beg: "/{group}/geolocation/ph_index_beg", count: "/{group}/geolocation/segment_ph_cnt" — are absolute templates. _planned_read_group resolves them with link["index_beg"].format(group=group) (processing.py:892); there is no path-join. Both shapes work today (absolute paths bypass the join short-circuit), but the inconsistency will trip the next person who copies issue #43's schema example literally — that example uses bare names for both coordinates and link fields:

link:
  to: photon
  index_beg: ph_index_beg          # start of each parent's children in `to`
  count: segment_ph_cnt            # children per parent

Two options, both small:

  • Pick relative everywhere on both ATL03 YAMLs: index_beg: ph_index_beg, count: segment_ph_cnt. Requires extending _planned_read_group (and the cross-level filter at processing.py:1007–1009) to apply the same join rule the coordinate resolver uses. That's where this is heading per Generic hierarchical filtering & targeted reads for multi-level point datasets #43, so it's the future-proof choice.
  • Pick absolute everywhere on both ATL03 YAMLs: rewrite the coords as latitude: "/{group}/geolocation/reference_photon_lat". Smaller change, but it freezes the asymmetry in for ATL03 specifically.

If you want neither now, please at minimum add a YAML comment block at the top of levels: noting which convention the link fields follow — otherwise the next mission's template will be written against the issue #43 spec and fail at link["index_beg"].format(group=group) when no / is there to absolute-prefix.


Generated by Claude Code

index_base: 1 # ATL03 ph_index_beg is 1-based per the v3 dict
read_plan:
spatial_index: segments
pad: 1 # one segment of padding on each side per #43

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🤖 from Claude (review)

pad: 1 choice is undocumented and not exercised by any test in this PR.

The YAML comment says # one segment of padding on each side per #43. #43's issue body in fact defines pad ("pad the matched parent range by N elements per side") without committing to a default; the prototype used 1 but doesn't justify it as the operational default for the grids zagg actually ships. The phase 5 planned-read tests use pad: 0 (tests/test_processing.py:1748) for tractability; TestPlanReadIndexBase uses pad: 0 for the index_base=1 case (tests/test_read_plan.py:131); the only pad: 1 exercise is implicit in TestPlannedReadGroup::test_planned_path_bounds_io_to_matched_segments where the wide segment spacing makes pad behavior irrelevant.

Concretely:

  • ATL03 segment spacing is ~20 m along-track. pad: 1 = one segment = ~20 m of padding on each side of a matched run.
  • The rectilinear grid in both YAMLs is resolution: 0.0001 (~10 m at the equator). A shard footprint of 10 m × 10 m sits inside a single 20 m segment — pad: 1 means up to ~20 m of always-discarded reads per shard, on top of the actual span. Not wasteful in absolute terms but worth knowing.
  • At high latitude (~88°S — the --hard AOI in the benchmark), a rectilinear cell of 0.0001° in longitude is meters across, but the segment footprint stretches differently. The 1-segment pad might genuinely under-pad there if the AOI clips a segment boundary at an angle.

This isn't blocking, but two requests:

  1. Add a YAML comment that justifies pad: 1 from the grid resolution, not just "per Generic hierarchical filtering & targeted reads for multi-level point datasets #43" — e.g. "1 segment ≈ 20 m, matches the 0.0001° rectilinear resolution to within an order of magnitude; tune if HEALPix order-19 lands."
  2. Either add an offline test that exercises pad: 1 at a segment boundary on the shipped template (would fit naturally into the integration test I'm flagging on tests/test_config.py:164), or at least add a unit test in tests/test_read_plan.py that confirms pad: 1 extends a single-segment match into a 3-segment slice with index_base=1 arithmetic. The current TestPlanReadIndexBase test uses pad: 0.

Generated by Claude Code


aggregation:
variables:
count:
function: len
source: h_ph
dtype: int32
fill_value: 0
h_min:
function: min
source: h_ph
dtype: float32
h_max:
function: max
source: h_ph
dtype: float32
h_mean:
function: mean
source: h_ph
dtype: float32
h_median:
function: median
source: h_ph
dtype: float32
h_variance:
function: var
source: h_ph
dtype: float32

output:
grid:
type: rectilinear
crs: EPSG:4326
resolution: 0.0001 # degrees (~10 m at the equator)
bounds: [-180, -90, 180, 90]
chunk_shape: [256, 256]
83 changes: 83 additions & 0 deletions src/zagg/configs/atl03_waveform_counts.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
# ATL03 per-cell waveform-counts template.
#
# Aggregates photon heights (h_ph) into a 128-bin fixed-width histogram at 2 m
# vertical spacing, centred on the per-cell median of h_ph itself. Each bin
# spans 2 m, giving a 256 m window [median(h_ph) - 128 m, median(h_ph) + 128 m).
# Out-of-range photons are dropped (not counted). The companion scalar
# bin_start records median(h_ph) so the absolute elevation of bin 0 can be
# recovered downstream.
#
# Note (issue #30 thread): an earlier draft of this template tried to centre
# on np.median(dem_h), but dem_h lives in /{group}/geophys_corr/dem_h at the
# segment level (one value per ~100 photons), not the heights group; centring
# on the cell's own median(h_ph) is more natural for a photon-height density
# template and avoids the cross-group / resolution-mismatch read.
#
# Caveat: median(h_ph) is a *per-cell* anchor, not a stable absolute reference.
# In a bimodal cell (e.g. ground + canopy returns) the median can land between
# the two clusters, splitting both into wings of the histogram; that is a real
# tradeoff of centring on the photon distribution itself. bin_start records the
# centring value so the bin-0 absolute elevation = (bin_start - 128 m) is still
# recoverable, but comparing waveforms across cells in the same absolute frame
# requires shifting each cell's bins by its own bin_start.
#
# Confidence filter is the same as atl03.yaml: drop only TEP photons by the
# structured ``filters:`` list (signal_conf_ph[:, 0] != -2; TEP is uniform
# across surface types). Multi-level form + ``read_plan`` also identical to
# atl03.yaml so the planned-IO benefits (#43 Phase C) apply to the waveform
# template equally. Grid is rectilinear, matching atl03.yaml.
data_source:
reader: h5coro
driver: s3
groups: [gt1l, gt1r, gt2l, gt2r, gt3l, gt3r]
coordinates:
latitude: "/{group}/heights/lat_ph"
longitude: "/{group}/heights/lon_ph"
variables:
h_ph: "/{group}/heights/h_ph"
filters:
- dataset: "/{group}/heights/signal_conf_ph"
column: 0 # land surface type; TEP is uniform across columns
op: ne
value: -2
base_level: photons
# See atl03.yaml for the rationale behind the bare-name coordinates +
# absolute link path templates + documentation-only level ``variables``.
levels:
photons:
path: "/{group}/heights"
coordinates: {latitude: lat_ph, longitude: lon_ph}
link: null
segments:
path: "/{group}/geolocation"
coordinates: {latitude: reference_photon_lat, longitude: reference_photon_lon}
link:
to: photons
index_beg: "/{group}/geolocation/ph_index_beg"
count: "/{group}/geolocation/segment_ph_cnt"
index_base: 1 # ATL03 ph_index_beg is 1-based per the v3 dict
read_plan:
spatial_index: segments
pad: 1 # one segment of padding on each side per #43

aggregation:
variables:
waveform_counts:
kind: vector
trailing_shape: 128
expression: "np.histogram(h_ph - np.median(h_ph), 128, (-128.0, 128.0))[0].astype('uint32')"
source: h_ph
dtype: uint32
fill_value: 0
bin_start:
function: median
source: h_ph
dtype: float32

output:
grid:
type: rectilinear
crs: EPSG:4326
resolution: 0.0001 # degrees (~10 m at the equator)
bounds: [-180, -90, 180, 90]
chunk_shape: [256, 256]
Loading
Loading