Skip to content

Exotic Model-Grid Adapters#

pyramids.grids turns model grids that are not row-major rasters — ORCA curvilinear ocean grids, octahedral reduced-Gaussian grids, and HEALPix sphere pixelizations — into a regular-grid Dataset.

Rather than shipping a new regridder, each adapter reshapes its grid into one of the two regridding bridges pyramids already has, then reuses it:

Adapter Input grid Bridge it reuses
from_orca curvilinear (ny, nx) lon/lat + data UGRID quad mesh → mesh_to_grid
from_octahedral ragged per-point lat/lon + values scattered points → grid_points
from_healpix per-pixel HEALPix values (nside) pixel centres → grid_points

Every adapter returns a single-band Dataset (array + geotransform + CRS) that renders with no further GIS code on the caller side.

No healpy dependency

from_healpix needs only the HEALPix pixel→centre mapping (pix2ang), which is a closed-form from the HEALPix paper. pyramids implements it in plain NumPy for both the RING and NESTED pixel orderings, so no healpy (or other HEALPix C library) is required.

See the exotic grids example notebook for runnable end-to-end usage of all three adapters.

Functions#

pyramids.grids.from_orca(lon2d, lat2d, data2d, *, cell_size, method='nearest', epsg=4326, nodata=-9999.0) #

Regrid an ORCA curvilinear field onto a regular-grid :class:Dataset.

The (ny, nx) coordinate arrays are treated as mesh nodes and stitched into (ny - 1) * (nx - 1) quadrilateral faces. Each face value is the NaN-aware mean of its four corner nodes (:func:numpy.nanmean), so every value in data2d (including the last row and column) contributes; NaN-masked nodes are ignored and a face is NaN only when all four of its corners are NaN. Use NaN (not a finite sentinel) to mark missing input values. The resulting UGRID mesh is interpolated to a regular grid via :meth:UgridDataset.to_dataset.

Parameters:

Name Type Description Default
lon2d ndarray

(ny, nx) array of node longitudes (x-coordinates).

 required
lat2d ndarray

(ny, nx) array of node latitudes (y-coordinates), same shape as lon2d.

 required
data2d ndarray

(ny, nx) array of field values aligned with the coordinate grid.

 required
cell_size float

Output pixel size in the target CRS units.

 required
method str

Interpolation method passed to mesh_to_grid ("nearest" or "linear").

 'nearest'
epsg int

Output EPSG code.

 4326
nodata float

No-data value stamped on cells the mesh does not cover.

 -9999.0

Returns:

Type Description
Dataset

A single-band :class:~pyramids.dataset.Dataset of the interpolated surface.

Raises:

Type Description
ValueError

lon2d, lat2d and data2d are not all the same 2-D shape, or the grid is smaller than 2 x 2 (no quad cells can be formed).

Examples:

  • Regrid a small curvilinear field and inspect the raster it produces: 
    >>> import numpy as np
>>> from pyramids.grids import from_orca
>>> lon2d = np.array([[0.0, 1.0, 2.0], [0.0, 1.0, 2.0]])
>>> lat2d = np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]])
>>> data2d = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
>>> ds = from_orca(lon2d, lat2d, data2d, cell_size=0.5)
>>> ds.band_count
1
>>> ds.epsg
4326
  • Mismatched coordinate/data shapes are rejected: 
    >>> import numpy as np
>>> from pyramids.grids import from_orca
>>> try:
...     from_orca(np.zeros((2, 3)), np.zeros((2, 2)), np.zeros((2, 3)), cell_size=0.5)
... except ValueError as exc:
...     print("same shape" in str(exc))
True
See Also
  • :func:pyramids.grids.from_octahedral: regrid ragged per-point fields.
  • :meth:pyramids.netcdf.ugrid.dataset.UgridDataset.to_dataset: the mesh→raster bridge this adapter delegates to.
Source code in src/pyramids/grids/orca.py 
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
def from_orca(
    lon2d: np.ndarray,
    lat2d: np.ndarray,
    data2d: np.ndarray,
    *,
    cell_size: float,
    method: str = "nearest",
    epsg: int = 4326,
    nodata: float = -9999.0,
) -> Dataset:
    """Regrid an ORCA curvilinear field onto a regular-grid :class:`Dataset`.

    The ``(ny, nx)`` coordinate arrays are treated as mesh nodes and stitched into
    ``(ny - 1) * (nx - 1)`` quadrilateral faces. Each face value is the NaN-aware mean
    of its four corner nodes (:func:`numpy.nanmean`), so every value in ``data2d``
    (including the last row and column) contributes; ``NaN``-masked nodes are ignored
    and a face is ``NaN`` only when all four of its corners are ``NaN``. Use ``NaN``
    (not a finite sentinel) to mark missing input values. The resulting UGRID mesh is
    interpolated to a regular grid via :meth:`UgridDataset.to_dataset`.

    Args:
        lon2d: ``(ny, nx)`` array of node longitudes (x-coordinates).
        lat2d: ``(ny, nx)`` array of node latitudes (y-coordinates), same shape as
            ``lon2d``.
        data2d: ``(ny, nx)`` array of field values aligned with the coordinate grid.
        cell_size: Output pixel size in the target CRS units.
        method: Interpolation method passed to ``mesh_to_grid`` ("nearest" or
            "linear").
        epsg: Output EPSG code.
        nodata: No-data value stamped on cells the mesh does not cover.

    Returns:
        A single-band :class:`~pyramids.dataset.Dataset` of the interpolated surface.

    Raises:
        ValueError: ``lon2d``, ``lat2d`` and ``data2d`` are not all the same 2-D shape,
            or the grid is smaller than ``2 x 2`` (no quad cells can be formed).

    Examples:
        - Regrid a small curvilinear field and inspect the raster it produces:
            ```python
            >>> import numpy as np
            >>> from pyramids.grids import from_orca
            >>> lon2d = np.array([[0.0, 1.0, 2.0], [0.0, 1.0, 2.0]])
            >>> lat2d = np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]])
            >>> data2d = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
            >>> ds = from_orca(lon2d, lat2d, data2d, cell_size=0.5)
            >>> ds.band_count
            1
            >>> ds.epsg
            4326

            ```
        - Mismatched coordinate/data shapes are rejected:
            ```python
            >>> import numpy as np
            >>> from pyramids.grids import from_orca
            >>> try:
            ...     from_orca(np.zeros((2, 3)), np.zeros((2, 2)), np.zeros((2, 3)), cell_size=0.5)
            ... except ValueError as exc:
            ...     print("same shape" in str(exc))
            True

            ```

    See Also:
        - :func:`pyramids.grids.from_octahedral`: regrid ragged per-point fields.
        - :meth:`pyramids.netcdf.ugrid.dataset.UgridDataset.to_dataset`: the mesh→raster
          bridge this adapter delegates to.
    """
    lon2d = np.asarray(lon2d, dtype=np.float64)
    lat2d = np.asarray(lat2d, dtype=np.float64)
    data2d = np.asarray(data2d, dtype=np.float64)
    if not (lon2d.shape == lat2d.shape == data2d.shape):
        raise ValueError(
            "lon2d, lat2d and data2d must share the same shape; got "
            f"{lon2d.shape}, {lat2d.shape}, {data2d.shape}."
        )
    if lon2d.ndim != 2:
        raise ValueError(f"ORCA inputs must be 2-D; got {lon2d.ndim}-D arrays.")
    ny, nx = lon2d.shape
    if ny < 2 or nx < 2:
        raise ValueError(
            f"ORCA grid must be at least 2 x 2 to form quad cells; got {ny} x {nx}."
        )

    node_x = lon2d.ravel()
    node_y = lat2d.ravel()
    idx = np.arange(ny * nx).reshape(ny, nx)
    faces = np.stack(
        [
            idx[:-1, :-1].ravel(),
            idx[:-1, 1:].ravel(),
            idx[1:, 1:].ravel(),
            idx[1:, :-1].ravel(),
        ],
        axis=1,
    )
    # Each face value is the NaN-aware mean of its four corner nodes, so every node in
    # data2d contributes (no last row/column dropped) and a NaN-masked node only blanks
    # a face when all four of its corners are NaN.
    corners = np.stack(
        [data2d[:-1, :-1], data2d[:-1, 1:], data2d[1:, 1:], data2d[1:, :-1]], axis=0
    )
    with warnings.catch_warnings():
        # All-NaN faces are intended (fully-masked cell -> NaN); suppress the warning.
        warnings.simplefilter("ignore", category=RuntimeWarning)
        face_values = np.nanmean(corners, axis=0)

    mesh = UgridDataset.create_from_arrays(
        node_x,
        node_y,
        faces,
        data={"z": face_values.ravel()},
        data_locations={"z": "face"},
        epsg=epsg,
    )
    result = mesh.to_dataset(
        "z", cell_size=cell_size, method=method, epsg=epsg, nodata=nodata
    )
    return result

pyramids.grids.from_octahedral(lats, lons, values, *, cell_size, algorithm='nearest', epsg=4326, bbox=None) #

Regrid an octahedral reduced-Gaussian field onto a regular-grid :class:Dataset.

The per-point lats/lons/values triples are wrapped in a point :class:~geopandas.GeoDataFrame and interpolated with gdal.Grid via :func:~pyramids.dataset.ops.interpolate.grid_points.

Parameters:

Name Type Description Default
lats ndarray

1-D array of point latitudes.

 required
lons ndarray

1-D array of point longitudes, same length as lats.

 required
values ndarray

1-D array of field values, same length as lats.

 required
cell_size float

Output pixel size in the target CRS units.

 required
algorithm str

A gdal.Grid algorithm string (e.g. "nearest", "invdist:power=2.0:smoothing=0.0", "linear").

 'nearest'
epsg int

Output EPSG code.

 4326
bbox tuple[float, float, float, float] | None

Optional (minx, miny, maxx, maxy) output extent in the target CRS. Defaults to the points' bounding box; pass e.g. (-180, -90, 180, 90) to pin a fixed global grid.

 None

Returns:

Type Description
Dataset

A single-band :class:~pyramids.dataset.Dataset of the interpolated surface.

Raises:

Type Description
ValueError

lats, lons and values are not 1-D arrays of equal length.

Examples:

  • Grid four corner observations with nearest-neighbour and inspect the result: 
    >>> import numpy as np
>>> from pyramids.grids import from_octahedral
>>> lats = np.array([0.0, 0.0, 5.0, 5.0])
>>> lons = np.array([0.0, 5.0, 0.0, 5.0])
>>> values = np.array([1.0, 2.0, 3.0, 4.0])
>>> ds = from_octahedral(lats, lons, values, cell_size=1.0, algorithm="nearest")
>>> (ds.rows, ds.columns, ds.band_count)
(5, 5, 1)
  • Arrays of unequal length are rejected: 
    >>> import numpy as np
>>> from pyramids.grids import from_octahedral
>>> try:
...     from_octahedral(np.zeros(4), np.zeros(3), np.zeros(4), cell_size=1.0)
... except ValueError as exc:
...     print("equal length" in str(exc))
True
See Also
  • :func:pyramids.grids.from_orca: regrid curvilinear (ny, nx) fields.
  • :func:pyramids.dataset.ops.interpolate.grid_points: the scattered-point interpolation this adapter delegates to.
Source code in src/pyramids/grids/octahedral.py 
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
def from_octahedral(
    lats: np.ndarray,
    lons: np.ndarray,
    values: np.ndarray,
    *,
    cell_size: float,
    algorithm: str = "nearest",
    epsg: int = 4326,
    bbox: tuple[float, float, float, float] | None = None,
) -> Dataset:
    """Regrid an octahedral reduced-Gaussian field onto a regular-grid :class:`Dataset`.

    The per-point ``lats``/``lons``/``values`` triples are wrapped in a point
    :class:`~geopandas.GeoDataFrame` and interpolated with ``gdal.Grid`` via
    :func:`~pyramids.dataset.ops.interpolate.grid_points`.

    Args:
        lats: 1-D array of point latitudes.
        lons: 1-D array of point longitudes, same length as ``lats``.
        values: 1-D array of field values, same length as ``lats``.
        cell_size: Output pixel size in the target CRS units.
        algorithm: A ``gdal.Grid`` algorithm string (e.g. ``"nearest"``,
            ``"invdist:power=2.0:smoothing=0.0"``, ``"linear"``).
        epsg: Output EPSG code.
        bbox: Optional ``(minx, miny, maxx, maxy)`` output extent in the target CRS.
            Defaults to the points' bounding box; pass e.g. ``(-180, -90, 180, 90)``
            to pin a fixed global grid.

    Returns:
        A single-band :class:`~pyramids.dataset.Dataset` of the interpolated surface.

    Raises:
        ValueError: ``lats``, ``lons`` and ``values`` are not 1-D arrays of equal
            length.

    Examples:
        - Grid four corner observations with nearest-neighbour and inspect the result:
            ```python
            >>> import numpy as np
            >>> from pyramids.grids import from_octahedral
            >>> lats = np.array([0.0, 0.0, 5.0, 5.0])
            >>> lons = np.array([0.0, 5.0, 0.0, 5.0])
            >>> values = np.array([1.0, 2.0, 3.0, 4.0])
            >>> ds = from_octahedral(lats, lons, values, cell_size=1.0, algorithm="nearest")
            >>> (ds.rows, ds.columns, ds.band_count)
            (5, 5, 1)

            ```
        - Arrays of unequal length are rejected:
            ```python
            >>> import numpy as np
            >>> from pyramids.grids import from_octahedral
            >>> try:
            ...     from_octahedral(np.zeros(4), np.zeros(3), np.zeros(4), cell_size=1.0)
            ... except ValueError as exc:
            ...     print("equal length" in str(exc))
            True

            ```

    See Also:
        - :func:`pyramids.grids.from_orca`: regrid curvilinear ``(ny, nx)`` fields.
        - :func:`pyramids.dataset.ops.interpolate.grid_points`: the scattered-point
          interpolation this adapter delegates to.
    """
    lats = np.asarray(lats, dtype=np.float64).ravel()
    lons = np.asarray(lons, dtype=np.float64).ravel()
    values = np.asarray(values, dtype=np.float64).ravel()
    if not (lats.size == lons.size == values.size):
        raise ValueError(
            "lats, lons and values must have equal length; got "
            f"{lats.size}, {lons.size}, {values.size}."
        )

    gdf = GeoDataFrame(
        {"z": values},
        geometry=points_from_xy(lons, lats),
        crs=epsg,
    )
    result = grid_points(
        gdf,
        "z",
        Dataset,
        algorithm=algorithm,
        cell_size=cell_size,
        bbox=bbox,
        epsg=epsg,
    )
    return result

pyramids.grids.from_healpix(values, *, nside=None, nest=False, cell_size, method='nearest', epsg=4326, bbox=None) #

Regrid a HEALPix field onto a regular-grid :class:Dataset.

Each HEALPix pixel's centre longitude/latitude is computed in plain NumPy (no healpy dependency) and the resulting points are interpolated with gdal.Grid via :func:~pyramids.dataset.ops.interpolate.grid_points.

Parameters:

Name Type Description Default
values ndarray

1-D array of per-pixel HEALPix values, length 12 * nside**2.

 required
nside int | None

HEALPix resolution parameter. Derived from len(values) when omitted.

 None
nest bool

True for NESTED pixel ordering, False (default) for RING.

 False
cell_size float

Output pixel size in the target CRS units (degrees for EPSG:4326).

 required
method str

A gdal.Grid algorithm string (e.g. "nearest", "linear", "invdist:power=2.0:smoothing=0.0").

 'nearest'
epsg int

Output EPSG code.

 4326
bbox tuple[float, float, float, float] | None

Optional (minx, miny, maxx, maxy) output extent in the target CRS. Defaults to the pixel-centres' bounding box; pass e.g. (-180, -90, 180, 90) to pin a fixed global grid.

 None

Returns:

Type Description
Dataset

A single-band :class:~pyramids.dataset.Dataset of the interpolated surface.

Raises:

Type Description
ValueError

values is not 1-D; len(values) is not a valid HEALPix pixel count (12 * nside**2); nside disagrees with len(values); or nest=True with an nside that is not a power of two.

Examples:

  • Regrid a synthetic nside=1 field (12 pixels) and inspect the raster: 
    >>> import numpy as np
>>> from pyramids.grids import from_healpix
>>> ds = from_healpix(np.arange(12.0), cell_size=30.0)
>>> ds.band_count
1
>>> ds.epsg
4326
  • NESTED ordering is supported and requires a power-of-two nside: 
    >>> import numpy as np
>>> from pyramids.grids import from_healpix
>>> ds = from_healpix(np.arange(48.0), nside=2, nest=True, cell_size=20.0)
>>> ds.band_count
1
  • A length that is not a valid HEALPix pixel count is rejected: 
    >>> import numpy as np
>>> from pyramids.grids import from_healpix
>>> try:
...     from_healpix(np.zeros(10), cell_size=30.0)
... except ValueError as exc:
...     print("valid HEALPix pixel count" in str(exc))
True
See Also
  • :func:pyramids.grids.from_octahedral: the sibling point-based adapter this function delegates to via grid_points.
Source code in src/pyramids/grids/healpix.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
def from_healpix(
    values: np.ndarray,
    *,
    nside: int | None = None,
    nest: bool = False,
    cell_size: float,
    method: str = "nearest",
    epsg: int = 4326,
    bbox: tuple[float, float, float, float] | None = None,
) -> Dataset:
    """Regrid a HEALPix field onto a regular-grid :class:`Dataset`.

    Each HEALPix pixel's centre longitude/latitude is computed in plain NumPy (no
    `healpy` dependency) and the resulting points are interpolated with `gdal.Grid`
    via :func:`~pyramids.dataset.ops.interpolate.grid_points`.

    Args:
        values: 1-D array of per-pixel HEALPix values, length `12 * nside**2`.
        nside: HEALPix resolution parameter. Derived from `len(values)` when omitted.
        nest: `True` for NESTED pixel ordering, `False` (default) for RING.
        cell_size: Output pixel size in the target CRS units (degrees for EPSG:4326).
        method: A `gdal.Grid` algorithm string (e.g. `"nearest"`, `"linear"`,
            `"invdist:power=2.0:smoothing=0.0"`).
        epsg: Output EPSG code.
        bbox: Optional `(minx, miny, maxx, maxy)` output extent in the target CRS.
            Defaults to the pixel-centres' bounding box; pass e.g.
            `(-180, -90, 180, 90)` to pin a fixed global grid.

    Returns:
        A single-band :class:`~pyramids.dataset.Dataset` of the interpolated surface.

    Raises:
        ValueError: `values` is not 1-D; `len(values)` is not a valid HEALPix pixel
            count (`12 * nside**2`); `nside` disagrees with `len(values)`; or
            `nest=True` with an `nside` that is not a power of two.

    Examples:
        - Regrid a synthetic `nside=1` field (12 pixels) and inspect the raster:
            ```python
            >>> import numpy as np
            >>> from pyramids.grids import from_healpix
            >>> ds = from_healpix(np.arange(12.0), cell_size=30.0)
            >>> ds.band_count
            1
            >>> ds.epsg
            4326

            ```
        - NESTED ordering is supported and requires a power-of-two `nside`:
            ```python
            >>> import numpy as np
            >>> from pyramids.grids import from_healpix
            >>> ds = from_healpix(np.arange(48.0), nside=2, nest=True, cell_size=20.0)
            >>> ds.band_count
            1

            ```
        - A length that is not a valid HEALPix pixel count is rejected:
            ```python
            >>> import numpy as np
            >>> from pyramids.grids import from_healpix
            >>> try:
            ...     from_healpix(np.zeros(10), cell_size=30.0)
            ... except ValueError as exc:
            ...     print("valid HEALPix pixel count" in str(exc))
            True

            ```

    See Also:
        - :func:`pyramids.grids.from_octahedral`: the sibling point-based adapter this
          function delegates to via `grid_points`.
    """
    values = np.asarray(values, dtype=np.float64).ravel()
    npix = values.size

    if nside is None:
        derived = int(math.isqrt(npix // 12))
        if derived < 1 or 12 * derived * derived != npix:
            raise ValueError(
                f"values length {npix} is not a valid HEALPix pixel count "
                "(12 * nside**2); pass an explicit nside if this is intentional."
            )
        nside = derived
    elif nside < 1 or 12 * nside * nside != npix:
        raise ValueError(
            f"nside={nside} implies {12 * (nside or 0) ** 2} pixels but values has "
            f"{npix}; they must be a valid HEALPix pair (12 * nside**2)."
        )

    if nest and (nside & (nside - 1)) != 0:
        raise ValueError(
            f"NESTED ordering requires nside to be a power of two; got nside={nside}."
        )

    lon, lat = _pix2lonlat(nside, np.arange(npix), nest)
    gdf = GeoDataFrame(
        {"z": values},
        geometry=points_from_xy(lon, lat),
        crs=epsg,
    )
    result = grid_points(
        gdf, "z", Dataset, algorithm=method, cell_size=cell_size, bbox=bbox, epsg=epsg
    )
    return result