Skip to content

Generic hierarchical filtering & targeted reads for multi-level point datasets #43

Description

@espg

Motivation

zagg aggregates sparse point observations onto grids, but data_source is today
flat and single-level with exactly one quality_filter that is
equality-only on a 1-D dataset (flag == value, keep matches — see
processing.py::_read_group). That is too limited for real products and is
leaving large performance on the table. Four related gaps:

  1. Memory floor of spatial filtering. To assign points to grid cells we must
    read the full coordinate arrays before we can subset. Measured on a real
    ATL03 v007 granule via h5coro: lat_ph+lon_ph full read is up to 245 MB
    per strong beam
    , ~660 MB peak per granule — paid before any filter or
    aggregation, and multiplied by granules-per-worker. This is the dominant
    driver of the Lambda OOMs on the ATL03 SERC benchmark.
  2. Filtering is underpowered. ATL03 noise removal is signal_conf_ph != 0 on
    a 2-D array (per-surface-type column) — the current equality-only / 1-D
    filter cannot express it. We also need multiple independent filters, and
    filters phrased as "keep these flag values" / "drop these".
  3. Cross-level (multi-rate) filters. Some quality flags live at a coarser
    rate than the observations. Clearest case: ATL03 podppd_flag (geolocation,
    20 m segment rate) marks geolocation failures that should drop the
    corresponding photons. We need to filter photon-rate data using segment-rate
    flags.
  4. Targeted/hierarchical reads. The same coarse level that carries those
    flags also lets us bound the read: geolocation is ~30× smaller — read it
    fully, find the segments intersecting the AOI, and read only the corresponding
    heights chunks.

These are facets of one structural idea: datasets have multiple linked levels
(resolutions), and both filtering and IO should be expressible across them.
We
want one clean, generic abstraction — built for ATL now, extending to other
missions/datasets later.

Evidence: hierarchical reads work (prototype)

benchmarks/geoloc_prefilter_proto.py (branch feature/claude-geoloc-prefilter,
off current main) implements box → segment range → ph_index_beg slice → read
only the overlapping heights chunks, and verifies it against the full-read path.

Verified index translation (over all 71,607 populated segments of a beam,
exact lat and lon match): the absolute heights index of a segment's reference
photon is ph_index_beg - 1 + reference_photon_index - 1 (both 1-based;
reference_photon_index is segment-relative, not an absolute heights
index). For bounding a read, only ph_index_beg + segment_ph_cnt are needed;
they tile heights exactly (sum(segment_ph_cnt) == n_photons).

Prototype result (one granule, ~10 km AOI, all 6 beams; in-box photon counts
match the full-read path exactly):

beam in-box photons base coord MB hierarchical MB heights chunks read
gt1l 35,834 59.3 4.27 2 / 38
gt1r 145,089 241.2 6.02 2 / 151
gt2l–gt3r (no AOI data) 0 57–245 3.69 (geoloc only) 0
total 856.8 25.0 (2.9%)

The empty-beam case is the common one for small AOIs: the geolocation index alone
(~3.7 MB) proves there is nothing to read, vs a 57–245 MB full coord read.

Proposed design direction

Generalize data_source from a flat single level to named levels linked by
index maps
, with filters scoped to a level, and a spatial pre-read
driven by a chosen level. Aggregation still happens at the finest ("base") level.

Levels + linkage

data_source:
  reader: h5coro
  driver: s3
  groups: [gt1l, gt1r, gt2l, gt2r, gt3l, gt3r]   # repeated sub-structures ("beams")
  base_level: photon            # where aggregation happens (finest)
  levels:
    photon:
      path: "/{group}/heights"
      coordinates: { latitude: lat_ph, longitude: lon_ph }
      variables:  { h_ph: h_ph }
    segment:
      path: "/{group}/geolocation"
      coordinates: { latitude: reference_photon_lat, longitude: reference_photon_lon }
      variables:  { podppd_flag: podppd_flag }
      link:                              # how this coarse level maps to a finer one
        to: photon
        index_beg: ph_index_beg          # start of each parent's children in `to`
        count: segment_ph_cnt            # children per parent
        index_base: 1                    # ph_index_beg is 1-based
        reference_index: reference_photon_index   # optional, 1-based within-parent

The link is the one piece of mission-specific knowledge, declared, not
hardcoded
. It powers both cross-level filtering and targeted reads. Generic
over index base (0/1) and start+count vs explicit start/end conventions. Other
missions declare their own (or a single level with no link).

Filters: a list, each scoped to a level, with real predicates

  filters:
    - name: drop_noise
      level: photon
      dataset: signal_conf_ph
      column: 0                          # N-D column selector (ATL03 surface-type: land)
      keep: { op: ne, value: 0 }         # or { op: in, values: [2, 3, 4] }
    - name: good_geolocation
      level: segment                     # coarse flag -> expanded to photons via link
      dataset: podppd_flag
      keep: { op: eq, value: 0 }
  • Multiple filters AND together.
  • keep operators: eq, ne, in, not_in, ge, le, lt, gt (single value, multi-value
    sets, ranges). column (or a named index) handles N-D flag arrays.
  • A coarse-level filter yields a coarse mask, expanded to the base level via the
    level's link
    (repeat each parent's verdict count times from index_beg),
    then AND-ed with base-level filters.
  • Today's quality_filter: {dataset, value} becomes sugar for one base-level
    {op: eq} filter (back-compat; migrate the atl06 configs).

Targeted (hierarchical) reads

  read_plan:
    spatial_index: segment    # bound the read using this level's coordinates
    pad: 1                    # pad the matched parent range by N elements per side

Per group: read the spatial_index level's coords (small) → AOI-match parents →
contiguous padded parent range → translate to a base-level slice via link
read only the overlapping base-level chunks (coords + needed variables + any
coarse-level flag ranges for cross-level filters) → apply combined filters →
aggregate. Falls back to a full read when no level/link is declared (today's
behavior).

Why this generalizes

Nothing above names ATL fields except inside the config. The engine only knows
levels, links (index_beg/count/base), filters (level, dataset,
column, predicate), and read via a spatial-index level. Single-level datasets
omit levels/link and use today's flat form; other instruments map their own
hierarchy (swath scan-line → pixel, beam → bin, …) onto the same primitives.

Open questions (for iteration)

  1. Schema surface vs back-compat. Keep flat coordinates/variables/
    quality_filter as sugar over a synthesized single photon level, or migrate
    all configs to levels? (Lean: support both; synthesize a base level from the
    flat form.)
  2. Column selectors by name. ATL03 signal_conf_ph columns are surface types
    (land/ocean/sea-ice/land-ice/inland-water). Allow column: land via a
    per-dataset enum map, or integer-only?
  3. Filter combination across levels. AND-only, or named groups / OR?
    (Lean: AND-only first.)
  4. Multi-run AOIs. A track can clip an AOI twice → non-contiguous parent runs.
    v1: min..max (conservative, reads the gap); later: split into runs.
  5. Where the read planner lives. Extend _read_group, or a new read_plan
    module the grid/runner call? (Affects the Lambda path; coordinate with the
    arrow refactor that just landed.)
  6. Decode the spatial-index coords once per granule and reuse across beams
    where shared (ATL: per-beam; other missions may share).

Phasing / dependencies

  • The grouping + arrow reducer (Refactor the per-cell aggregation handoff: sort/hash grouping + Arrow path (additive, benchmarked) #30 phases 1/2/2b) is now on main. First
    measure
    whether arrow alone gets the ATL03 SERC run under the 2 GB Lambda
    threshold — that calibrates urgency. The read floor is independent of arrow, so
    hierarchical reads help regardless.
  • Phase A: multiple filters + predicate operators + N-D column selector
    (base-level only). Unblocks signal_conf_ph != 0. Aggregation-time only; no
    read-path restructure.
  • Phase B: levels + link abstraction + cross-level (coarse→fine) filter
    expansion (podppd_flag).
  • Phase C: targeted hierarchical reads driven by the spatial-index level (the
    prototype), folded into the read path on top of B.
  • Each phase is independently testable; B and C share the link machinery.

Refs: prototype benchmarks/geoloc_prefilter_proto.py; memory-floor +
index-translation findings above; builds on #30 (grouping/arrow).

Metadata

Metadata

Assignees

No one assigned

    Labels

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions