Regridding

This section covers the core remapping functionality of grid_doctor. The functions below are the main entry points for generating reusable weight files and transforming data onto HEALPix grids.

Weight file generation

These functions prepare reusable remapping weights so that later transformations can be applied efficiently without recomputing the full mapping each time.

grid_doctor.utils.cached_weights(ds, level=None, *, method='conservative', nest=True, source_units='auto', cache_path=None, **kwargs)

Compute or load a cached HEALPix weight file.

Parameters:

Name	Type	Description	Default
ds	Dataset	Source dataset whose grid geometry defines the weights.	required
level	int | None	HEALPix level.	None
method	RemapMethod	Weight-generation method. Supported values are "nearest" and "conservative".	'conservative'
nest	bool	Use nested HEALPix ordering when True.	True
source_units	SourceUnits	Unit convention of the source coordinates.	'auto'
cache_path	str | Path | None	Cache directory or explicit file name. When omitted, the default package cache directory is used.	None
**kwargs	Any	Any additional keyword arguments for compute_healpix_weights	{}

Returns:

Type	Description
Path	Path to the cached NetCDF weight file.

Examples:

weight_file = cached_weights(ds, level=8, method="conservative")

Source code in grid_doctor/utils.py

def cached_weights(
    ds: xr.Dataset,
    level: int | None = None,
    *,
    method: RemapMethod = "conservative",
    nest: bool = True,
    source_units: SourceUnits = "auto",
    cache_path: str | Path | None = None,
    **kwargs: Any,
) -> Path:
    """Compute or load a cached HEALPix weight file.

    Parameters
    ----------
    ds:
        Source dataset whose grid geometry defines the weights.
    level:
        HEALPix level.
    method:
        Weight-generation method. Supported values are `"nearest"` and
        `"conservative"`.
    nest:
        Use nested HEALPix ordering when `True`.
    source_units:
        Unit convention of the source coordinates.
    cache_path:
        Cache directory or explicit file name. When omitted,
        the default package cache directory is used.
    **kwargs:
        Any additional keyword arguments for
        [`compute_healpix_weights`][grid_doctor.remap.compute_healpix_weights]

    Returns
    -------
    pathlib.Path
        Path to the cached NetCDF weight file.

    Examples
    --------
    ```python
    weight_file = cached_weights(ds, level=8, method="conservative")
    ```
    """
    from .helpers import get_latlon_resolution, resolution_to_healpix_level
    from .remap import compute_healpix_weights

    digest = hashlib.sha256()
    for candidate in (
        "clon_vertices",
        "clat_vertices",
        "lon_vertices",
        "lat_vertices",
        "clon",
        "clat",
        "lon",
        "lat",
        "longitude",
        "latitude",
        "rlon",
        "rlat",
        "X",
        "Y",
    ):
        if candidate in ds:
            digest.update(
                np.ascontiguousarray(np.asarray(ds[candidate].values)).tobytes()
            )
    digest.update(
        f"level={level};method={method};nest={nest};units={source_units}".encode()
    )
    key = digest.hexdigest()[:16]

    if cache_path is None:
        weight_file = cache_dir() / f"weights_{key}.nc"
    else:
        path = Path(cache_path)
        weight_file = path / f"weights_{key}.nc" if path.is_dir() else path
    if weight_file.exists():
        logger.info("Using cached weight file %s", weight_file)
        return weight_file

    logger.info("Generating HEALPix weight file %s", weight_file)
    if level is None:
        level = resolution_to_healpix_level(get_latlon_resolution(ds))

    return compute_healpix_weights(
        ds,
        level,
        method=method,
        nest=nest,
        source_units=source_units,
        weights_path=weight_file,
        **kwargs,
    )

compute_healpix_weights(ds, level, *, method='nearest', nest=True, source_units='auto', weights_path=None, grid=None, source_kind='auto', ignore_unmapped=None, large_file=True, prefer_offline=None, nproc=1, esmf_regrid_weightgen='ESMF_RegridWeightGen', keep_intermediates=False, workdir=None, spectral_transform_command=None)

Generate a reusable NetCDF weight file for HEALPix remapping.

Parameters:

Name	Type	Description	Default
ds	Dataset | str | Path	Source dataset or a source path.	required
level	int	HEALPix refinement level. nside = 2**level.	required
method	RemapMethod	"nearest" or "conservative".	'nearest'
nest	bool	Use nested HEALPix ordering when True.	True
source_units	SourceUnits	Unit convention of the source coordinates.	'auto'
weights_path	str | Path | None	Output NetCDF file. A temporary file is created when omitted.	None
grid	Dataset | None	Optional external geometry dataset.	None
source_kind	SourceKind	Explicit source representation or "auto".	'auto'
ignore_unmapped	bool | None	Ignore unmapped destination cells.	None
large_file	bool	Forwarded to the in-memory ESMPy workflow.	True
prefer_offline	bool | None	Force or disable the offline ESMF path.	None
nproc	int	Number of MPI ranks for the offline path.	1
esmf_regrid_weightgen	str	Name or path of the offline ESMF executable.	'ESMF_RegridWeightGen'
keep_intermediates	bool	Keep intermediate mesh files.	False
workdir	str | Path | None	Working directory for offline intermediate files.	None
spectral_transform_command	list[str] | tuple[str, ...] | None	External command for source_kind="spectral".	None

Returns:

Type	Description
Path	Path to the generated NetCDF weight file.

See Also

apply_weight_file: Apply a previously generated weight file.

Source code in grid_doctor/remap.py

def compute_healpix_weights(
    ds: xr.Dataset | str | Path,
    level: int,
    *,
    method: RemapMethod = "nearest",
    nest: bool = True,
    source_units: SourceUnits = "auto",
    weights_path: str | Path | None = None,
    grid: xr.Dataset | None = None,
    source_kind: SourceKind = "auto",
    ignore_unmapped: bool | None = None,
    large_file: bool = True,
    prefer_offline: bool | None = None,
    nproc: int = 1,
    esmf_regrid_weightgen: str = "ESMF_RegridWeightGen",
    keep_intermediates: bool = False,
    workdir: str | Path | None = None,
    spectral_transform_command: list[str] | tuple[str, ...] | None = None,
) -> Path:
    """Generate a reusable NetCDF weight file for HEALPix remapping.

    Parameters
    ----------
    ds:
        Source dataset or a source path.
    level:
        HEALPix refinement level.  ``nside = 2**level``.
    method:
        ``"nearest"`` or ``"conservative"``.
    nest:
        Use nested HEALPix ordering when `True`.
    source_units:
        Unit convention of the source coordinates.
    weights_path:
        Output NetCDF file.  A temporary file is created when omitted.
    grid:
        Optional external geometry dataset.
    source_kind:
        Explicit source representation or ``"auto"``.
    ignore_unmapped:
        Ignore unmapped destination cells.
    large_file:
        Forwarded to the in-memory ESMPy workflow.
    prefer_offline:
        Force or disable the offline ESMF path.
    nproc:
        Number of MPI ranks for the offline path.
    esmf_regrid_weightgen:
        Name or path of the offline ESMF executable.
    keep_intermediates:
        Keep intermediate mesh files.
    workdir:
        Working directory for offline intermediate files.
    spectral_transform_command:
        External command for ``source_kind="spectral"``.

    Returns
    -------
    pathlib.Path
        Path to the generated NetCDF weight file.

    See Also
    --------
    [`apply_weight_file`][grid_doctor.remap.apply_weight_file]:
        Apply a previously generated weight file.
    """
    offline = OfflineWeightConfig(
        enabled=bool(prefer_offline) if prefer_offline is not None else False,
        nproc=nproc,
        esmf_regrid_weightgen=esmf_regrid_weightgen,
        keep_intermediates=keep_intermediates,
        workdir=Path(workdir) if workdir is not None else None,
    )
    return compute_healpix_weights_backend(
        ds,
        level,
        method=method,
        nest=nest,
        source_units=source_units,
        weights_path=weights_path,
        grid=grid,
        source_kind=source_kind,
        ignore_unmapped=ignore_unmapped,
        large_file=large_file,
        offline=offline,
        spectral_transform_command=spectral_transform_command,
    )

Applying remapping

The following functions use either direct remapping logic or precomputed weights to move data onto the target HEALPix representation.

regrid_to_healpix(ds, level, *, nest=True, method='conservative', source_units='auto', weights_path=None, missing_policy='renormalize', backend='auto', grid=None, source_kind='auto', ignore_unmapped=None, large_file=True, prefer_offline=None, nproc=1, esmf_regrid_weightgen='ESMF_RegridWeightGen', keep_intermediates=False, workdir=None, spectral_transform_command=None, cell_chunks='auto')

Regrid ds to a HEALPix target grid.

Parameters:

Name	Type	Description	Default
ds	Dataset	Source dataset on a regular, curvilinear, or unstructured grid.	required
level	int	HEALPix refinement level.	required
nest	bool	Use nested HEALPix ordering when True.	True
method	RemapMethod	"nearest" or "conservative".	'conservative'
source_units	SourceUnits	Coordinate units for the source grid.	'auto'
weights_path	str | Path | None	Optional existing or target weight file.	None
missing_policy	MissingPolicy	How NaN values in the source grid are treated. "renormalize" (default) excludes NaN source cells and rescales the remaining weights so that a target cell is valid whenever at least one source cell contributes — appropriate for most continuous fields. "propagate" sets a target cell to NaN if any contributing source cell is NaN — useful for strict budget-closure checks but very aggressive along coastlines and data edges. See apply_weight_file for a detailed description.	'renormalize'
backend	ApplyBackend	Application backend ("auto", "scipy", "numba").	'auto'
grid	Dataset | None	Optional external geometry dataset.	None
source_kind	SourceKind	Explicit source representation or "auto".	'auto'
ignore_unmapped	bool | None	Ignore unmapped destination cells.	None
large_file	bool	Forwarded to the in-memory ESMPy workflow.	True
prefer_offline	bool | None	Force or disable the offline ESMF path.	None
nproc	int	Number of MPI ranks for the offline path.	1
esmf_regrid_weightgen	str	Name or path of the offline ESMF executable.	'ESMF_RegridWeightGen'
keep_intermediates	bool	Keep intermediate mesh files.	False
workdir	str | Path | None	Working directory for offline intermediate files.	None
spectral_transform_command	list[str] | tuple[str, ...] | None	External command for source_kind="spectral".	None
cell_chunks	int | Literal['auto'] | None	Output chunk size along cell for dask-backed input, see apply_weight_file. Coarser pyramid levels keep the same chunk size.	'auto'

Returns:

Type	Description
Dataset	Regridded dataset on the HEALPix grid.

Source code in grid_doctor/remap.py

def regrid_to_healpix(
    ds: xr.Dataset,
    level: int,
    *,
    nest: bool = True,
    method: RemapMethod = "conservative",
    source_units: SourceUnits = "auto",
    weights_path: str | Path | None = None,
    missing_policy: MissingPolicy = "renormalize",
    backend: ApplyBackend = "auto",
    grid: xr.Dataset | None = None,
    source_kind: SourceKind = "auto",
    ignore_unmapped: bool | None = None,
    large_file: bool = True,
    prefer_offline: bool | None = None,
    nproc: int = 1,
    esmf_regrid_weightgen: str = "ESMF_RegridWeightGen",
    keep_intermediates: bool = False,
    workdir: str | Path | None = None,
    spectral_transform_command: list[str] | tuple[str, ...] | None = None,
    cell_chunks: int | Literal["auto"] | None = "auto",
) -> xr.Dataset:
    """Regrid *ds* to a HEALPix target grid.

    Parameters
    ----------
    ds:
        Source dataset on a regular, curvilinear, or unstructured grid.
    level:
        HEALPix refinement level.
    nest:
        Use nested HEALPix ordering when ``True``.
    method:
        ``"nearest"`` or ``"conservative"``.
    source_units:
        Coordinate units for the source grid.
    weights_path:
        Optional existing or target weight file.
    missing_policy:
        How NaN values in the source grid are treated.
        ``"renormalize"`` (default) excludes NaN source cells and
        rescales the remaining weights so that a target cell is valid
        whenever at least one source cell contributes — appropriate
        for most continuous fields.  ``"propagate"`` sets a target
        cell to NaN if any contributing source cell is NaN — useful
        for strict budget-closure checks but very aggressive along
        coastlines and data edges.  See
        [`apply_weight_file`][grid_doctor.remap.apply_weight_file]
        for a detailed description.
    backend:
        Application backend (``"auto"``, ``"scipy"``, ``"numba"``).
    grid:
        Optional external geometry dataset.
    source_kind:
        Explicit source representation or ``"auto"``.
    ignore_unmapped:
        Ignore unmapped destination cells.
    large_file:
        Forwarded to the in-memory ESMPy workflow.
    prefer_offline:
        Force or disable the offline ESMF path.
    nproc:
        Number of MPI ranks for the offline path.
    esmf_regrid_weightgen:
        Name or path of the offline ESMF executable.
    keep_intermediates:
        Keep intermediate mesh files.
    workdir:
        Working directory for offline intermediate files.
    spectral_transform_command:
        External command for ``source_kind="spectral"``.
    cell_chunks:
        Output chunk size along ``cell`` for dask-backed input, see
        [`apply_weight_file`][grid_doctor.remap.apply_weight_file].
        Coarser pyramid levels keep the same chunk size.

    Returns
    -------
    xarray.Dataset
        Regridded dataset on the HEALPix grid.
    """
    weight_file = Path(weights_path) if weights_path is not None else None
    if weight_file is None or not weight_file.exists():
        weight_file = compute_healpix_weights(
            ds,
            level,
            method=method,
            nest=nest,
            source_units=source_units,
            weights_path=weight_file,
            grid=grid,
            source_kind=source_kind,
            prefer_offline=prefer_offline,
            nproc=nproc,
            esmf_regrid_weightgen=esmf_regrid_weightgen,
            keep_intermediates=keep_intermediates,
            workdir=workdir,
            spectral_transform_command=spectral_transform_command,
            ignore_unmapped=ignore_unmapped,
            large_file=large_file,
        )
    return apply_weight_file(
        ds,
        weight_file,
        missing_policy=missing_policy,
        backend=backend,
        cell_chunks=cell_chunks,
    )

apply_weight_file(ds, weights_path, *, missing_policy='renormalize', grid=None, source_dims=None, source_units='auto', backend='auto', cell_chunks='auto')

Apply a previously generated ESMF weight file to ds.

Uses the batched engine in grid_doctor.remap_apply which replaces per-slice vectorize=True with a single batched sparse matmul for a typical 10–50× speedup.

Parameters:

Name	Type	Description	Default
ds	Dataset	Source dataset to remap.	required
weights_path	str | Path	Path to the NetCDF weight file.	required
missing_policy	MissingPolicy	Strategy for handling NaN (missing) values in the source grid during weight application. "renormalize" (default) Source cells that are NaN are excluded from the weighted sum and the weights of the remaining valid source cells are rescaled so they sum to 1. A target cell receives a valid value as long as at least one contributing source cell is valid. This is the right choice for most climate fields: for instance, a HEALPix cell on a coastline that partly overlaps land (NaN) and partly ocean still receives a valid SST from the ocean fraction. "propagate" If any contributing source cell is NaN, the target cell is set to NaN — no partial averages. Use this when completeness of the source coverage matters more than spatial coverage of the output, for example when computing area-integrated flux budgets where a partial average would bias the integral. Note that this is very aggressive: a single NaN source cell at the edge of a target cell's footprint is enough to discard it. For conservative remapping where a HEALPix cell can overlap dozens of source cells, large portions of the output along coastlines, swath edges, and orbit gaps will be NaN.	'renormalize'
grid	Dataset | None	Optional grid dataset with source geometry.	None
source_dims	tuple[str, ...] | None	Optional explicit source dimensions.	None
source_units	SourceUnits	Unit convention for geometry-based inference.	'auto'
backend	ApplyBackend	Application backend ("auto", "scipy", "numba").	'auto'
cell_chunks	int | Literal['auto'] | None	Chunk size of the output along cell for dask-backed input. Each output chunk is computed from its own slice of weight-matrix rows, so no task holds the whole HEALPix field. "auto" uses DEFAULT_CELL_CHUNK cells (rounded to a multiple of 4), None produces a single chunk. NumPy-backed input is always computed in one piece.	'auto'

Returns:

Type	Description
Dataset	Dataset on the HEALPix target grid with a cell dimension.

Examples:

ds_hp = apply_weight_file(ds_icon, "icon_to_healpix_d8.nc")

Source code in grid_doctor/remap.py

def apply_weight_file(
    ds: xr.Dataset,
    weights_path: str | Path,
    *,
    missing_policy: MissingPolicy = "renormalize",
    grid: xr.Dataset | None = None,
    source_dims: tuple[str, ...] | None = None,
    source_units: SourceUnits = "auto",
    backend: ApplyBackend = "auto",
    cell_chunks: int | Literal["auto"] | None = "auto",
) -> xr.Dataset:
    """Apply a previously generated ESMF weight file to *ds*.

    Uses the batched engine in
    `grid_doctor.remap_apply` which replaces per-slice
    ``vectorize=True`` with a single batched sparse matmul for a
    typical 10–50× speedup.

    Parameters
    ----------
    ds:
        Source dataset to remap.
    weights_path:
        Path to the NetCDF weight file.
    missing_policy:
        Strategy for handling NaN (missing) values in the source
        grid during weight application.

        ``"renormalize"`` (default)
            Source cells that are NaN are excluded from the
            weighted sum and the weights of the remaining valid
            source cells are rescaled so they sum to 1.  A target
            cell receives a valid value as long as **at least one**
            contributing source cell is valid.  This is the right
            choice for most climate fields: for instance, a
            HEALPix cell on a coastline that partly overlaps land
            (NaN) and partly ocean still receives a valid SST from
            the ocean fraction.

        ``"propagate"``
            If **any** contributing source cell is NaN, the target
            cell is set to NaN — no partial averages.  Use this
            when completeness of the source coverage matters more
            than spatial coverage of the output, for example when
            computing area-integrated flux budgets where a partial
            average would bias the integral.  Note that this is
            very aggressive: a single NaN source cell at the edge
            of a target cell's footprint is enough to discard it.
            For conservative remapping where a HEALPix cell can
            overlap dozens of source cells, large portions of the
            output along coastlines, swath edges, and orbit gaps
            will be NaN.
    grid:
        Optional grid dataset with source geometry.
    source_dims:
        Optional explicit source dimensions.
    source_units:
        Unit convention for geometry-based inference.
    backend:
        Application backend (``"auto"``, ``"scipy"``, ``"numba"``).

    cell_chunks:
        Chunk size of the output along ``cell`` for dask-backed input.
        Each output chunk is computed from its own slice of weight-matrix
        rows, so no task holds the whole HEALPix field.  ``"auto"`` uses
        ``DEFAULT_CELL_CHUNK`` cells (rounded to a multiple of 4), ``None``
        produces a single chunk.  NumPy-backed input is always computed
        in one piece.

    Returns
    -------
    xarray.Dataset
        Dataset on the HEALPix target grid with a ``cell`` dimension.

    Examples
    --------
    ```python
    ds_hp = apply_weight_file(ds_icon, "icon_to_healpix_d8.nc")
    ```
    """
    if missing_policy not in {"renormalize", "propagate"}:
        raise ValueError("missing_policy must be 'renormalize' or 'propagate'.")

    matrix, n_target, n_source, level, order, method, stored_sd = _read_weight_file(
        weights_path
    )

    resolved_sd = _resolve_source_dims_for_weight_application(
        ds,
        n_source=n_source,
        grid=grid,
        source_dims=source_dims,
        source_units=source_units,
        stored_source_dims=stored_sd,
    )
    n_src_dims = len(resolved_sd)
    cell_chunk = _resolve_cell_chunk(cell_chunks, n_target)
    row_blocks: list[Any] | None = None

    regridded: dict[str, xr.DataArray] = {}
    for name, data in ds.data_vars.items():
        if not set(resolved_sd).issubset(map(str, data.dims)):
            regridded[str(name)] = data
            continue

        if cell_chunk is not None and data.chunks is not None:
            if row_blocks is None:
                row_blocks = _row_blocks(matrix, cell_chunk)
            regridded[str(name)] = _apply_weights_rowblocked(
                data,
                resolved_sd,
                row_blocks,
                missing_policy=missing_policy,
                backend=backend,
                parallel=False,
            )
            continue

        regridded[str(name)] = cast(
            xr.DataArray,
            xr.apply_ufunc(
                apply_weights_nd,
                data,
                input_core_dims=[list(resolved_sd)],
                output_core_dims=[["cell"]],
                exclude_dims=set(resolved_sd),
                dask="parallelized",
                kwargs={
                    "matrix": matrix,
                    "n_source_dims": n_src_dims,
                    "missing_policy": missing_policy,
                    "backend": backend,
                    # dask already runs tasks in parallel threads; nested
                    # thread teams exhaust the process thread limit.
                    "parallel": data.chunks is None,
                },
                output_dtypes=[np.float64],
                dask_gufunc_kwargs={"output_sizes": {"cell": n_target}},
            ),
        )

    result = xr.Dataset(regridded, attrs=ds.attrs.copy())
    if level >= 0:
        result = _attach_healpix_coords(
            result,
            level=level,
            nest=(order == "nested"),
            method=method,
        )
    return result

Pyramid post-processing

After remapping, these helpers can be used to derive coarser representations for building multiresolution HEALPix pyramids.

coarsen_healpix(ds, target_level, coarsen_mode='auto', min_valid_fraction=0.5, valid_fraction=False)

Coarsen a HEALPix dataset to a lower-resolution level.

The coarsening is performed as a single reshape + reduction over all batch dimensions simultaneously — no per-slice Python loops.

Parameters:

Name	Type	Description	Default
ds	Dataset	HEALPix dataset containing a cell dimension and the attributes healpix_nside and healpix_order.	required
target_level	int	Target HEALPix level (must be lower than the current level).	required
coarsen_mode	CoarsenMode	"mean" — NaN-aware averaging for continuous fields. "mode" — most-frequent-value for categorical fields. "auto" — infer from grid_doctor_method in dataset attributes: "nearest" → mode, everything else → mean.	'auto'
min_valid_fraction	float	Minimum fraction of valid (non-NaN) children required to produce a valid parent cell. Parents with fewer valid children are set to NaN. Default 0.5 (at least half of the children must be valid).	0.5
valid_fraction	ValidFraction	Add <name>_valid_fraction variables: the fraction of valid cells of ds under each coarse cell (see create_healpix_pyramid). Fractions are relative to the level of ds, so coarsen from the finest level to get weights for exact global means.	False

Returns:

Type	Description
Dataset	Coarsened dataset.

Notes

Nested HEALPix indices have a direct parent-child relationship: pixel i at level L contains children 4*i to 4*i+3 at level L+1. Coarsening therefore reduces to grouping contiguous blocks of 4**delta_level child cells and averaging (or taking the mode for categorical data).

Ring-ordered datasets do not have contiguous parent-child layout and must be remapped directly at each target level.

With coarsen_mode="mean", always coarsen from the finest level rather than chaining level by level. A chained mean weights every valid parent equally, regardless of how many valid finest-level cells it was built from, so data with NaNs (land, sea ice, observation gaps) drifts between levels.

Raises:

Type	Description
ValueError	When the ordering is not nested, or target_level is not lower than the current level.

Source code in grid_doctor/helpers.py

def coarsen_healpix(
    ds: xr.Dataset,
    target_level: int,
    coarsen_mode: CoarsenMode = "auto",
    min_valid_fraction: float = 0.5,
    valid_fraction: ValidFraction = False,
) -> xr.Dataset:
    """Coarsen a HEALPix dataset to a lower-resolution level.

    The coarsening is performed as a single reshape + reduction over
    all batch dimensions simultaneously — no per-slice Python loops.

    Parameters
    ----------
    ds:
        HEALPix dataset containing a ``cell`` dimension and the
        attributes ``healpix_nside`` and ``healpix_order``.
    target_level:
        Target HEALPix level (must be lower than the current level).
    coarsen_mode:
        ``"mean"`` — NaN-aware averaging for continuous fields.
        ``"mode"`` — most-frequent-value for categorical fields.
        ``"auto"`` — infer from ``grid_doctor_method`` in dataset
        attributes: ``"nearest"`` → mode, everything else → mean.
    min_valid_fraction:
        Minimum fraction of valid (non-NaN) children required to
        produce a valid parent cell.  Parents with fewer valid
        children are set to NaN.  Default ``0.5`` (at least half
        of the children must be valid).
    valid_fraction:
        Add ``<name>_valid_fraction`` variables: the fraction of valid
        cells of *ds* under each coarse cell (see
        [`create_healpix_pyramid`][grid_doctor.helpers.create_healpix_pyramid]).
        Fractions are relative to the level of *ds*, so coarsen from the
        finest level to get weights for exact global means.

    Returns
    -------
    xarray.Dataset
        Coarsened dataset.

    Notes
    -----
    Nested HEALPix indices have a direct parent-child relationship:
    pixel *i* at level *L* contains children ``4*i`` to ``4*i+3`` at
    level *L+1*.  Coarsening therefore reduces to grouping contiguous
    blocks of ``4**delta_level`` child cells and averaging (or taking
    the mode for categorical data).

    Ring-ordered datasets do not have contiguous parent-child layout
    and must be remapped directly at each target level.

    With ``coarsen_mode="mean"``, always coarsen from the finest level
    rather than chaining level by level.  A chained mean weights every
    valid parent equally, regardless of how many valid finest-level
    cells it was built from, so data with NaNs (land, sea ice,
    observation gaps) drifts between levels.

    Raises
    ------
    ValueError
        When the ordering is not nested, or *target_level* is not
        lower than the current level.
    """
    current_nside = int(ds.attrs["healpix_nside"])
    target_nside = 2**target_level
    if target_nside >= current_nside:
        raise ValueError("target_level must be lower than the current HEALPix level.")

    is_nested = str(ds.attrs.get("healpix_order", "nested")) in {"nested", "nest"}
    if not is_nested:
        raise ValueError(
            "coarsen_healpix only supports nested HEALPix ordering. "
            "Use create_healpix_pyramid(..., nest=False) to "
            "regenerate ring levels directly."
        )

    current_level = int(ds.attrs.get("healpix_level", int(np.log2(current_nside))))
    delta_level = current_level - target_level
    if delta_level <= 0:
        raise ValueError("target_level must be lower than the current HEALPix level.")

    return coarsen_dataset(
        ds,
        source_level=current_level,
        target_level=target_level,
        coarsen_mode=resolve_coarsen_mode(ds, coarsen_mode),
        min_valid_fraction=min_valid_fraction,
        mode_kernel=_coarsen_array_mode,
        valid_fraction=valid_fraction,
    )