Coastal & Marine Spatial Analysis Pipelines

Filtering Multibeam by Beam Angle in Python

Trimming the outer swath is the first filter in almost every multibeam pipeline, and the trim value is almost always inherited: 60°, because the last project used 60°. That number is a reasonable default for a hull-mounted system in moderate depths and it is wrong at both ends — too generous for a system with a poor sound-velocity profile, too aggressive for a well-calibrated survey where the outer beams are the only coverage over a feature. This page derives the trim from the survey’s own residuals, within the point-cloud filtering workflow.

Root cause: three separate things degrade with angle

Beam angle is a proxy for three distinct problems that happen to grow together, which is why a single trim can address all of them and why the right value differs between surveys.

Geometry. The beam footprint grows with the secant of the angle, so an outer beam averages the seabed over an area several times larger than a nadir beam and reports a single depth for it.

Sound-velocity sensitivity. Refraction error scales with the horizontal path length through the water column, which is near zero at nadir and maximal at the swath edge. An error of 2 m/s in the surface layer is invisible at nadir and produces a smile or frown of tens of centimetres at 65°, as removing multibeam refraction artifacts sets out.

Signal strength. Backscatter falls with grazing angle, so the outer beams have the worst signal-to-noise ratio and the highest bottom-detection failure rate.

Only the first is purely geometric and therefore predictable. The other two depend on the water column and the seabed on the day, which is why the trim belongs to the survey rather than to the pipeline.

import numpy as np

def beam_angle_deg(across_track_m: np.ndarray, depth_m: np.ndarray) -> np.ndarray:
    """Reconstruct beam angle from across-track offset and depth.

    Positive to starboard, negative to port. Use only when the file does not
    store the angle — a reconstruction assumes a flat seabed beneath the
    transducer and is wrong on a slope by roughly the slope angle.
    """
    d = np.abs(depth_m)
    if np.any(d <= 0):
        raise ValueError("depths must be non-zero to reconstruct an angle")
    return np.degrees(np.arctan2(across_track_m, d))

The docstring caveat is the trap. On a 10° slope the reconstructed angle is wrong by about 10°, in opposite directions on the two sides of the swath, so a symmetric trim becomes asymmetric on the seabed. Where the format stores the angle — most raw formats do, most converted XYZ files do not — use the stored value.

Read the Trim Off the Residual Curve A plot of cross-line depth residual, expressed as a standard deviation in centimetres, against absolute beam angle in degrees. The curve is flat near four centimetres from nadir out to about forty-five degrees, bends upward through a knee near fifty degrees, and rises steeply beyond, reaching about twenty centimetres at seventy degrees. A horizontal line marks the survey's accuracy budget at nine centimetres and crosses the curve at fifty-eight degrees. A vertical marker at the conventional sixty-degree trim sits slightly beyond the crossing, showing that the habitual value retains a band of data that is already outside budget. A second, fainter curve from a survey with a fresher sound-velocity cast stays inside budget out to sixty-six degrees, illustrating that the correct trim is a property of the survey rather than of the sonar. absolute beam angle (degrees) residual σ (cm) accuracy budget, 9 cm crosses at 58 deg the habitual 60 deg keeps a band already out of budget fresher SVP cast: in budget to 66 deg 0305875 Two surveys, one sonar, two correct trims. The number belongs to the water column, not to the equipment.

Step-by-step fix with production code

1. Build the residual-versus-angle curve

The measurement needs cross-line data: places where two survey lines cover the same seabed from different angles. Those crossings are the only independent check available without external control.

import logging

import numpy as np
import pandas as pd
from scipy.spatial import cKDTree

logger = logging.getLogger(__name__)


def residual_by_angle(
    df: pd.DataFrame, *, radius_m: float = 1.0, bin_deg: float = 5.0
) -> pd.DataFrame:
    """Cross-line depth residual, binned by absolute beam angle.

    `df` needs columns x, y, z, angle_deg, line_id — one row per sounding,
    already corrected for tide and sound velocity.
    """
    required = {"x", "y", "z", "angle_deg", "line_id"}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"missing columns: {sorted(missing)}")

    tree = cKDTree(df[["x", "y"]].to_numpy())
    pairs = tree.query_pairs(radius_m, output_type="ndarray")
    if pairs.size == 0:
        raise ValueError(f"no soundings within {radius_m} m of each other — "
                         "cross-line residuals need overlapping lines")

    a, b = pairs[:, 0], pairs[:, 1]
    cross = df["line_id"].to_numpy()[a] != df["line_id"].to_numpy()[b]
    a, b = a[cross], b[cross]
    logger.info("%d cross-line sounding pairs within %.1f m", len(a), radius_m)

    resid = np.abs(df["z"].to_numpy()[a] - df["z"].to_numpy()[b])
    angle = np.maximum(np.abs(df["angle_deg"].to_numpy()[a]),
                       np.abs(df["angle_deg"].to_numpy()[b]))
    bins = np.arange(0, 80 + bin_deg, bin_deg)
    out = (pd.DataFrame({"angle": angle, "resid": resid})
             .groupby(pd.cut(angle, bins), observed=True)["resid"]
             .agg(["count", "std", "median"])
             .reset_index())
    return out

Attributing each pair to the larger of the two angles is deliberate: a nadir sounding disagreeing with a 65° sounding is evidence about the 65° beam, and averaging the two angles would smear that evidence into the middle of the swath.

2. Cut where the budget is exceeded, then price the cut 

def trim_angle_for_budget(curve: pd.DataFrame, budget_m: float) -> float:
    """Largest angle whose binned residual σ still fits inside the budget."""
    ok = curve[(curve["std"] <= budget_m) & (curve["count"] >= 100)]
    if ok.empty:
        raise ValueError(
            f"no angle bin meets a {budget_m:.3f} m budget — the problem is not "
            "the outer beams, it is the sound-velocity profile or the tide model")
    return float(ok["angle"].apply(lambda iv: iv.right).max())


def coverage_cost(df: pd.DataFrame, trim_deg: float) -> dict[str, float]:
    """What the trim removes: soundings, swath width, and line overlap."""
    keep = df["angle_deg"].abs() <= trim_deg
    depth = df["z"].abs().median()
    full_w = 2.0 * depth * np.tan(np.radians(df["angle_deg"].abs().max()))
    kept_w = 2.0 * depth * np.tan(np.radians(trim_deg))
    return {
        "soundings_removed_pct": float((~keep).mean() * 100.0),
        "swath_width_lost_pct": float((1 - kept_w / full_w) * 100.0),
        "kept_swath_multiple_of_depth": float(kept_w / depth),
    }

The coverage numbers are what turn a filtering decision into a survey-planning one. Trimming from 70° to 55° removes roughly a third of the swath width, which means the line spacing chosen for the acquisition no longer delivers the planned overlap — and a survey that was planned at 70° and processed at 55° has holidays that nobody planned for.

The Trim Is a Survey-Planning Decision in Disguise Two cross-sections of adjacent survey lines. In the upper case, each line's swath extends to seventy degrees either side of nadir and the two swaths overlap by a comfortable margin, matching the line spacing the survey was planned around. In the lower case, the same two lines are trimmed at fifty-five degrees; each swath is narrower and the two no longer meet, leaving an unsurveyed strip down the middle between every pair of lines. A caption notes that the data for the strip exists but has been rejected, so the choice is between a gap and coverage that is outside the accuracy budget, and that this is a decision for the survey lead rather than for the filter. planned and processed at 70 degrees swaths overlap as planned same lines, trimmed at 55 degrees unsurveyed strip, every line pair The soundings for that strip exist and have been rejected. The choice is a gap or data outside budget — which is the survey lead's call, not the filter's.

3. Apply the trim where the angle still exists

The filter itself is trivial; where it sits in the pipeline is not. Beam angle is available in the raw format and in most intermediate formats, and it disappears the moment the cloud is exported to a plain XYZ or a generic LAS without an extra dimension carrying it. A pipeline that converts early has to reconstruct the angle from geometry, with the flat-seabed error described above.

def trim_by_angle(df: pd.DataFrame, trim_deg: float, *, keep_rejected: bool = True):
    """Classify by beam angle. Rejected soundings are marked, not deleted."""
    if not 0 < trim_deg <= 90:
        raise ValueError("trim_deg must be in (0, 90]")
    if "angle_deg" not in df.columns:
        raise ValueError(
            "no angle_deg column: the angle was lost in an earlier conversion. "
            "Reconstructing it from across-track offset assumes a flat seabed "
            "and is wrong on a slope by roughly the slope angle.")

    out = df.copy()
    out["rejected_by"] = out.get("rejected_by", pd.Series("", index=out.index))
    over = out["angle_deg"].abs() > trim_deg
    out.loc[over & (out["rejected_by"] == ""), "rejected_by"] = f"angle>{trim_deg:g}"
    logger.info("angle trim at %g deg: %d of %d soundings marked",
                trim_deg, int(over.sum()), len(out))
    return out if keep_rejected else out[out["rejected_by"] == ""].copy()

Marking rather than deleting follows the rule from the parent topic: the rejected set is the diagnostic. A rejection rate that jumps between adjacent lines is a sound-velocity problem or a sea-state problem, and it is only visible if the rejections survive to the QC report.

The rejected_by column also makes filter interactions legible. When the statistical outlier filter runs afterwards, the soundings it removes are ones the angle trim kept, so the two counts are independent and comparable. Overwriting a single boolean flag would conflate them and make the ordering question from the parent topic unanswerable after the fact.

The Rejection Rate Is a Diagnostic, Not Waste A bar chart of angle-trim rejection rate per survey line across one acquisition day. Most lines sit near four per cent, forming a steady baseline. Two adjacent lines rise to eleven and fourteen per cent. A time axis beneath marks when sound-velocity casts were taken, and the two elevated lines fall in the ninety-minute window after the last cast expired, before the next one was taken. A caption states that the rejection counter is the earliest available signal that the profile has gone stale, and that it is only available if rejected soundings are marked rather than deleted. rejected % baseline about 4% 11% and 14% cast cast cast The two bad lines sit in the window after a cast expired. The counter saw it before anyone looked at the surface.

Verification and acceptance test 

def test_trim_leaves_the_planned_overlap_intact(df, trim_deg, line_spacing_m):
    depth = float(df["z"].abs().median())
    kept_half_width = depth * np.tan(np.radians(trim_deg))
    assert 2 * kept_half_width > line_spacing_m, (
        f"at {trim_deg}° the swath is {2 * kept_half_width:.0f} m wide but the "
        f"lines are {line_spacing_m:.0f} m apart — this trim opens holidays")


def test_trim_actually_improves_the_residual(df, trim_deg):
    before = residual_by_angle(df)["std"].max()
    after = residual_by_angle(df[df["angle_deg"].abs() <= trim_deg])["std"].max()
    assert after < before * 0.8, (
        "trimming barely changed the worst residual — the error is not "
        "angle-dependent, so this filter is removing data for nothing")

The second test is the one that stops the trim becoming a ritual. If the residual is flat with angle, the outer beams are fine and the problem is elsewhere; trimming in that case removes coverage and improves nothing.

Edge cases and gotchas

  • Asymmetric trims are sometimes correct. A vessel with a heel bias, or a system with one damaged array half, produces genuinely different quality port and starboard. Measure both sides separately before assuming symmetry.
  • Angle is not stored in XYZ. Converting to XYZ before trimming throws away the field the trim needs, which is one more reason the LAS-to-XYZ conversion belongs downstream of filtering rather than upstream of it.
  • Trim before, not after, statistical filters. A density- or distribution-based outlier test assumes one population; a full swath is two, and running the test first deletes the outer beams for being outer beams. The ordering argument is set out in the parent topic.

Up: Point Cloud Filtering for Multibeam Sonar