Processing one parent cell at a time (subset_for_parent_cell)#

KNeighborsResampler.__init__ computes a healpix_geo indexing pass and a KNN neighbourhood search over every input sample up front, and comp_matrix() materializes sparse (N, K)/(K, N) operators sized to that same N. For a global dataset, that doesn’t fit in memory (and isn’t necessary) if you only want results for one region at a time.

healpix_resample.subset_for_parent_cell restricts both sides of the problem to one coarse HEALPix “parent cell”: the output cells (an exact, cheap lookup) and the input samples actually relevant to them (the expensive part, and the actual point of this helper).

Why out_cell_ids alone doesn’t solve this#

Every KNN-mode resampler already accepts an out_cell_ids= kwarg to restrict its output. But that restriction is applied only after the expensive part – indexing and the neighbourhood search still run over the full, unfiltered input array. subset_for_parent_cell adds the missing piece: filtering the input down to a local subset before a resampler is ever constructed, so N in the constructor call is the local sample count, not the global one.

Usage#

from healpix_resample import PSFResampler, subset_for_parent_cell

level_parent = 6   # coarse: nside = 64
level = 20          # fine target resolution

sample_idx, out_ids = subset_for_parent_cell(
    lon, lat,
    parent_cell_id=parent_id,
    level_parent=level_parent,
    level=level,
)

op = PSFResampler(lon_deg=lon[sample_idx], lat_deg=lat[sample_idx], level=level, out_cell_ids=out_ids)
result = op.resample(val[..., sample_idx])

Loop this over every parent cell at level_parent to process a full global dataset in bounded-memory chunks. Reassembling results from multiple parent cells into one global map is not something this function does for you – it only makes one parent cell’s computation tractable.

What it returns#

Deliberately not a filtered copy of val – see “Why this returns an index” below.

  • sample_idx: an integer index into the sample axis – samples whose own level_parent cell is parent_cell_id itself, plus a margin_rings-neighbour buffer around it (see below). Index lon/lat and any value array yourself: lon[sample_idx], val[..., sample_idx] (the last-axis indexing works whether val is (N,) or batched (B, N)).

  • out_cell_ids: the level-resolution cells contained in parent_cell_id (healpix_geo.*.zoom_to), ready to pass straight into a KNN-mode resampler’s out_cell_ids=.

Why this returns an index rather than filtered arrays#

subset_for_parent_cell doesn’t take a val and doesn’t hand back lon_sub/lat_sub/val_sub. A single (lon, lat) grid is very commonly shared by several different value arrays – multiple variables, a time series over the same station network – and which samples are relevant to parent_cell_id only depends on lon/lat, never on val. Returning a plain index lets you compute it once per grid and reuse it across every co-located array:

sample_idx, out_ids = subset_for_parent_cell(lon, lat, parent_cell_id=parent_id, level_parent=6, level=20)

lon_sub, lat_sub = lon[sample_idx], lat[sample_idx]
temperature_sub = temperature[..., sample_idx]
pressure_sub = pressure[..., sample_idx]

instead of recomputing the same geometric membership test once per variable.

The margin is a correctness guarantee, not a performance knob#

A sample can sit just outside parent_cell_id’s boundary, in a sibling level_parent cell, while still being close enough in the fine-level kernel sense (sigma_m/threshold) to legitimately contribute to a fine cell near the parent’s edge. Filtering input strictly to “same parent cell id” would silently starve edge cells of legitimate neighbours – a subtle bug that only shows up as slightly-degraded results near parent-cell boundaries, not as an error.

margin_rings (default 1) is measured in fine-level HEALPix rings, not level_parent rings: a sample is kept if its own level cell is within margin_rings rings of any fine cell inside parent_cell_id (i.e. of out_cell_ids, expanded via healpix_geo.*.kth_neighbourhood at level). This matters because a fine cell is typically vastly smaller than a parent cell – buffering by whole neighbouring parent-level cells instead (one ring already means “pull in every neighbouring parent cell in full”) would keep a hugely disproportionate number of irrelevant samples relative to the resampler’s actual kernel reach at level. This is a correctness guarantee, not just a performance knob: the margin must be wide enough that the resampler’s actual kernel reach (sigma_m and the effective radius implied by threshold) can never extend past it, but sized at the fine scale the kernel actually operates at. margin_rings=1 is very likely generous relative to a fine cell’s own width, but this depends on your sigma_m/threshold choice, not on this function’s defaults – if you shrink sigma_m or loosen threshold enough that the kernel reaches past a handful of fine cells, increase margin_rings accordingly. As a rule of thumb, compare a fine cell’s angular width (sqrt(4*pi/(12*4**level)) radians) against the kernel’s actual reach for your resampler, and size the margin so the kernel can never reach past it.

A separate, smaller caveat: out_cell_ids gap-filling fallbacks are local, not global#

Even with a correctly-sized margin, don’t expect a per-parent-cell reassembly to match a global run to floating-point exactness on every cell. NearestResampler (and, in its own out_cell_ids fallback path, PSFResampler) can end up with cells inside out_cell_ids that the KNN ring search doesn’t reach on its own; both patch these in via a fallback that searches for the nearest/best sample within whatever sample set that resampler instance was constructed from. Constructed from a margin-filtered local subset, that fallback can pick a different (still reasonable, but not necessarily identical) answer than the same fallback would running against the full global dataset. This is a distinct, smaller effect from the margin discussed above – it doesn’t grow if you widen margin_rings further once the margin already comfortably covers the kernel’s real reach – so treat a handful of cells disagreeing at parent-cell boundaries as expected background noise from this fallback interaction, not evidence the margin itself is wrong. The tutorial (docs/tutorials/5parent_cell_subsetting.md) reports the mismatch rate directly so you can see this in practice.

Resamplers that use group_by=True#

ConservativeResampler, GroupByResampler, and CellPointResampler derive their cells purely from the samples actually present in the (already margin-filtered) input – they have no KNN neighbourhood search and no out_cell_ids intersection logic to hook into:

  • ConservativeResampler raises NotImplementedError if you pass out_cell_ids at all.

  • GroupByResampler/CellPointResampler silently store but never use it (no error, no effect).

For these classes, use only sample_idx from subset_for_parent_cell to slice lon/lat/val yourself, and do not pass its out_cell_ids to their constructors – the set of cells they produce is simply whatever the filtered input actually hits.

from healpix_resample import ConservativeResampler, subset_for_parent_cell

sample_idx, _out_ids_unused = subset_for_parent_cell(
    lon, lat, parent_cell_id=parent_id, level_parent=6, level=20,
)
op = ConservativeResampler(
    lon_deg=lon[sample_idx], lat_deg=lat[sample_idx], level=20, area=area[sample_idx]
)
result = op.resample(val[..., sample_idx])