HEALPix-to-HEALPix resampling#

resample_healpix converts NESTED HEALPix data between Grid4Earth levels and between full or partial cell domains. It is implemented as a sequence of the package’s local HealPixDown or HealPixUp operators.

HealPixResampler provides the same operation as a reusable PyTorch module. Use the class when several maps share the same input and output geometries, so the sparse operators are constructed only once.

Quick start#

from healpix_analyse import resample_healpix

out_data, returned_ids = resample_healpix(
    in_data,
    in_level=10,
    out_level=7,
    in_cell_ids=in_cell_ids,
    out_cell_ids=out_cell_ids,
    ellipsoid="sphere",
)

returned_ids always has the same order as out_cell_ids. If out_cell_ids=None, it contains every output cell in canonical NESTED order.

Full and partial domains#

The two cell-ID arguments are independent:

Argument	`None`	Array provided
`in_cell_ids`	Input is the complete sphere at `in_level`	Last input dimension matches these IDs
`out_cell_ids`	Return the complete sphere at `out_level`	Return exactly these cells in this order

in_level and out_level are always required. A full-sphere input is selected with in_cell_ids=None, not with in_level=None, because the resolution must remain known.

Input IDs may be in any order. Output IDs may also be in any order, and that order is preserved exactly. Both arrays must contain unique valid NESTED IDs.

The three resolution cases#

Same level#

When in_level == out_level, no interpolation is performed. Values are joined directly by NESTED identifier. Requested cells absent from a partial input are filled with NaN.

Resolution reduction#

When in_level > out_level, the operator constructs

level L -> Down -> level L-1 -> ... -> requested output level

Every step uses HealPixDown(mode="smooth"). For partial maps, the exact coarse IDs are derived from the available fine cells. The final result is then joined to the requested output domain.

Resolution increase#

When in_level < out_level, every requested output ID is first mapped to its input-level NESTED ancestor:

ancestor_ids = out_cell_ids // 4 ** (out_level - in_level)

Only ancestors present in the input are retained. HealPixUp is then applied successively, keeping the complete descendant tree of every retained ancestor. The final data are gathered using the exact requested IDs and order. If an ancestor is absent, all requested descendants relying on it remain NaN.

Missing data and coastlines#

Non-finite input values are treated as missing samples. At every Up or Down step, the implementation applies the same sparse operator to:

the data with missing samples replaced by zero;
a finite-data mask;
an all-valid reference mask, precomputed when the resampler is built.

The result is renormalised by the available weighted support. Consequently:

a coastal cell remains calculable from the available ocean samples;
a channel may have a different validity mask from another channel;
a target with no non-zero support is NaN;
when all inputs are finite, the result is the unchanged output of the underlying HealPixDown or HealPixUp operator.

This behaviour follows the local masked-map philosophy of healpix-analyse.

Reusable operator#

import torch
from healpix_analyse import HealPixResampler

resampler = HealPixResampler(
    in_level=11,
    out_level=8,
    in_cell_ids=ocean_ids_level_11,
    out_cell_ids=comparison_ids_level_8,
    ellipsoid="sphere",
    dtype=torch.float64,
    device="cuda" if torch.cuda.is_available() else "cpu",
)

temperature_8, ids_8 = resampler(temperature_11)
velocity_8, ids_8 = resampler(velocity_11)

The geometry, sparse matrices, ID joins and full-support normalisations are reused by both calls.

Shapes, NumPy and Torch#

The last dimension is always the HEALPix pixel dimension. Arbitrary leading dimensions are retained:

[N]       -> [N_out]
[C, N]    -> [C, N_out]
[B, C, N] -> [B, C, N_out]

NumPy input returns a NumPy array. Torch input returns a differentiable tensor on the configured device. Output values use the configured floating-point dtype, since integer arrays cannot represent NaN.

Filtering options#

The following parameters are forwarded to the underlying operators:

Parameter	Default	Meaning
`ellipsoid`	`"WGS84"`	HEALPix geometry
`radius_deg`	level-dependent	Gaussian support radius
`sigma_deg`	level-dependent	Gaussian width
`weight_norm`	`"l1"`	Down normalisation
`up_norm`	`"col_l1"`	Up normalisation

Leaving radius_deg and sigma_deg unset lets every step select a filter appropriate to its own native level. Supplying them uses the same angular filter at every step.

Function signature#

out_data, out_ids = resample_healpix(
    in_data,
    in_level=...,
    out_level=...,
    in_cell_ids=None,
    out_cell_ids=None,
    ellipsoid="WGS84",
    radius_deg=None,
    sigma_deg=None,
    weight_norm="l1",
    up_norm="col_l1",
    dtype=torch.float32,
    device=None,
)

The arguments after in_data are keyword-only to prevent accidental exchange of input/output levels or cell-ID arrays.

Performance considerations#

Requesting the complete output sphere necessarily allocates 12 * 4**out_level identifiers and output values. Full-sphere Down operations also construct full-sphere sparse operators. For repeated comparisons, build a single HealPixResampler and reuse it.

For partial upsampling, only the input-level ancestors required by the output cells and their descendant trees are propagated. This follows NESTED hierarchy semantics and avoids constructing unrelated output regions.