pyinterp.regrid.conservative

Contents

pyinterp.regrid.conservative#

pyinterp.regrid.conservative(data, x, y, target_x, target_y, *, axes=(0, 1), geodetic=True, spheroid=None, min_valid_fraction=0.0, normalize=True, num_threads=0)[source]#

Regrid a field between rectilinear grids, conserving its integral.

Each target cell receives the mean of the source values weighted by the area of their overlap with the cell (first order conservative regridding). This is the NumPy layer of pyinterp.backends.xarray.DataArrayAccessor.regrid().

The cells of an axis are given by their centers, a 1D array, or by their bounds, an (n, 2) array as the CF bounds variables. The bounds of cells given by their centers are the midpoints between consecutive centers (see cell_bounds()), so an axis of a single cell must be given by its bounds. Edges of contiguous cells, an array of n + 1 values, are converted with np.column_stack((edges[:-1], edges[1:])).

The areas are computed on a sphere unless spheroid is given; the cells are longitude/latitude boxes, bounded by meridians and parallels. ESMF and xESMF draw the edges of the cells as great circles instead, which bulge towards the poles by an amount proportional to the square of the width of the cell and to sin 2φ: the results differ most at mid-latitudes, not at the poles. Coarsening a smooth global field from 0.25° to 1°, they differ by up to 2·10⁻⁵ of its amplitude between 15° and 75° of latitude, 3·10⁻⁶ between 75° and 85° and 3·10⁻⁷ beyond; to 5°, by up to 4·10⁻⁴ between 15° and 75°. Next to undefined values, the slivers of the great circles reach defined values that the boxes do not: some target cells get a value from xESMF and NaN from pyinterp, and the means of partly undefined cells differ by up to 3·10⁻³ at 5°.

Parameters:

  • data (ArrayLike) – Field to regrid. Its other dimensions (time, depth…) are regridded independently.

  • x (ArrayLike) – Cells of the source grid along the first axis: longitudes in degrees if geodetic.

  • y (ArrayLike) – Cells of the source grid along the second axis: latitudes in degrees if geodetic.

  • target_x (ArrayLike) – Cells of the target grid along the first axis.

  • target_y (ArrayLike) – Cells of the target grid along the second axis.

  • axes (tuple[int, int]) – Dimensions of data along the first and second axes.

  • geodetic (bool) – If true, the first axis is a longitude, periodic with a period of 360 degrees: the source and target grids can use different conventions (e.g. [0, 360) and [-180, 180)). The second axis is a latitude, clamped to [-90, 90]. Otherwise, both axes are Cartesian and the cells are weighted by their surface.

  • spheroid (Spheroid | None) – Ellipsoid on which the areas of the cells are computed. Defaults to a sphere, as ESMF. Pass pyinterp.geometry.geographic.Spheroid () for WGS-84: the areas are then exact on the ellipsoid, using the authalic latitude. Coarsening to 1°, the results differ from those on the sphere by about 10⁻⁶ of the field; to 5°, by about 10⁻⁵.

  • min_valid_fraction (float) – Minimum fraction of the area of the target cell covered by defined source values (and by the source grid); below, the result is NaN. With the default, a target cell gets a value as soon as it overlaps a defined source value.

  • normalize (bool) – If true (conservative_normed in xESMF), the undefined source values (NaN) and the parts of the target cell outside the source grid are ignored: the result is the weighted mean of the defined values. Otherwise (conservative), the weighted sum is divided by the area of the whole target cell, as if the undefined values were zero.

  • num_threads (int) – Number of threads, 0 for all CPUs.

Returns:

The regridded field, of the same dimensions as data, with the lengths of the target axes. Its type is float32 if data is float32, float64 otherwise.

Raises:

ValueError – If the cells of an axis are invalid, or if the shape of data does not match the source grid.

Return type:

NDArray[np.floating]

Examples

Coarsen a global 0.25° field, stored as (lon, lat), to 1°:

>>> lon = np.arange(0, 360, 0.25) + 0.125
>>> lat = np.arange(-90, 90, 0.25) + 0.125
>>> field = np.ones((lon.size, lat.size))
>>> result = conservative(
...     field, lon, lat, np.arange(-179.5, 180), np.arange(-89.5, 90)
... )
>>> result.shape
(360, 180)

A single target cell, given by its bounds, on a (time, lat, lon) array:

>>> mean = conservative(
...     np.ones((3, lat.size, lon.size)),
...     lon,
...     lat,
...     [[-10.0, 10.0]],
...     [[40.0, 50.0]],
...     axes=(2, 1),
... )
>>> mean.shape
(3, 1, 1)
Contents