Coastal & Marine Spatial Analysis Pipelines

Removing Bathymetric Artifacts and Noise

This workflow sits inside the Bathymetric Processing & Terrain Modeling pipeline as the deterministic preprocessing stage that runs immediately after gridding and before any terrain derivative is computed. Raw multibeam echosounder (MBES), singlebeam (SBES), and airborne LiDAR-derived bathymetry routinely contain acquisition spikes, nadir voids, cross-track striping, and edge ringing. These anomalies originate from vessel motion compensation errors, sound velocity profile miscalibration, multipath reflections, and gridding interpolation failures. If unfiltered, they propagate into every derivative terrain product — corrupting slope calculations, rugosity metrics, benthic habitat classification, and hydrodynamic boundary conditions. This page defines a production-grade, memory-constrained pipeline for isolating and suppressing surface-level anomalies prior to downstream terrain modeling. It assumes pre-aligned, CRS-validated raster or dense point-derived grids as input; rigorous point cloud filtering for multibeam sonar upstream reduces raster-level artifact frequency by 40–60%.

Bathymetric Artifact Removal Pipeline Flow Horizontal data-flow diagram. Upstream point-cloud filtering feeds a gridded depth raster into the artifact-removal stage, which runs a CRS validation gate, overlap-buffered chunked tiling, 5×5 median plus 3.5-sigma spike detection, a secondary 3×3 median pass for cross-track striping, and artifact-mask dilation, then passes through three validation gates before writing a cleaned cloud-optimized GeoTIFF and a binary artifact mask that hand off to downstream DEM interpolation. Artifact Removal Pipeline — Deterministic Preprocessing Stage Upstream Point-cloud filtering (MBES / LiDAR) Gridded depth input COG / BAG v2 / XYZ float32, NoData flagged CRS validation gate reject EPSG:4326 require metric UTM Chunked tiling 512×512 px blocks 16–32 px overlap per-tile filtering (overlap-buffered, idempotent) Spike detection 5×5 median + 3.5σ Secondary median 3×3 — striping Artifact mask composite + dilation Gate 1 mask coverage ≤ 8–12% Gate 2 slope variance ≥ 35% reduction Gate 3 cross-track continuity no terracing Cleaned COG (float32) Binary mask (uint8) Downstream DEM interpolation Each stage writes audit metadata; a failed run is diagnosed without re-executing the whole pipeline.

The end-to-end stage sequence — CRS validation, overlap-buffered chunking, statistical spike detection, secondary median suppression of striping, mask composition, and validation gating — is detailed in the sections below. Each stage is idempotent and writes intermediate audit metadata so a failed run can be diagnosed without re-executing the whole pipeline.

Reference Configuration and Specification Table

Parameter Recommended Value Notes
Input formats Cloud-Optimized GeoTIFF (COG), IHO BAG v2, tiled XYZ Must carry explicit NoData flags and spatial reference metadata
CRS requirement Metric projected (e.g., UTM — EPSG:326xx / EPSG:327xx) Geographic EPSG:4326 invalidates gradient thresholds; reject or reproject
Chunk size 512×512 px (dense grids) or 1024×1024 px (sparse grids) Align to native tile boundaries
Overlap buffer 16–32 px per tile edge Mandatory to prevent discontinuities at tile seams
Depth dtype float32 float64 doubles I/O overhead without bathymetric precision gain
Spike σ-threshold 3.5σ above local 5×5 median Increase to 4.0σ for high-energy reef surveys with genuine relief variance
Max artifact mask coverage ≤ 8–12% of grid area Higher values indicate upstream acquisition failures
Slope variance reduction ≥ 35% high-frequency reduction Without flattening geomorphic features (reef crests, channel margins)
scipy version ≥ 1.11 generic_filter behaviour stable since 1.9; pin in requirements
rasterio version ≥ 1.3 Block-window API and COG write support required
IHO S-44 depth accuracy Order 1a / Special Order as applicable Artifact thresholds must not erode IHO allowable uncertainty bands

Artifact Taxonomy and Detection Logic

Bathymetric Artifact Taxonomy Four artifact classes — acquisition spikes, nadir gaps and cross-track striping, edge ringing, and interpolation blockyness — each shown with its origin and primary detection method. Bathymetric Artifact Taxonomy & Detection Artifact Class Primary Origin Detection Method Acquisition Spikes Single-cell / few-cell depth excursions Multipath reflections, air bubbles, biological scatterers 5×5 local median + 3.5σ threshold Replace with local median estimate → automated-spike-removal page Nadir Gaps & Cross-Track Striping Systematic voids / linear depressions Poor swath overlap, beam geometry / SV miscalibration Directional gradient analysis, periodicity transforms along track-perpendicular axis Edge Ringing Gibbs-like oscillations at survey boundaries Interpolation seam discontinuities, boundary condition errors High-freq Fourier decomposition, localized Laplacian variance mapping near seams Interpolation Artifacts Blocky / terraced surfaces from sparse grids Nearest-neighbor / linear gridding on low-density point clouds Adaptive median filtering, curvature-constrained smoothing → surface-smoothing algorithms

Artifact isolation relies on multi-scale statistical thresholds and morphological operations rather than manual digitization. The four primary noise classes are:

  1. Acquisition spikes — single-cell or few-cell depth excursions exceeding local standard deviation by >3.5σ. Typically caused by multipath reflections, air bubbles, or biological scatterers. Detailed detection and handling is covered in automated spike removal in sonar datasets.
  2. Nadir gaps and cross-track striping — systematic voids or linear depressions along the vessel trackline where beam geometry yields poor coverage or swath overlap is insufficient. Detected via directional gradient analysis and periodicity transforms along the track-perpendicular axis.
  3. Edge ringing — Gibbs-like oscillations at survey boundaries or interpolation seams. Identifiable through high-frequency Fourier decomposition and localized Laplacian variance mapping near tile edges.
  4. Interpolation artifacts — blocky or terraced surfaces generated by aggressive nearest-neighbor or linear gridding on sparse point clouds. Suppressed via adaptive median filtering and curvature-constrained surface smoothing.

Pipeline Initialization and Memory Architecture

Large-scale seafloor grids routinely exceed available RAM when loaded as monolithic arrays. Pipeline stability requires explicit chunking, lazy evaluation, and strict CRS enforcement before any filtering operation executes.

CRS enforcement: All inputs must be projected to a metric coordinate system (UTM zone, EPSG:326xx for northern hemisphere). Geographic coordinates (EPSG:4326) introduce non-linear distance distortion that invalidates gradient calculations, morphological kernel sizing, and statistical distance thresholds. The validate_crs function in the implementation below must reject or raise on non-metric projections — it must not silently proceed. See CRS alignment for coastal GIS projects for datum and compound CRS decisions that precede this step.

Chunking strategy: Processing executes on 512×512 or 1024×1024 pixel blocks aligned to native tile boundaries. Overlap buffers of 16–32 pixels per edge are mandatory to prevent discontinuities during neighborhood operations. Overlap regions are discarded after filtering to maintain spatial continuity across tile seams.

Memory budget: Peak RAM must remain below 75% of available system memory. Depth values are cast to float32 at ingestion — float64 doubles I/O overhead without bathymetric precision gain. For datasets exceeding 50 GB, replace the rasterio block iterator with dask.array.from_array and map_blocks for fully lazy execution.

Memory-Constrained Python Implementation

The following implementation demonstrates a chunked, memory-safe pipeline using rasterio, numpy, and scipy.ndimage. It enforces CRS validation, applies overlap-aware filtering, and outputs a cleaned depth raster alongside a binary artifact mask.

import logging
import numpy as np
import rasterio
from rasterio.windows import Window
from scipy.ndimage import median_filter, generic_filter, binary_dilation

log = logging.getLogger(__name__)


def validate_crs(src: rasterio.io.DatasetReader) -> None:
    """Enforce metric projection. Raise on EPSG:4326 or undefined CRS."""
    if src.crs is None:
        raise ValueError(
            f"Input '{src.name}' has no CRS. A metric projected CRS is required."
        )
    if src.crs.is_geographic:
        raise ValueError(
            f"Input '{src.name}' is in geographic CRS {src.crs.to_epsg()}. "
            "Reproject to a metric CRS (e.g., UTM, EPSG:326xx) before artifact removal. "
            "Geographic coordinates invalidate gradient thresholds and kernel sizing."
        )
    log.info("CRS validated: %s (metric projected)", src.crs.to_string())


def process_chunk(
    chunk: np.ndarray,
    nodata: float | None,
    sigma_thresh: float = 3.5,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Apply statistical and morphological filtering to a single raster block.

    Returns
    -------
    cleaned : np.ndarray (float32)
        Depth values with spikes replaced by local median estimates.
    artifact_mask : np.ndarray (uint8)
        Binary mask — 1 where artifact detected or NoData.
    """
    # --- 1. Build NoData mask ---
    if nodata is not None and not np.isnan(nodata):
        nodata_mask: np.ndarray = np.isclose(chunk, nodata)
    else:
        nodata_mask = np.isnan(chunk)

    valid: np.ndarray = np.where(nodata_mask, 0.0, chunk.astype(np.float32))

    # --- 2. Spike detection: 5×5 local median + sigma threshold ---
    local_median: np.ndarray = median_filter(valid, size=5).astype(np.float32)
    local_std: np.ndarray = generic_filter(
        valid, lambda x: np.nanstd(x), size=5
    ).astype(np.float32)

    # Avoid division by zero in flat regions; floor std at 0.001 m
    local_std = np.maximum(local_std, 0.001)

    spike_mask: np.ndarray = (
        ~nodata_mask
        & (np.abs(valid - local_median) > (sigma_thresh * local_std))
    )
    log.debug("Spike pixels detected: %d", int(spike_mask.sum()))

    # Replace spikes with local median estimate
    cleaned: np.ndarray = np.where(spike_mask, local_median, valid)

    # --- 3. Secondary 3×3 median pass to suppress cross-track striping ---
    cleaned = median_filter(
        np.where(nodata_mask, 0.0, cleaned), size=3
    ).astype(np.float32)
    cleaned = np.where(nodata_mask, np.nan, cleaned)

    # --- 4. Composite artifact mask with 1-pixel dilation for edge coverage ---
    artifact_mask: np.ndarray = nodata_mask | spike_mask
    artifact_mask = binary_dilation(artifact_mask, structure=np.ones((3, 3)))

    return cleaned, artifact_mask.astype(np.uint8)


def run_pipeline(
    input_path: str,
    output_path: str,
    mask_path: str,
    chunk_size: int = 512,
    overlap: int = 16,
    sigma_thresh: float = 3.5,
) -> dict[str, float]:
    """
    Execute chunked bathymetric artifact removal pipeline.

    Parameters
    ----------
    input_path : str
        Path to input COG, BAG v2, or tiled GeoTIFF depth raster.
    output_path : str
        Destination path for cleaned depth raster (float32, NoData=NaN).
    mask_path : str
        Destination path for binary artifact mask (uint8).
    chunk_size : int
        Target tile dimension in pixels; actual chunks follow block windows.
    overlap : int
        Overlap buffer width (pixels) added to each tile edge.
    sigma_thresh : float
        Sigma multiplier for spike detection. Default 3.5.

    Returns
    -------
    dict with keys 'mask_coverage_pct' and 'tiles_processed'.
    """
    total_pixels: int = 0
    masked_pixels: int = 0
    tiles: int = 0

    with rasterio.open(input_path) as src:
        validate_crs(src)
        log.info(
            "Opening %s — size %d×%d, resolution %.3f m/px",
            input_path, src.width, src.height, src.res[0],
        )

        profile = src.profile.copy()
        profile.update(dtype="float32", nodata=float("nan"))
        mask_profile = profile.copy()
        mask_profile.update(dtype="uint8", nodata=0)

        with (
            rasterio.open(output_path, "w", **profile) as out,
            rasterio.open(mask_path, "w", **mask_profile) as mask_out,
        ):
            for _ji, window in src.block_windows(1):
                w = window

                # Expand window to include overlap buffer
                col_off = max(0, w.col_off - overlap)
                row_off = max(0, w.row_off - overlap)
                w_width = min(src.width - col_off, w.width + 2 * overlap)
                w_height = min(src.height - row_off, w.height + 2 * overlap)

                read_window = Window(col_off, row_off, w_width, w_height)
                data = src.read(1, window=read_window).astype(np.float32)

                cleaned, artifact_mask = process_chunk(
                    data, src.nodata, sigma_thresh=sigma_thresh
                )

                # Trim overlap to restore original window bounds
                trim_x = w.col_off - col_off
                trim_y = w.row_off - row_off
                sl_x = slice(trim_x, trim_x + w.width)
                sl_y = slice(trim_y, trim_y + w.height)

                out.write(cleaned[sl_y, sl_x], 1, window=w)
                mask_out.write(artifact_mask[sl_y, sl_x], 1, window=w)

                total_pixels += w.width * w.height
                masked_pixels += int(artifact_mask[sl_y, sl_x].sum())
                tiles += 1

        coverage_pct = 100.0 * masked_pixels / max(total_pixels, 1)
        log.info(
            "Pipeline complete — %d tiles, mask coverage %.2f%%", tiles, coverage_pct
        )

    return {"mask_coverage_pct": coverage_pct, "tiles_processed": tiles}

Implementation notes:

  • float32 casting at read time prevents memory bloat. np.nan propagation ensures NoData regions do not contaminate neighborhood statistics.
  • The local_std floor at 0.001 m prevents division-by-zero in perfectly flat areas (e.g., artificially filled NoData regions), which would otherwise flag every valid cell as a spike.
  • Overlap trimming guarantees seamless tile stitching. The returned dict provides audit metadata for the validation step below.
  • For datasets exceeding 50 GB, replace the rasterio block iterator with dask.array.from_array and map_blocks for fully lazy, distributed execution.

Validation Gates and Quality Control

Artifact suppression must be validated before outputs enter terrain modeling workflows. Three quantitative gates are mandatory:

Gate 1 — Mask coverage ratio. Artifact mask coverage should remain below 8–12% of total grid area. Higher percentages indicate upstream acquisition failures or overly aggressive sigma thresholds. Inspect the mask spatially before adjusting thresholds — systematic patterns (e.g., mask concentrated at trackline centres) identify origin rather than threshold miscalibration.

import rasterio
import numpy as np

def check_mask_coverage(mask_path: str, threshold_pct: float = 12.0) -> None:
    with rasterio.open(mask_path) as src:
        mask = src.read(1)
    coverage = 100.0 * mask.sum() / mask.size
    if coverage > threshold_pct:
        raise RuntimeError(
            f"Artifact mask coverage {coverage:.1f}% exceeds threshold {threshold_pct}%. "
            "Check upstream acquisition quality or raise sigma_thresh."
        )
    print(f"Mask coverage: {coverage:.2f}% — PASS")

Gate 2 — Residual slope variance. Compute slope variance before and after filtering. A successful run reduces high-frequency slope variance by ≥35% without flattening geomorphic features such as reef crests or channel margins. Compute slope using the standard 3×3 Horn gradient estimator on the cleaned grid and compare to the raw input.

Gate 3 — Cross-track continuity. Cross-track profiles must show no artificial terracing or step functions. Use directional variograms to confirm isotropic spatial continuity. Persistent directional anisotropy in the EW or NS direction after filtering signals unresolved swath striping that the secondary median pass did not suppress — revisit point cloud filtering for multibeam sonar to address the problem at the source.

Cross-link validation outputs to the automated spike removal in sonar datasets debug page for per-beam spike frequency histograms, which reveal whether spikes concentrate at outer swath angles (beam-angle artefact) or appear uniformly distributed (SV profile error).

Common Failure Modes and Diagnosis

Failure 1: ValueError: Input is in geographic CRS

Symptom: validate_crs raises immediately on pipeline start.

Root cause: Input raster stored in EPSG:4326 or another geographic CRS. This is common when files are exported from desktop GIS without reprojection, or when the CRS metadata is absent (no .prj sidecar for legacy formats).

Remediation: Reproject before passing to the pipeline:

gdalwarp -t_srs EPSG:32618 -r bilinear input_geo.tif input_utm.tif

Then re-run validate_crs to confirm before executing the full pipeline. See CRS alignment for coastal GIS projects for guidance on selecting the correct UTM zone for your survey area.

Failure 2: Mask coverage >15% on valid survey data

Symptom: check_mask_coverage raises; visual inspection shows artifact mask covering the entire nadir region.

Root cause: Sigma threshold (3.5σ) is too aggressive for surveys over high-relief terrain such as rocky reef or steep canyon walls, where genuine depth variance exceeds the flat-seafloor baseline assumption of the 5×5 median estimator.

Remediation: Raise sigma_thresh to 4.0–4.5σ for high-relief surveys. Alternatively, segment the grid by terrain roughness (e.g., using a pre-computed curvature raster) and apply zone-specific thresholds:

# Apply relaxed threshold only where curvature magnitude exceeds 0.5 m/m²
sigma_map = np.where(curvature > 0.5, 4.5, 3.5)
# Pass sigma_map into a vectorised process_chunk variant

Failure 3: Visible tile seams in cleaned output

Symptom: Horizontal or vertical linear discontinuities at multiples of 512 px in the output raster.

Root cause: Overlap buffer was too small (< 16 px) for the filter kernel sizes in use, causing the 5×5 median and secondary 3×3 pass to operate on edge-incomplete neighbourhoods.

Remediation: Increase overlap to 32 px (filter half-width × 3 is a safe minimum for combined 5×5 + 3×3 passes). Verify by checking the output with gdalinfo -stats for suspiciously aligned variance jumps, or by computing a difference raster between passes at multiple overlap settings.

Failure 4: generic_filter running extremely slowly on large tiles

Symptom: Processing a 1024×1024 tile takes >60 s; pipeline stalls on local-std computation.

Root cause: scipy.ndimage.generic_filter with a Python lambda (np.nanstd) is pure-Python and not vectorised — it calls the Python function once per pixel, making it O(N²) for large tiles.

Remediation: Replace the generic_filter call with a vectorised rolling-std approximation:

from scipy.ndimage import uniform_filter

def rolling_std(arr: np.ndarray, size: int = 5) -> np.ndarray:
    """Vectorised local std approximation via E[X²] - E[X]² identity."""
    mean = uniform_filter(arr, size=size)
    mean_sq = uniform_filter(arr ** 2, size=size)
    variance = np.clip(mean_sq - mean ** 2, 0, None)
    return np.sqrt(variance).astype(np.float32)

This runs 20–50× faster than generic_filter with a Python callable and is numerically stable for depth values in the −6000 m to 0 m range.

Pipeline Integration and Downstream Handoff

Validated outputs are staged as a paired set: the cleaned depth raster (*_cleaned.tif, float32, NoData=NaN) and the binary artifact mask (*_artifact_mask.tif, uint8). Both files are written as Cloud-Optimized GeoTIFFs with internal tiling and LZW compression:

profile.update(
    driver="GTiff",
    tiled=True,
    blockxsize=512,
    blockysize=512,
    compress="lzw",
    predictor=3,          # floating-point horizontal differencing
    bigtiff="IF_SAFER",
)

A JSON manifest accompanies each output pair:

{
  "input": "survey_20240901_utm18n.tif",
  "cleaned": "survey_20240901_cleaned.tif",
  "artifact_mask": "survey_20240901_artifact_mask.tif",
  "crs_epsg": 32618,
  "sigma_thresh": 3.5,
  "overlap_px": 16,
  "mask_coverage_pct": 4.7,
  "iho_s44_order": "1a",
  "processed_utc": "2024-09-02T08:14:33Z"
}

The iho_s44_order field records the applicable IHO S-44 standard accuracy class, ensuring artifact thresholds are traceable to recognised hydrographic survey requirements.

Cleaned grids feed directly into DEM interpolation techniques for seafloor mapping, where they serve as primary inputs for TIN generation, kriging, or spline-based surface reconstruction. The artifact mask travels with the cleaned raster and is consumed by downstream interpolation routines to exclude masked cells from the conditioning data rather than treating them as valid observations. For surveys where Gaussian smoothing is required after artifact removal, outputs pass next to the surface smoothing algorithms stage.

For operational deployment, run this pipeline inside a containerised environment with pinned scipy (≥1.11) and rasterio (≥1.3) versions. Monitor memory utilisation via resource.getrusage and enforce automatic retry on chunk-level I/O failures using exponential backoff before raising.