Coastal & Marine Spatial Analysis Pipelines

Choosing a Grid Cell Size for Multibeam Data

The cell size of a bathymetric grid is chosen in about four seconds, usually by picking a round number, and it determines more about the product than any later processing decision. Too fine and most cells contain a single sounding, so the grid is a scatter plot with interpolation between the dots and every statistic computed per cell is meaningless. Too coarse and the survey’s real resolution — paid for in vessel days — is averaged away and cannot be recovered. This page derives the number instead of picking it, within the DEM interpolation workflow.

Root cause: density is not uniform, and neither is the beam footprint

Two things vary across a multibeam survey and both bear directly on the cell size.

Sounding density varies with depth and across the swath. A multibeam sonar emits a fixed number of beams per ping regardless of depth, so the swath widens with depth and the same beams cover more ground: density falls roughly with the square of depth. Across a single swath, outer beams are more widely spaced than nadir beams for the same reason.

Beam footprint grows with depth and angle. Each beam illuminates an ellipse on the seabed whose size is set by the beam width and the slant range. A 1° × 1° system at 50 m nadir has a footprint under a metre; the same system at 500 m has one nearly nine metres across. A grid finer than the footprint is asking the data for a resolution the physics did not deliver.

import numpy as np

def beam_footprint_m(depth_m: float, beam_width_deg: float, angle_deg: float = 0.0) -> float:
    """Along-track footprint diameter of one beam, at nadir or off to the side."""
    theta = np.radians(beam_width_deg)
    slant = depth_m / np.cos(np.radians(angle_deg))
    return 2.0 * slant * np.tan(theta / 2.0) / np.cos(np.radians(angle_deg))

for d in (20, 50, 200, 500):
    print(f"{d:>4} m: nadir {beam_footprint_m(d, 1.0):5.2f} m  "
          f"60 deg {beam_footprint_m(d, 1.0, 60.0):6.2f} m")
#   20 m: nadir  0.35 m  60 deg   1.40 m
#   50 m: nadir  0.87 m  60 deg   3.49 m
#  200 m: nadir  3.49 m  60 deg  13.96 m
#  500 m: nadir  8.73 m  60 deg  34.91 m

A survey spanning 20 m to 500 m has a physical resolution limit varying by a factor of twenty-five across its own extent. No single cell size is correct for all of it, which is the fact that a variable-resolution grid exists to handle and a fixed grid has to compromise around.

Two Bounds on Cell Size, Moving in Opposite Directions A chart with water depth on the horizontal axis from 20 to 500 metres. Two curves are plotted against a vertical cell-size axis. The lower bound curve, set by the nadir beam footprint, rises linearly from under a metre at 20 metres depth to nearly nine metres at 500 metres. The upper bound curve, set by the requirement that a cell contain at least a target number of soundings, also rises because density falls with depth, but faster; it starts near one metre and reaches roughly twelve metres. The band between the two curves is shaded as the range of defensible cell sizes, and it is wide in shallow water and narrow in deep water. A horizontal line at one metre marks a fixed grid choice, showing that it sits inside the band nearshore and below the footprint bound beyond about 60 metres, where it asks the data for resolution the sonar did not deliver. water depth (m) cell size (m) lower bound: beam footprint upper bound: soundings per cell defensible range a fixed 1 m grid inside the band below the footprint — resolution that is not there 20150300500 The band narrows with depth: in shallow water almost any sensible number works, in deep water very few do.

Step-by-step derivation with production code

1. Measure the density you actually have

The nominal figure from the acquisition log is an average over the whole survey and is wrong nearly everywhere. Measure it:

import logging
from dataclasses import dataclass

import numpy as np
from scipy.spatial import cKDTree

logger = logging.getLogger(__name__)


@dataclass(frozen=True)
class DensityStats:
    per_m2_p10: float          # sparse tail — the depth and swath-edge cells
    per_m2_median: float
    per_m2_p90: float


def sounding_density(xy: np.ndarray, *, sample: int = 20_000, k: int = 9) -> DensityStats:
    """Local sounding density from k-nearest-neighbour spacing, in soundings/m²."""
    if xy.ndim != 2 or xy.shape[1] != 2:
        raise ValueError("xy must be an (n, 2) array of projected coordinates")
    tree = cKDTree(xy)
    rng = np.random.default_rng(0)                       # deterministic sample
    idx = rng.choice(len(xy), size=min(sample, len(xy)), replace=False)
    dist, _ = tree.query(xy[idx], k=k)
    # k points inside a disc of radius r_k gives density k / (pi r_k^2).
    r_k = dist[:, -1]
    dens = k / (np.pi * np.maximum(r_k, 1e-6) ** 2)
    stats = DensityStats(*(float(np.percentile(dens, p)) for p in (10, 50, 90)))
    logger.info("density soundings/m²: p10 %.2f, median %.2f, p90 %.2f",
                stats.per_m2_p10, stats.per_m2_median, stats.per_m2_p90)
    return stats

Using the tenth percentile rather than the median for the cell-size decision is the choice that matters. Sizing on the median produces a grid that is well populated over half the survey and empty over the sparse parts, which are the outer beams and the deep water — the areas most likely to be the reason the survey was commissioned.

2. Turn a soundings-per-cell target into a cell size 

MIN_SOUNDINGS_PER_CELL = 5      # below this, a per-cell median is not a statistic


def cell_size_from_density(dens_per_m2: float, target: int = MIN_SOUNDINGS_PER_CELL) -> float:
    """Smallest cell that still holds `target` soundings at this density."""
    if dens_per_m2 <= 0:
        raise ValueError("density must be positive")
    return float(np.sqrt(target / dens_per_m2))

Five is a defensible default and not a law. It is the point below which a per-cell median stops being robust to a single outlier, which is the whole reason the gridder computes a median rather than a mean. A survey feeding a navigation product might want more; one feeding a regional geomorphology map can accept fewer.

3. Apply the physical and standards floors 

def choose_cell_size(
    xy: np.ndarray, depth_m: np.ndarray, *, beam_width_deg: float = 1.0,
    target: int = MIN_SOUNDINGS_PER_CELL,
) -> float:
    """Cell size satisfying density, beam footprint and a rounded reporting step."""
    dens = sounding_density(xy)
    from_density = cell_size_from_density(dens.per_m2_p10, target)
    deepest = float(np.nanpercentile(np.abs(depth_m), 95))
    from_physics = beam_footprint_m(deepest, beam_width_deg)

    chosen = max(from_density, from_physics)
    # Round up to a tidy step so adjacent surveys share a lattice.
    steps = np.array([0.25, 0.5, 1, 2, 5, 10, 25, 50])
    chosen = float(steps[np.searchsorted(steps, chosen)])
    logger.info("cell size: density wants >= %.2f m, physics wants >= %.2f m, using %.2f m",
                from_density, from_physics, chosen)
    return chosen

Taking the maximum rather than a compromise is the correct operation: both are lower bounds, and violating either produces a grid making a claim the data does not support. Rounding up to a tidy step is not cosmetic — adjacent surveys gridded at 1.0 m and 1.13 m cannot be mosaicked without resampling one of them.

4. Where depth varies widely, stop using one number

A survey from 20 m to 500 m has no single defensible cell size. The standard answer is a variable-resolution grid: nested tiles at 0.5 m inshore, 2 m mid-shelf, 8 m deep, each region gridded at a size its own data supports. BAG (Bathymetric Attributed Grid) supports this natively through refinement grids; a COG does not, and the usual COG-based approach is one file per resolution band with a coverage index.

The cost is that every consumer must handle multiple resolutions. The benefit is that no part of the product claims a resolution its data cannot support, and the shallow part is not throttled to the deep part’s limit. For a survey spanning one order of magnitude in depth, a single grid throws away roughly three quarters of the shallow resolution or fabricates most of the deep detail — there is no setting that does neither.

One Bedform, Three Cell Sizes, Two Wrong Answers Three profiles across the same sand wave, gridded at different cell sizes. At 0.25 metres, most cells contain a single sounding, so the profile is a noisy sequence of spikes in which the bedform is present but buried in per-cell noise, and the per-cell median is computed from one value. At 1 metre, cells hold five to eight soundings, the noise averages out and the sand wave appears clearly with its crest and troughs intact. At 8 metres, the cell is wider than the bedform's wavelength, so the crest and trough average together and the profile is a smooth slope with no wave in it at all. A caption records that only the middle grid can be recovered from: the fine grid can be coarsened later, but the detail lost in the coarse grid cannot be recovered without regridding from the soundings. 0.25 m — 1 sounding per cell the wave is in there, under the per-cell noise 1 m — 5 to 8 per cell crest and troughs resolved 8 m — cell wider than the wavelength dashed: what was there. solid: what was published. A fine grid can be coarsened later. A coarse one cannot be sharpened — it needs a regrid from the soundings. Which is why the soundings are the archive. Only the middle grid is defensible, and it is the one nobody would have picked by eye from the first panel.

5. Publish the count, not just the depth

Every argument above rests on soundings per cell, and that number is already computed by the gridder on its way to producing a depth. Writing it out costs one integer band and answers, per cell, the question the cell-size analysis answered only on average.

def grid_with_counts(xy, z, transform, shape):
    """Median depth and sounding count per cell, in one pass."""
    inv = ~transform
    col, row = inv * (xy[:, 0], xy[:, 1])
    r, c = np.floor(row).astype("int64"), np.floor(col).astype("int64")
    inside = (r >= 0) & (r < shape[0]) & (c >= 0) & (c < shape[1])
    flat = r[inside] * shape[1] + c[inside]

    counts = np.bincount(flat, minlength=shape[0] * shape[1]).reshape(shape)
    order = np.argsort(flat, kind="stable")
    # Median per cell via the sorted-group midpoint; empty cells stay NaN.
    depth = np.full(shape[0] * shape[1], np.nan)
    starts = np.searchsorted(flat[order], np.arange(shape[0] * shape[1]))
    ends = np.searchsorted(flat[order], np.arange(shape[0] * shape[1]), side="right")
    zs = z[inside][order]
    for cell in np.flatnonzero(counts.ravel()):
        depth[cell] = np.median(zs[starts[cell]:ends[cell]])
    logger.info("gridded %d soundings into %d populated cells",
                int(inside.sum()), int((counts > 0).sum()))
    return depth.reshape(shape), counts

A consumer who receives both bands can compute their own confidence mask, exclude thin cells from a volume calculation, and see immediately where the survey ran out of data — none of which is possible from a depth band alone. It is the same argument as the uncertainty layer in the parent section: the number that qualifies the measurement is nearly free at the point it is produced and unrecoverable afterwards.

Verification and acceptance test 

import pytest


def test_grid_is_populated_where_it_claims_to_be(counts, cell_size):
    """Most cells inside the survey footprint must hold enough soundings."""
    inside = counts > 0
    thin = (counts[inside] < MIN_SOUNDINGS_PER_CELL).mean()
    assert thin < 0.20, (
        f"{thin:.0%} of populated cells hold fewer than {MIN_SOUNDINGS_PER_CELL} "
        f"soundings at {cell_size} m — the grid is finer than the data")


def test_cell_size_is_not_below_the_beam_footprint(cell_size, depth_m):
    deepest = float(np.nanpercentile(np.abs(depth_m), 95))
    assert cell_size >= beam_footprint_m(deepest, 1.0) * 0.99, (
        "cell size is below the beam footprint at the 95th-percentile depth")

Publish the per-cell sounding count as a companion band whether or not these tests are wired in. It is one integer per cell, it is the direct evidence for both assertions, and it lets a consumer decide for themselves which parts of the grid they trust.

One Depth Range, Three Resolutions, No Compromise A plan view of a survey running from a shoreline at the left to the shelf break at the right. In the nested layout, the inshore third is tiled at half a metre, the mid-shelf third at two metres and the offshore third at eight metres, each matched to the sounding density and beam footprint at that depth. Beneath it, two single-grid alternatives are shown: a uniform eight-metre grid, which is defensible offshore and discards fifteen-fold detail inshore, and a uniform half-metre grid, which is right inshore and offshore consists mostly of cells holding no soundings at all. A caption notes that a variable-resolution product costs the consumer a resolution index and costs the producer nothing that a single grid did not already cost somewhere. nested, matched to depth 0.5 m · 20–60 m depth 2 m · 60–200 m 8 m · 200–500 m one grid at 8 m defensible offshore · discards 15× detail inshore one grid at 0.5 m right inshore · offshore is mostly empty cells with interpolation between them The nested product costs the consumer a resolution index. It costs the producer nothing a single grid did not.

Edge cases and gotchas

  • Overlapping survey lines inflate density locally. Where two lines overlap, density doubles and the derived cell size shrinks. If the overlap is a small fraction of the survey, use the tenth percentile as above and it is handled; if the whole survey is 100% overlapped by design, measure density on one line rather than on the merged cloud.
  • A gridder’s “cell size” and a COG’s are not always the same thing. Some tools take a cell size in the source CRS and others in the output CRS. On a geographic grid a “0.0001 degree” cell is 11 m east-west at the equator and 8 m at 45° — and rectangular, not square.
  • Rounding up hides a marginal decision. When the derived size is 1.04 m and the tidy step is 2 m, the grid is twice as coarse as the data supports. Log both numbers so the gap is visible; where it is large, the nested-resolution route is the honest answer.

Up: DEM Interpolation Techniques for Seafloor Mapping