Coastal & Marine Spatial Analysis Pipelines

Clipping Marine Rasters to a Coastline Mask

Every coastal raster product gets clipped to the water at some point: a bathymetric grid trimmed at the shoreline, a current field masked to the model domain, a habitat suitability surface cut to the survey footprint. The operation looks like a one-liner, and three of its defaults are wrong for marine data in ways that are invisible in a plot. This page sits under coastal vector geometry processing and covers the masking side of the same boundary.

Root cause: a polygon boundary does not fall on cell edges

A raster is a partition of space into cells; a mask polygon is an arbitrary curve. Where the curve crosses a cell, that cell is partly in and partly out, and the masking function has to decide. rasterio.mask exposes the decision as all_touched: False (the default) keeps a cell only if the polygon covers its centre; True keeps a cell if the polygon touches it at all.

For an inland catchment with a 30 m grid, the difference is a rounding detail. For a coastal product it is a systematic bias of one cell ring along the entire coastline, and coastlines are long. On a 1 m bathymetric grid over 400 km of shoreline, the two settings differ by roughly 0.4 km² — and the difference is not noise, it is entirely on one side.

import geopandas as gpd
import rasterio
from rasterio.mask import mask

shore = gpd.read_file("shoreline_mhw.gpkg").to_crs("EPSG:32618")
with rasterio.open("bathy_1m.tif") as src:
    strict, _ = mask(src, shore.geometry, crop=True, all_touched=False)
    loose, _ = mask(src, shore.geometry, crop=True, all_touched=True)

print((strict != src.nodata).sum(), (loose != src.nodata).sum())
# 1_884_221_004   1_884_612_889   -> 391,885 cells, all of them at the coast

Which is right depends on what the product is for. A navigation surface wants all_touched=False on the water polygon, so a cell that is only partly water is excluded and the surface never claims depth where there may be land. An area or volume statistic wants neither: it wants fractional coverage weights, which rasterio.features.rasterize can approximate by supersampling but which neither flag provides. What matters is that the choice is deliberate and recorded, because it is not recoverable from the output.

Three Answers to “Is This Cell Water?” A raster grid of square cells is crossed diagonally by a shoreline curve. Under all_touched set to false, only cells whose centres fall on the water side are shaded, leaving a stepped edge inside the true shoreline. Under all_touched set to true, every cell the curve passes through is shaded, leaving a stepped edge outside the true shoreline. A third panel shades each boundary cell in proportion to the fraction of its area that is water, which is what an area or volume statistic actually needs and which neither flag provides. The three coastal areas differ by a full ring of cells along the whole coast. all_touched=False centre must be wet → edge sits inside the shore all_touched=True any overlap counts → edge sits outside the shore fractional coverage what an area statistic needs → neither flag gives it Over 400 km of 1 m coastline the first two differ by about 0.4 km², all of it on one side. That is not rounding. Record which policy produced a given product; it cannot be recovered from the output raster.

Step-by-step fix with production code

1. Reproject the mask, never the raster

Reprojecting a raster resamples it: values are interpolated onto a new grid and the original cell values no longer exist. Reprojecting the mask polygon is exact — it moves vertices, and the raster grid is untouched. Always move the vector.

def align_mask(mask_gdf: gpd.GeoDataFrame, raster_crs) -> gpd.GeoDataFrame:
    """Bring a mask into the raster's CRS. Never the other way round."""
    if mask_gdf.crs is None:
        raise ValueError("mask has no CRS; assign the correct one before clipping")
    if mask_gdf.crs == raster_crs:
        return mask_gdf
    logger.info("reprojecting mask %s → %s", mask_gdf.crs.to_string(), raster_crs.to_string())
    return mask_gdf.to_crs(raster_crs)

2. Read the window, not the file

rasterio.mask.mask with crop=True still reads the full extent before cropping unless it is given a windowed dataset. On a terabyte mosaic clipped to a single estuary, computing the window first turns an impossible read into a small one.

import numpy as np
import rasterio
from rasterio.features import geometry_mask
from rasterio.windows import from_bounds


def clip_windowed(
    raster_path: str,
    mask_gdf: gpd.GeoDataFrame,
    *,
    all_touched: bool = False,
) -> tuple[np.ndarray, dict]:
    """Clip to a mask, reading only the mask's window from the raster."""
    with rasterio.open(raster_path) as src:
        if src.nodata is None:
            raise ValueError(
                f"{raster_path} declares no nodata value; masking would have to "
                "invent one, and 0 is a valid depth"
            )
        aligned = align_mask(mask_gdf, src.crs)
        window = from_bounds(*aligned.total_bounds, transform=src.transform)
        window = window.round_offsets().round_lengths()
        data = src.read(1, window=window)
        transform = src.window_transform(window)

        inside = geometry_mask(
            aligned.geometry, out_shape=data.shape, transform=transform,
            invert=True, all_touched=all_touched,
        )
        out = np.where(inside, data, src.nodata).astype(data.dtype)

        profile = src.profile | {
            "height": data.shape[0], "width": data.shape[1], "transform": transform,
        }
        logger.info(
            "clipped %s: %d of %d cells retained", raster_path,
            int(inside.sum()), inside.size,
        )
        return out, profile

The nodata is None guard is the third default that bites. Masking has to write something into the excluded cells, and libraries that default to zero produce a bathymetric grid in which every land cell reads as a depth of exactly 0.0 m — a valid, plausible, catastrophic value that will be averaged into every statistic and contoured as a shoal.

3. Keep the grid, and prove you kept it

The other silent corruption in clipping is grid drift. Any operation that reprojects, resamples or “snaps” the output moves the cell centres, and a mosaic assembled from tiles clipped at slightly different origins will not align — the seams appear as a one-cell offset that looks like a real terrain step.

The defence is to derive the output transform from the source transform and a window, as the function above does, rather than from the mask’s bounds. A window is an integer offset into the existing grid, so the output’s cell centres are a subset of the input’s by construction. Building a transform from mask.total_bounds instead produces an origin at an arbitrary real coordinate and every tile lands on its own lattice.

def assert_grid_preserved(src_transform, out_transform, cell_tol: float = 1e-9) -> None:
    """Output cells must lie on the source lattice, not merely near it."""
    if abs(src_transform.a - out_transform.a) > cell_tol or \
       abs(src_transform.e - out_transform.e) > cell_tol:
        raise ValueError("cell size changed: the clip resampled the raster")
    dx = (out_transform.c - src_transform.c) / src_transform.a
    dy = (out_transform.f - src_transform.f) / src_transform.e
    if abs(dx - round(dx)) > 1e-6 or abs(dy - round(dy)) > 1e-6:
        raise ValueError(
            f"output origin is {dx - round(dx):.6f}, {dy - round(dy):.6f} cells off "
            "the source lattice; tiles clipped this way will not mosaic"
        )

This assertion costs microseconds and catches a class of defect whose usual discovery route is a visible seam in a published mosaic.

Window-Derived Tiles Mosaic; Bounds-Derived Tiles Do Not Upper row: two adjacent clipped tiles whose transforms were derived from integer windows into the source raster. Their cell boundaries line up exactly across the join and the mosaic shows a continuous grid. Lower row: the same two tiles with transforms built from each mask's own bounding box. Each tile's grid starts at an arbitrary offset, so at the join the two lattices are shifted by a fraction of a cell; a magnified inset shows the resulting sliver of overlap and gap along the seam, and a note explains that a resampled mosaic fills it by interpolation and renders it as a low ridge that looks like terrain. window-derived — same lattice tile A ends and tile B begins on a cell edge → mosaic is exact, no resampling needed bounds-derived — two lattices tile B’s grid starts 0.39 of a cell late the mosaic must resample to join them → a sub-cell ridge along every seam The ridge is smooth, continuous and the width of one cell — which is exactly what a real sand wave looks like.

Verification and acceptance test

The retained cell count and the mask polygon’s area are two independent measurements of the same thing, and they should agree to within the boundary ring.

def test_clip_area_matches_mask_area():
    shore = gpd.read_file("shoreline_mhw.gpkg").to_crs("EPSG:32618")
    data, profile = clip_windowed("bathy_1m.tif", shore, all_touched=False)

    cell_area = abs(profile["transform"].a * profile["transform"].e)
    retained = int((data != profile["nodata"]).sum())
    raster_area = retained * cell_area
    polygon_area = float(shore.geometry.area.sum())

    # The boundary ring is the expected disagreement: perimeter × cell size.
    perimeter = float(shore.geometry.length.sum())
    ring = perimeter * abs(profile["transform"].a)
    assert abs(raster_area - polygon_area) < ring, (
        f"raster {raster_area:,.0f} m² vs polygon {polygon_area:,.0f} m² "
        f"differs by more than one boundary ring ({ring:,.0f} m²)"
    )
    assert not np.any(data[~np.isnan(data)] == 0.0), "zero-depth cells: nodata leaked"
Clip Last: a Kernel Needs Neighbours the Clip Has Already Removed Two orderings of the same two operations on a coastal grid. In the clip-then-slope ordering, the raster is masked to the water first, so when the three-by-three slope kernel is applied the outermost ring of water cells has nodata on its landward side; the computed slope there is derived from fewer samples and is wrong, and the diagram marks that ring as a band of incorrect values following the entire coastline. In the slope-then-clip ordering, the slope is computed on the full grid where every water cell still has land neighbours available, and the mask is applied afterwards; the same ring now holds correct values. A note adds that the affected ring is the coastal margin, which is usually the part of the product the analysis exists to examine. clip → slope shaded ring: kernel saw nodata on the landward side → a band of wrong slope along the whole coast slope → clip same ring, computed while land neighbours still existed → correct, then masked away or kept as chosen The affected ring is the coastal margin — which is the part of the product most analyses exist to look at.

Edge cases and gotchas

  • A mask polygon that extends beyond the raster. from_bounds returns a window partly outside the dataset; round_lengths does not clip it, and the read returns a smaller array than the window implies, silently misaligning the mask. Intersect the window with the dataset window before reading.
  • Multipart masks with many small parts. geometry_mask rasterizes every part; a mask of 40,000 islets is slow and, worse, allocates per part. Dissolve to a single geometry first with union_all, which is one call and typically an order of magnitude faster.
  • Float rasters and integer nodata. A float32 bathymetric grid whose profile declares nodata=-9999 is fine; the same grid written as int16 with nodata=-9999 silently clips real depths below −9999 cm. Check that the nodata sentinel is outside the data’s actual range, not merely outside its plausible range.

Pipeline placement: clip last, and only for publication

The instinct is to clip early — trim the raster to the study area, then process the smaller grid. For coastal rasters that is usually the wrong order, for two reasons.

Neighbourhood operations need a halo. Slope, curvature, smoothing and gap-filling all read a window around each cell, so a grid clipped to the analysis boundary produces wrong values in the outermost ring of cells: the kernel finds nodata where it needed neighbours. Clipping first and then computing terrain derivatives puts a rim of subtly wrong values around every product, and the rim is exactly the coastal margin the analysis is usually about.

Clipping is also lossy in a way that is not recoverable. Once land cells are nodata, a later decision to move the shoreline — a new datum, a re-survey, a corrected tide model — cannot be applied to the clipped product; it needs the unclipped one. Since the shoreline is a modelled surface that moves whenever the model improves, a pipeline that clips early has baked a provisional decision into its archive.

The working order is: process on the full grid with a generous margin, clip as the final publication step, and keep the unclipped product as the archival copy. Storage is cheap relative to a reprocess, and the clipped view can be regenerated from the archive in seconds whenever the shoreline changes.

The one legitimate early clip is a coarse one: trimming a continental mosaic to a generous bounding box around the survey before processing, with a margin several kernel widths wider than any operation that will run. That is a read-volume optimisation rather than a masking step, and it uses a rectangle rather than the coastline, so none of the above applies to it.

Up: Vector Shoreline & Coastal Geometry Processing