Surface Smoothing Algorithms in Python

Surface smoothing is the post-interpolation regularization stage within the Bathymetric Processing & Terrain Modeling pipeline: it suppresses high-frequency sonar acquisition noise and interpolation ringing without corrupting the validated depth values that precede it. It sits immediately after gridding and before terrain-derivative extraction, so any error it introduces propagates into every slope, aspect, and Bathymetric Position Index layer computed downstream. This page specifies the algorithm-selection logic, a memory-bounded Python implementation built on dask, scipy, and rasterio, the mandatory validation gates that confirm the filter behaved deterministically, and the downstream handoff contract for hydrodynamic modeling and habitat classification workflows. The objective is a smoothed, georeferenced surface whose every depth excursion stays inside the survey’s stated vertical uncertainty.

Reference Configuration #

Algorithm Selection for Marine Terrain #

Marine bathymetry exhibits spatial characteristics that constrain kernel choice: steep slope breaks at shelf edges and canyon walls, low-relief abyssal plains, and high-frequency acquisition artifacts from multipath returns, water column backscatter, and vessel motion residuals. No single kernel handles all three simultaneously without trade-offs.

Isotropic Gaussian convolution applies a weighted average in a configurable radius, attenuating high-frequency noise across the full spectral range. It is computationally inexpensive, integrates cleanly with dask.array.map_overlap, and is appropriate when the DEM will feed hydrodynamic models that require smooth gradients. The trade-off is progressive blunting of sharp morphological edges — shelf breaks, fault scarps, and anthropogenic features such as cable routes narrow in amplitude with increasing sigma.

Median filtering suppresses salt-and-pepper spikes while preserving edges more faithfully than Gaussian convolution, at the cost of a larger computational kernel and slightly different boundary behaviour. It is the preferred choice when the upstream DEM interpolation used IDW or nearest-neighbour methods that generate discontinuous void-fill artifacts.

Anisotropic diffusion and total variation minimization preserve edges even more aggressively but introduce non-linear computation that is difficult to parallelize efficiently across dask chunks. Reserve them for single-survey QC workflows, not regional-scale tiling pipelines.

The diagram below maps algorithm choice against survey resolution and downstream application:

Smoothing must follow point cloud filtering for multibeam sonar, not precede it. Applying a convolution kernel to unfiltered soundings propagates systematic errors — multipath spikes, refraction residuals, vessel heave artifacts — into the regularized surface in a way that subsequent processing cannot recover. Once verified soundings are gridded, smoothing acts on the interpolated surface as a spatial low-pass filter.

Memory-Constrained Python Implementation #

The implementation below handles Cloud Optimized GeoTIFFs and NetCDF bathymetry through xarray’s lazy backend, validates CRS metadata before any computation, and uses dask.array.map_overlap to manage convolution kernel overlap at chunk boundaries. The output is a tiled COG with embedded processing metadata.

import logging
from typing import Literal

import dask.array as da
import numpy as np
import rasterio
import xarray as xr
from rasterio.crs import CRS
from rasterio.transform import from_bounds
from scipy.ndimage import gaussian_filter, median_filter

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
logger = logging.getLogger(__name__)

SmoothMethod = Literal["gaussian", "median"]


def _smooth_block(
    block: np.ndarray,
    method: SmoothMethod,
    sigma: float,
    kernel_size: int,
) -> np.ndarray:
    """
    Apply convolution to one dask chunk, preserving NaN acoustic voids.

    Zero-filling is used for convolution stability: convolution treats
    NaN as 0, and the mask is restored after filtering so voids are
    never interpolated across.
    """
    mask = np.isnan(block)
    if np.all(mask):
        return block  # All-void chunk — skip computation

    filled = np.where(mask, 0.0, block)

    if method == "gaussian":
        smoothed = gaussian_filter(filled, sigma=sigma, mode="nearest")
    elif method == "median":
        smoothed = median_filter(filled, size=kernel_size, mode="nearest")
    else:
        raise ValueError(f"Unsupported smoothing method '{method}'. Use 'gaussian' or 'median'.")

    return np.where(mask, np.nan, smoothed)


def smooth_marine_dem(
    input_path: str,
    output_path: str,
    method: SmoothMethod = "gaussian",
    sigma: float = 1.5,
    kernel_size: int = 3,
    chunk_px: int = 2048,
) -> None:
    """
    Production-grade out-of-core smoothing for marine digital elevation models.

    Reads Cloud Optimized GeoTIFFs or NetCDF bathymetry via xarray,
    validates projected CRS metadata, applies overlap-aware chunked
    convolution, and writes a tiled COG with embedded processing provenance.

    Args:
        input_path:  Path to input GeoTIFF or NetCDF file.
        output_path: Destination path for the smoothed COG.
        method:      'gaussian' for spectral noise suppression; 'median' for
                     spike and void-fill artifact removal.
        sigma:       Gaussian sigma in pixels (ignored for median).
        kernel_size: Median filter window in pixels (ignored for gaussian).
        chunk_px:    Dask chunk size in pixels per side. 2048 balances
                     RAM usage and dask task-graph overhead on 16 GB systems.
    """
    logger.info("Opening %s with %d×%d px chunks", input_path, chunk_px, chunk_px)
    ds: xr.DataArray = xr.open_dataarray(
        input_path, chunks={"y": chunk_px, "x": chunk_px}
    )

    # --- CRS validation -------------------------------------------------------
    # CRS absence is a silent corruption vector: downstream hydrodynamic models
    # will misplace the DEM without raising an error. Fail loudly here.
    crs_str: str | None = ds.attrs.get("crs") or (
        ds.coords.get("spatial_ref", xr.DataArray()).attrs.get("crs_wkt")
    )
    if not crs_str:
        raise RuntimeError(
            "CRS metadata is absent. Assign a valid EPSG code or WKT string "
            "to the 'crs' attribute before smoothing."
        )
    crs = CRS.from_user_input(crs_str)
    if crs.is_geographic:
        raise RuntimeError(
            f"CRS {crs.to_epsg()} is geographic (degrees). Re-project to a "
            "metric projected CRS (e.g. UTM) before applying pixel-based kernels."
        )
    logger.info("CRS validated: %s", crs.to_string())

    # --- Overlap depth ---------------------------------------------------------
    # The overlap depth ensures the convolution kernel has access to neighbour
    # values across chunk boundaries. Insufficient depth introduces visible
    # seam lines in the output DEM.
    depth = int(np.ceil(sigma * 3)) if method == "gaussian" else kernel_size // 2
    logger.info("Using %s filter | depth=%d | sigma=%.2f | kernel=%d",
                method, depth, sigma, kernel_size)

    arr: da.Array = ds.data

    filtered: da.Array = da.map_overlap(
        _smooth_block,
        arr,
        method=method,
        sigma=sigma,
        kernel_size=kernel_size,
        depth=depth,
        boundary="nearest",
        dtype=arr.dtype,
    )

    logger.info("Executing smoothing graph (this may take several minutes for large surveys)...")
    result: np.ndarray = filtered.compute()

    # --- COG output ------------------------------------------------------------
    x_min = float(ds.x.min())
    x_max = float(ds.x.max())
    y_min = float(ds.y.min())
    y_max = float(ds.y.max())

    out_arr: np.ndarray = result if result.ndim == 2 else result[0]

    profile = {
        "driver": "GTiff",
        "dtype": "float32",
        "count": 1,
        "height": out_arr.shape[0],
        "width": out_arr.shape[1],
        "crs": crs,
        "transform": from_bounds(x_min, y_min, x_max, y_max,
                                  out_arr.shape[1], out_arr.shape[0]),
        "compress": "deflate",
        "predictor": 3,       # floating-point predictor improves ratio
        "tiled": True,
        "blockxsize": 256,
        "blockysize": 256,
        "nodata": float("nan"),
    }

    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(out_arr.astype(np.float32), 1)
        dst.update_tags(
            smoothing_method=method,
            sigma=str(sigma),
            kernel_size=str(kernel_size),
            chunk_px=str(chunk_px),
            processing_pipeline="coastal_marine_spatial_v2",
        )

    logger.info("Smoothed DEM written to %s", output_path)

NaN-Aware Convolution at Void Boundaries #

The _smooth_block function above zero-fills acoustic voids before convolving and restores the original NaN mask afterward. This is stable and fast, but it has a known edge artifact: cells adjacent to a void are pulled toward the zero fill value, biasing them toward unrealistically shallow depths along survey gaps and shoreline cutlines. For surveys with extensive coverage holes — interleaved swaths, turn-arounds, or masked wreck exclusion zones — replace the zero-fill kernel with a normalized convolution that ignores void cells entirely.

The normalized form convolves both the data (with voids set to zero) and a validity mask, then divides one by the other so that each output cell is the kernel-weighted mean of only its valid neighbours:

import numpy as np
from scipy.ndimage import gaussian_filter


def _smooth_block_nan_aware(
    block: np.ndarray,
    sigma: float,
) -> np.ndarray:
    """
    NaN-aware Gaussian smoothing via normalized convolution.

    Each output cell is the Gaussian-weighted mean of its *valid*
    neighbours only, so void edges are not biased toward the zero
    fill value. Cells that were void in the input remain void.
    """
    mask = np.isnan(block)
    if np.all(mask):
        return block

    valid = (~mask).astype(np.float32)
    filled = np.where(mask, 0.0, block).astype(np.float32)

    # Convolve data and validity weights with the identical kernel.
    num = gaussian_filter(filled, sigma=sigma, mode="constant", cval=0.0)
    den = gaussian_filter(valid, sigma=sigma, mode="constant", cval=0.0)

    # Avoid divide-by-zero where no valid neighbour exists in range.
    with np.errstate(invalid="ignore", divide="ignore"):
        smoothed = np.where(den > 1e-6, num / den, np.nan)

    # Restore original voids — never interpolate across a true gap.
    return np.where(mask, np.nan, smoothed)

Swap _smooth_block for _smooth_block_nan_aware in the map_overlap call when the input coverage mask is fragmented. The normalized variant roughly doubles per-chunk cost (two convolutions instead of one) but eliminates the shallow-bias halo around every void, which otherwise survives into the slope and curvature derivatives used for benthic habitat classification. astropy.convolution.convolve implements the same behaviour with a configurable nan_treatment if you prefer a maintained dependency over a hand-rolled kernel.

Validation Gates and Quality Control #

After smoothing, three mandatory checks confirm the filter behaved deterministically and did not corrupt depth values.

1. Variance ratio check #

Smoothing must reduce variance; if variance increases, the filter has introduced ringing or the overlap depth was too small. Compute the ratio before writing to disk:

import numpy as np

def check_variance_ratio(
    original: np.ndarray,
    smoothed: np.ndarray,
    max_ratio: float = 0.95,
) -> None:
    """Variance must decrease after smoothing. Ratios above max_ratio indicate instability."""
    orig_var = float(np.nanvar(original))
    smooth_var = float(np.nanvar(smoothed))
    ratio = smooth_var / orig_var if orig_var > 0 else 1.0
    if ratio > max_ratio:
        raise ValueError(
            f"Smoothing variance ratio {ratio:.4f} exceeds threshold {max_ratio}. "
            "Increase overlap depth or reduce sigma."
        )

2. Maximum depth excursion check #

Smoothing must not push any grid cell beyond the survey uncertainty bounds. For most multibeam systems this is ± 0.5 % of depth (IHO S-44 Order 1a) but use the actual survey specification:

def check_depth_excursion(
    original: np.ndarray,
    smoothed: np.ndarray,
    tolerance_m: float = 0.5,
) -> None:
    """Maximum cell-wise depth change must stay within survey uncertainty."""
    diff = np.abs(smoothed - original)
    max_diff = float(np.nanmax(diff))
    if max_diff > tolerance_m:
        raise ValueError(
            f"Maximum depth excursion {max_diff:.3f} m exceeds {tolerance_m} m tolerance. "
            "Reduce sigma or switch to median filtering."
        )

3. CRS round-trip check #

After writing, re-open the COG and confirm the CRS and pixel resolution survived the rasterio write:

import rasterio
from rasterio.crs import CRS

def verify_cog_crs(output_path: str, expected_epsg: int) -> None:
    with rasterio.open(output_path) as src:
        written_epsg = src.crs.to_epsg()
        if written_epsg != expected_epsg:
            raise RuntimeError(
                f"CRS mismatch: expected EPSG:{expected_epsg}, "
                f"found EPSG:{written_epsg} in output COG."
            )

Run all three gates inside the smoothing function before the COG is written, and fail the run rather than emitting a surface that breaches them. The variance and depth-excursion thresholds overlap with the acceptance criteria used during bathymetric artifact removal, so the same uncertainty budget governs both stages and a value rejected there will not silently reappear here. For parameter tuning specific to high-resolution surveys, the dedicated guide on applying Gaussian filters to marine DEMs covers spectral response analysis and sigma selection by resolution class.

Common Failure Modes and Diagnosis #

Seam lines in the output DEM #

Symptom: Visible horizontal or vertical banding in the smoothed raster, most obvious when rendered with a hillshade.

Root cause: The depth parameter in map_overlap is smaller than the effective kernel radius. For a Gaussian with sigma=2.0, the kernel spans roughly ceil(2.0 × 3) = 6 pixels; a depth of 2 leaves four pixels of the kernel without valid neighbours at each chunk boundary.

Remediation: Use depth = int(np.ceil(sigma * 3)) for Gaussian; for median use depth = kernel_size // 2 + 1. Rebuild the dask graph and recompute.

Geographic CRS passed to a pixel-unit kernel #

Symptom: RuntimeError: CRS ... is geographic (degrees).

Root cause: The input raster’s CRS is stored in geographic degrees (e.g. EPSG:4326). Sigma values in pixel units correspond to different metric distances depending on latitude, making the filter spatially non-uniform and geomorphically misleading.

Remediation: Re-project to a metric CRS before smoothing. For coastal North America, EPSG:32617 (UTM Zone 17N) or a regional State Plane / NAD83 CRS is appropriate. The CRS alignment pipeline page covers datum selection and re-projection patterns in detail.

NaN voids expanding across acoustic data #

Symptom: After smoothing, void regions (no-data) are larger than in the input. Bathymetric coverage has shrunk.

Root cause: The zero-fill approach in _smooth_block pulls the smoothed value at void-edge cells towards zero (no-data fill value), and the downstream np.where(mask, np.nan, smoothed) call only masks cells that were NaN in the original chunk — but the convolution has already blended edge values with the zero fill, pulling cells near voids toward unrealistic shallow values rather than expanding the NaN mask.

Remediation: Switch to the normalized _smooth_block_nan_aware kernel shown in the NaN-aware convolution subsection above, which divides the smoothed data by a smoothed validity mask so void edges are never blended with the zero fill. The astropy.convolution.convolve function handles the same case natively. Alternatively, erode the no-data boundary by one cell before smoothing using scipy.ndimage.binary_erosion on the mask array.

Out-of-memory crash during compute() #

Symptom: MemoryError or Linux OOM-killer termination during filtered.compute().

Root cause: chunk_px=2048 produces 2048 × 2048 × 4 bytes ≈ 16 MB per chunk at float32, but map_overlap materialises input chunk plus overlap on both sides — effectively a (2048 + 2 × depth)² region. For sigma=2.5 and depth=8, this is 2064 × 2064 × 4 ≈ 17 MB; that is safe. The problem arises when multiple workers materialise their chunks simultaneously. Use dask.config.set({"num_workers": 2}) or reduce chunk_px to 1024 on memory-limited systems.

Pipeline Integration and Downstream Handoff #

Surface smoothing is the final terrain regularization stage before outputs branch into two downstream paths:

1. Hydrodynamic model ingestion: Smoothed DEMs are consumed by models such as ADCIRC, Delft3D, or FVCOM. These systems import raster bathymetry via NetCDF or ASCII grid formats. Export the smoothed COG to NetCDF using xarray.Dataset.to_netcdf() with CF-1.8 conventions; include the standard_name: sea_floor_depth_below_sea_surface attribute and a grid_mapping variable referencing the CRS.

2. Habitat classification and benthic mapping: Terrain derivatives — slope, aspect, Bathymetric Position Index (BPI) — are computed from the smoothed surface. These feed removing bathymetric artifacts and noise validation checks and classification workflows.

The metadata manifest accompanying each smoothed raster must include:

Log the manifest to a sidecar .json file alongside the COG. Downstream teams use this manifest to reproduce the exact smoothing parameters during regulatory audits and multi-agency data exchanges without re-processing the full survey.

For authoritative convolution mathematics and mode parameter behaviour, consult the SciPy ndimage.gaussian_filter documentation. Output rasters must conform to Cloud Optimized GeoTIFF specifications to enable parallelized streaming in web-GIS environments and cloud-native analysis pipelines.

---

Related

• Applying Gaussian Filters to Marine DEMs — sigma selection, spectral response analysis, and resolution-class tuning
• DEM Interpolation Techniques for Seafloor Mapping — the upstream interpolation stage whose output this workflow regularizes
• Point Cloud Filtering for Multibeam Sonar — sounding-level cleaning that must precede gridding and smoothing
• Removing Bathymetric Artifacts and Noise — artifact classification and automated spike removal on the pre-smoothing surface
• Up: Bathymetric Processing & Terrain Modeling