Skip to content

High Level HEALPix API

This section contains the main functionality for creating and working with HEALPix-based dataset representations, including pyramid construction and resolution handling.

Pyramid creation

These functions build HEALPix pyramids from gridded source data and provide the main higher-level workflow for preparing multi-resolution outputs.

create_healpix_pyramid(ds, max_level=None, min_level=0, *, coarsen_mode='auto', min_valid_fraction=0.5, valid_fraction=False, **kwargs)

Create a multi-resolution HEALPix pyramid.

The input dataset is first remapped to max_level with regrid_to_healpix. For nested output ordering, lower levels are derived efficiently with coarsen_healpix. For ring ordering, each lower level is regenerated directly from the source dataset.

For dask-backed input the result is lazy and chunked along cell (cell_chunks, forwarded to regrid_to_healpix); coarser levels keep that chunk size. Every level builds on the finest one, so write the pyramid with save_pyramid, which computes all levels in one pass. Calling .compute() per level instead regrids the source once per level.

Parameters:

Name Type Description Default
ds Dataset

Source dataset.

required
max_level int | None

Finest HEALPix level.

None
min_level int

Coarsest HEALPix level to keep.

0
coarsen_mode CoarsenMode

Coarsening strategy for building lower pyramid levels. "mean" for continuous data, "mode" for categorical, "auto" to infer from the remapping method.

'auto'
min_valid_fraction float

Minimum fraction of valid children for a parent cell to be valid. Default 0.5.

0.5
valid_fraction ValidFraction

Store <name>_valid_fraction next to the selected variables on every level: the fraction of valid finest-level cells under each cell, as float32, linked via the CF ancillary_variables attribute. Cell values are means over their valid area, so global or regional means over a coarse level need these weights, e.g. ds.sst.weighted(ds.sst_valid_fraction.fillna(0)).mean("cell").

False (default) stores nothing. True gives every cell variable a fraction of its full shape (correct for masks that change over time or height). "static" stores one cell-only fraction from the first slice along all other dimensions -- small, but only correct when the mask never changes (e.g. land/sea). A list of names selects variables (full shape); a mapping such as {"sst": "static", "ice": True} sets the shape per variable. Requires nested ordering.

False
**kwargs Any

Forwarded to regrid_to_healpix.

{}

Returns:

Type Description
dict[int, Dataset]

Pyramid keyed by level.

Source code in grid_doctor/helpers.py
def create_healpix_pyramid(
    ds: xr.Dataset,
    max_level: int | None = None,
    min_level: int = 0,
    *,
    coarsen_mode: CoarsenMode = "auto",
    min_valid_fraction: float = 0.5,
    valid_fraction: ValidFraction = False,
    **kwargs: Any,
) -> dict[int, xr.Dataset]:
    """Create a multi-resolution HEALPix pyramid.

    The input dataset is first remapped to *max_level* with
    [`regrid_to_healpix`][grid_doctor.remap.regrid_to_healpix].
    For nested output ordering, lower levels are derived efficiently
    with
    [`coarsen_healpix`][grid_doctor.helpers.coarsen_healpix].
    For ring ordering, each lower level is regenerated directly from
    the source dataset.

    For dask-backed input the result is lazy and chunked along ``cell``
    (``cell_chunks``, forwarded to
    [`regrid_to_healpix`][grid_doctor.remap.regrid_to_healpix]); coarser
    levels keep that chunk size.  Every level builds on the finest one,
    so write the pyramid with
    [`save_pyramid`][grid_doctor.helpers.save_pyramid], which computes
    all levels in one pass.  Calling ``.compute()`` per level instead
    regrids the source once per level.

    Parameters
    ----------
    ds:
        Source dataset.
    max_level:
        Finest HEALPix level.
    min_level:
        Coarsest HEALPix level to keep.
    coarsen_mode:
        Coarsening strategy for building lower pyramid levels.
        ``"mean"`` for continuous data, ``"mode"`` for categorical,
        ``"auto"`` to infer from the remapping method.
    min_valid_fraction:
        Minimum fraction of valid children for a parent cell to be
        valid.  Default ``0.5``.
    valid_fraction:
        Store ``<name>_valid_fraction`` next to the selected variables on
        every level: the fraction of valid finest-level cells under each
        cell, as float32, linked via the CF ``ancillary_variables``
        attribute.  Cell values are means over their valid area, so
        global or regional means over a coarse level need these weights,
        e.g. ``ds.sst.weighted(ds.sst_valid_fraction.fillna(0)).mean("cell")``.

        ``False`` (default) stores nothing.  ``True`` gives every cell
        variable a fraction of its full shape (correct for masks that
        change over time or height).  ``"static"`` stores one ``cell``-only
        fraction from the first slice along all other dimensions -- small,
        but only correct when the mask never changes (e.g. land/sea).
        A list of names selects variables (full shape); a mapping such as
        ``{"sst": "static", "ice": True}`` sets the shape per variable.
        Requires nested ordering.
    **kwargs:
        Forwarded to
        [`regrid_to_healpix`][grid_doctor.remap.regrid_to_healpix].

    Returns
    -------
    dict[int, xarray.Dataset]
        Pyramid keyed by level.
    """
    if max_level is None:
        max_level = resolution_to_healpix_level(get_latlon_resolution(ds))

    is_nested = bool(kwargs.get("nest", True))
    if not is_nested and valid_fraction is not False:
        raise ValueError("valid_fraction requires nested ordering (nest=True).")

    pyramid: dict[int, xr.Dataset] = {}
    finest = regrid_to_healpix(ds, max_level, **kwargs)
    pyramid[max_level] = finest

    if is_nested:
        fractions = resolve_valid_fraction(valid_fraction, finest)
        pyramid[max_level] = with_finest_fractions(finest, fractions, max_level)
        pyramid.update(
            coarse_levels(
                finest,
                max_level=max_level,
                min_level=min_level,
                coarsen_mode=resolve_coarsen_mode(finest, coarsen_mode),
                min_valid_fraction=min_valid_fraction,
                mode_kernel=_coarsen_array_mode,
                valid_fraction=fractions,
            )
        )
        return pyramid

    for level in range(max_level - 1, min_level - 1, -1):
        pyramid[level] = regrid_to_healpix(ds, level, **kwargs)
    return pyramid

Resolution helpers

These utilities help translate between source grid resolution and the corresponding HEALPix level used for remapping and pyramid generation.

get_latlon_resolution(ds)

Estimate the horizontal resolution of ds in degrees.

Parameters:

Name Type Description Default
ds Dataset

Dataset on a regular lon/lat grid, a curvilinear grid, or an unstructured grid.

required

Returns:

Type Description
float

Approximate grid spacing in degrees.

Examples:

resolution = get_latlon_resolution(ds)
level = resolution_to_healpix_level(resolution)
Source code in grid_doctor/helpers.py
def get_latlon_resolution(ds: xr.Dataset) -> float:
    """Estimate the horizontal resolution of *ds* in degrees.

    Parameters
    ----------
    ds:
        Dataset on a regular lon/lat grid, a curvilinear grid, or an
        unstructured grid.

    Returns
    -------
    float
        Approximate grid spacing in degrees.

    Examples
    --------
    ```python
    resolution = get_latlon_resolution(ds)
    level = resolution_to_healpix_level(resolution)
    ```
    """
    if _is_unstructured(ds):
        cell_dim = _get_unstructured_dim(ds)
        n_cells = ds.sizes[cell_dim]
        return float(np.degrees(np.sqrt(4.0 * np.pi / float(n_cells))))

    lat, lon = _get_latlon_arrays(ds)
    if lat.ndim == 1:
        lat_res = float(np.nanmin(np.abs(np.diff(lat))))
        lon_res = float(np.nanmin(np.abs(np.diff(lon))))
    else:
        lat_res = float(np.nanmean(np.abs(np.diff(lat, axis=0))))
        lon_res = float(np.nanmean(np.abs(np.diff(lon, axis=1))))
    return min(lat_res, lon_res)

resolution_to_healpix_level(resolution_deg)

Convert an approximate source resolution to a HEALPix level.

The heuristic uses a characteristic HEALPix pixel spacing of about 58.6° / 2**level.

Parameters:

Name Type Description Default
resolution_deg float

Approximate source resolution in degrees.

required

Returns:

Type Description
int

Suggested HEALPix level.

Examples:

level = resolution_to_healpix_level(0.25)
Source code in grid_doctor/helpers.py
def resolution_to_healpix_level(resolution_deg: float) -> int:
    """Convert an approximate source resolution to a HEALPix level.

    The heuristic uses a characteristic HEALPix pixel spacing of about
    ``58.6° / 2**level``.

    Parameters
    ----------
    resolution_deg:
        Approximate source resolution in degrees.

    Returns
    -------
    int
        Suggested HEALPix level.

    Examples
    --------
    ```python
    level = resolution_to_healpix_level(0.25)
    ```
    """
    if resolution_deg <= 0.0:
        raise ValueError("resolution_deg must be positive.")
    level = int(np.round(np.log2(58.6 / resolution_deg)))
    return max(0, level)