Grid Operations

bathy.clip(data, geometry=None, region=None)

Clip a grid to a polygon boundary or region preset.

Parameters:

Name	Type	Description	Default
data	DataArray	Elevation grid.	required
geometry	GeoDataFrame, str, or Path	Clipping geometry. Can be a GeoDataFrame, or a path to a vector file (Shapefile, GeoPackage, GeoJSON, etc.).	None
region	str	Named region preset (see :func:bathy.list_regions). Cannot be used together with geometry.	None

Returns:

Type	Description
DataArray	Clipped grid with dims normalised to bathy conventions.

Raises:

Type	Description
ValueError	If both or neither of geometry and region are specified.

Examples:

Clip to a region preset:

>>> clipped = bathy.clip(data, region="mediterranean")

Clip to a custom polygon:

>>> from shapely.geometry import box
>>> gdf = gpd.GeoDataFrame(geometry=[box(-8, 51, -7, 53)], crs="EPSG:4326")
>>> clipped = bathy.clip(data, geometry=gdf)

Source code in src/bathy/grid.py

def clip(
    data: xr.DataArray,
    geometry: gpd.GeoDataFrame | str | Path | None = None,
    region: str | None = None,
) -> xr.DataArray:
    """Clip a grid to a polygon boundary or region preset.

    Parameters
    ----------
    data : xr.DataArray
        Elevation grid.
    geometry : GeoDataFrame, str, or Path, optional
        Clipping geometry. Can be a GeoDataFrame, or a path to a vector
        file (Shapefile, GeoPackage, GeoJSON, etc.).
    region : str, optional
        Named region preset (see :func:`bathy.list_regions`).
        Cannot be used together with *geometry*.

    Returns
    -------
    xr.DataArray
        Clipped grid with dims normalised to bathy conventions.

    Raises
    ------
    ValueError
        If both or neither of *geometry* and *region* are specified.

    Examples
    --------
    Clip to a region preset:

    >>> clipped = bathy.clip(data, region="mediterranean")

    Clip to a custom polygon:

    >>> from shapely.geometry import box
    >>> gdf = gpd.GeoDataFrame(geometry=[box(-8, 51, -7, 53)], crs="EPSG:4326")
    >>> clipped = bathy.clip(data, geometry=gdf)
    """
    if (geometry is None) == (region is None):
        raise ValueError("Specify exactly one of 'geometry' or 'region'.")

    if region is not None:
        gdf = _region_to_geodataframe(region)
    elif isinstance(geometry, (str, Path)):
        gdf = gpd.read_file(geometry)
    else:
        # geometry is guaranteed non-None by the check above
        gdf: gpd.GeoDataFrame = geometry  # ty: ignore[invalid-assignment]

    data = _prepare_spatial(data)

    # Reproject clipping geometry to match data CRS if needed
    if gdf.crs is not None and not gdf.crs.equals(data.rio.crs):
        gdf = gdf.to_crs(data.rio.crs)

    clipped = data.rio.clip(gdf.geometry, drop=True)
    return _normalize_dims(clipped)

bathy.resample(data, resolution_degrees=None, resolution_m=None, method='bilinear')

Resample a grid to a new resolution.

Parameters:

Name	Type	Description	Default
data	DataArray	Elevation grid.	required
resolution_degrees	float	Target resolution in degrees. For projected CRS this is converted to metres. Exactly one of resolution_degrees or resolution_m must be given.	None
resolution_m	float	Target resolution in metres. For geographic CRS this is converted to degrees using the grid's centre latitude.	None
method	str	Resampling method. One of 'nearest', 'bilinear', 'cubic', 'cubic_spline', 'lanczos', 'average', 'min', 'max'. Default 'bilinear'.	'bilinear'

Returns:

Type	Description
DataArray	Resampled grid with dims normalised to bathy conventions.

Examples:

Resample to approximately 500 m:

>>> resampled = bathy.resample(data, resolution_m=500)

Resample to 0.01°:

>>> coarser = bathy.resample(data, resolution_degrees=0.01)

Source code in src/bathy/grid.py

def resample(
    data: xr.DataArray,
    resolution_degrees: float | None = None,
    resolution_m: float | None = None,
    method: str = "bilinear",
) -> xr.DataArray:
    """Resample a grid to a new resolution.

    Parameters
    ----------
    data : xr.DataArray
        Elevation grid.
    resolution_degrees : float, optional
        Target resolution in degrees. For projected CRS this is
        converted to metres. Exactly one of *resolution_degrees* or
        *resolution_m* must be given.
    resolution_m : float, optional
        Target resolution in metres. For geographic CRS this is
        converted to degrees using the grid's centre latitude.
    method : str, optional
        Resampling method. One of ``'nearest'``, ``'bilinear'``,
        ``'cubic'``, ``'cubic_spline'``, ``'lanczos'``, ``'average'``,
        ``'min'``, ``'max'``. Default ``'bilinear'``.

    Returns
    -------
    xr.DataArray
        Resampled grid with dims normalised to bathy conventions.

    Examples
    --------
    Resample to approximately 500 m:

    >>> resampled = bathy.resample(data, resolution_m=500)

    Resample to 0.01°:

    >>> coarser = bathy.resample(data, resolution_degrees=0.01)
    """
    if (resolution_degrees is None) == (resolution_m is None):
        raise ValueError(
            "Specify exactly one of 'resolution_degrees' or 'resolution_m'."
        )

    if resolution_m is not None:
        resolution = _metres_to_resolution(data, resolution_m)
    else:
        resolution = resolution_degrees

    data = _prepare_spatial(data)
    resampled = data.rio.reproject(
        data.rio.crs,
        resolution=resolution,
        resampling=_get_resampling(method),
    )
    return _normalize_dims(resampled)

bathy.reproject(data, target_crs, resolution=None, method='bilinear')

Reproject a grid to a different coordinate reference system.

Dim names are automatically updated: lon/lat for geographic CRS, x/y for projected CRS.

Parameters:

Name	Type	Description	Default
data	DataArray	Elevation grid.	required
target_crs	str, int, or pyproj.CRS	Target CRS (e.g. "EPSG:32629").	required
resolution	float	Target resolution in the units of target_crs. If None, rioxarray picks an appropriate resolution automatically.	None
method	str	Resampling method (default 'bilinear'). See :func:resample.	'bilinear'

Returns:

Type	Description
DataArray	Reprojected grid with dim names matching bathy conventions.

Examples:

Reproject to UTM zone 29N:

>>> utm = bathy.reproject(data, target_crs="EPSG:32629")
>>> utm.dims  # ('y', 'x') — projected convention

Source code in src/bathy/grid.py

def reproject(
    data: xr.DataArray,
    target_crs: str | int | CRS,
    resolution: float | None = None,
    method: str = "bilinear",
) -> xr.DataArray:
    """Reproject a grid to a different coordinate reference system.

    Dim names are automatically updated: ``lon``/``lat`` for geographic
    CRS, ``x``/``y`` for projected CRS.

    Parameters
    ----------
    data : xr.DataArray
        Elevation grid.
    target_crs : str, int, or pyproj.CRS
        Target CRS (e.g. ``"EPSG:32629"``).
    resolution : float, optional
        Target resolution in the units of *target_crs*. If ``None``,
        rioxarray picks an appropriate resolution automatically.
    method : str, optional
        Resampling method (default ``'bilinear'``). See :func:`resample`.

    Returns
    -------
    xr.DataArray
        Reprojected grid with dim names matching bathy conventions.

    Examples
    --------
    Reproject to UTM zone 29N:

    >>> utm = bathy.reproject(data, target_crs="EPSG:32629")
    >>> utm.dims  # ('y', 'x') — projected convention
    """
    data = _prepare_spatial(data)

    kwargs: dict = {"resampling": _get_resampling(method)}
    if resolution is not None:
        kwargs["resolution"] = resolution

    reprojected = data.rio.reproject(target_crs, **kwargs)
    return _normalize_dims(reprojected)

bathy.merge(datasets, method='mean', resolution=None)

Merge multiple grids into a single grid.

Overlapping cells are combined using method. When grids have different resolutions they are resampled to the finest resolution (or to resolution if given) before merging.

Parameters:

Name	Type	Description	Default
datasets	list of xr.DataArray	Grids to merge. Must share the same CRS.	required
method	str	How to combine overlapping cells: 'mean', 'min', 'max', or 'first'. Default 'mean'.	'mean'
resolution	float	Target resolution. If None, the finest input resolution is used.	None

Returns:

Type	Description
DataArray	Merged grid with dims normalised to bathy conventions.

Raises:

Type	Description
ValueError	If fewer than two datasets are given, CRS do not match, or method is unknown.

Examples:

Merge two adjacent tiles, averaging overlapping cells:

>>> merged = bathy.merge([tile_west, tile_east])

Merge keeping the minimum depth in overlaps:

>>> merged = bathy.merge([gebco, emodnet], method="min")

Source code in src/bathy/grid.py

def merge(
    datasets: list[xr.DataArray],
    method: str = "mean",
    resolution: float | None = None,
) -> xr.DataArray:
    """Merge multiple grids into a single grid.

    Overlapping cells are combined using *method*. When grids have
    different resolutions they are resampled to the finest resolution
    (or to *resolution* if given) before merging.

    Parameters
    ----------
    datasets : list of xr.DataArray
        Grids to merge. Must share the same CRS.
    method : str, optional
        How to combine overlapping cells: ``'mean'``, ``'min'``,
        ``'max'``, or ``'first'``. Default ``'mean'``.
    resolution : float, optional
        Target resolution. If ``None``, the finest input resolution is
        used.

    Returns
    -------
    xr.DataArray
        Merged grid with dims normalised to bathy conventions.

    Raises
    ------
    ValueError
        If fewer than two datasets are given, CRS do not match, or
        *method* is unknown.

    Examples
    --------
    Merge two adjacent tiles, averaging overlapping cells:

    >>> merged = bathy.merge([tile_west, tile_east])

    Merge keeping the minimum depth in overlaps:

    >>> merged = bathy.merge([gebco, emodnet], method="min")
    """
    if len(datasets) < 2:
        raise ValueError("Need at least 2 datasets to merge.")

    valid_methods = {"mean", "min", "max", "first"}
    if method not in valid_methods:
        raise ValueError(
            f"Unknown method '{method}'. Choose from: {sorted(valid_methods)}"
        )

    # Determine common CRS
    crs_list = [get_crs(d) for d in datasets]
    ref_crs = crs_list[0] or CRS.from_epsg(4326)
    for i, c in enumerate(crs_list[1:], 1):
        effective = c or CRS.from_epsg(4326)
        if not effective.equals(ref_crs):
            raise ValueError(
                f"CRS mismatch: dataset 0 has {ref_crs}, dataset {i} has "
                f"{effective}. Reproject grids to a common CRS first."
            )

    # Determine target resolution (finest among inputs)
    if resolution is None:
        resolutions = []
        for d in datasets:
            x_dim, _ = get_dim_names(d)
            resolutions.append(float(np.abs(np.diff(d[x_dim].values).mean())))
        resolution = min(resolutions)

    # Resample all datasets to common resolution
    use_projected = is_projected(datasets[0])
    aligned = []
    for d in datasets:
        d = _prepare_spatial(d)
        x_dim, _ = get_dim_names(d)
        current_res = float(np.abs(np.diff(d[x_dim].values).mean()))
        if not np.isclose(current_res, resolution, rtol=0.01):
            if use_projected:
                d = resample(d, resolution_m=resolution)
            else:
                d = resample(d, resolution_degrees=resolution)
        aligned.append(d)

    # Build common coordinate grid spanning all inputs
    x_dim, y_dim = get_dim_names(aligned[0])
    all_x = np.concatenate([d[x_dim].values for d in aligned])
    all_y = np.concatenate([d[y_dim].values for d in aligned])
    new_x = np.arange(
        float(all_x.min()), float(all_x.max()) + resolution / 2, resolution
    )
    new_y = np.arange(
        float(all_y.min()), float(all_y.max()) + resolution / 2, resolution
    )

    # Reindex each dataset onto the common grid, then combine
    reindexed = []
    for d in aligned:
        d_x, d_y = get_dim_names(d)
        d = d.reindex(
            {d_x: new_x, d_y: new_y}, method="nearest", tolerance=resolution * 0.6
        )
        reindexed.append(d)

    stacked = xr.concat(reindexed, dim="__source__")

    if method == "mean":
        result = stacked.mean(dim="__source__")
    elif method == "min":
        result = stacked.min(dim="__source__")
    elif method == "max":
        result = stacked.max(dim="__source__")
    else:  # first
        result = stacked.isel(__source__=0)
        for i in range(1, len(reindexed)):
            mask = result.isnull()
            result = result.where(~mask, stacked.isel(__source__=i))

    result.name = "elevation"
    if ref_crs is not None:
        result = result.rio.write_crs(ref_crs)
    return _normalize_dims(result)

bathy.fill_gaps(data, method='nearest')

Fill NaN gaps in a grid using interpolation.

Parameters:

Name	Type	Description	Default
data	DataArray	Elevation grid with NaN gaps.	required
method	str	Interpolation method: 'nearest' or 'linear'. Default 'nearest'.	'nearest'

Returns:

Type	Description
DataArray	Grid with NaN values filled.

Raises:

Type	Description
ValueError	If method is not recognised.

Examples:

>>> filled = bathy.fill_gaps(data, method="linear")

Source code in src/bathy/grid.py

def fill_gaps(
    data: xr.DataArray,
    method: str = "nearest",
) -> xr.DataArray:
    """Fill NaN gaps in a grid using interpolation.

    Parameters
    ----------
    data : xr.DataArray
        Elevation grid with NaN gaps.
    method : str, optional
        Interpolation method: ``'nearest'`` or ``'linear'``.
        Default ``'nearest'``.

    Returns
    -------
    xr.DataArray
        Grid with NaN values filled.

    Raises
    ------
    ValueError
        If *method* is not recognised.

    Examples
    --------
    >>> filled = bathy.fill_gaps(data, method="linear")
    """
    valid = {"nearest", "linear"}
    if method not in valid:
        raise ValueError(f"Unknown method '{method}'. Choose from: {sorted(valid)}")

    x_dim, y_dim = get_dim_names(data)
    z = data.values.astype(float)
    mask = np.isnan(z)

    if not mask.any():
        return data

    if not (~mask).any():
        raise ValueError(
            "Cannot fill gaps: all values are NaN, no valid data to interpolate from"
        )

    xs = data[x_dim].values
    ys = data[y_dim].values
    xx, yy = np.meshgrid(xs, ys)

    known = ~mask
    filled = griddata(
        (xx[known], yy[known]),
        z[known],
        (xx, yy),
        method=method,
    )

    return data.copy(data=filled)