Coastal & Marine Spatial Analysis Pipelines

Applying Gaussian Filters to Marine DEMs

Gaussian smoothing of a marine Digital Elevation Model (DEM) looks like a one-line scipy.ndimage.gaussian_filter call until the grid contains acoustic voids. The moment a single NaN enters the kernel footprint, the convolution spreads that NaN across every neighbouring cell, and a survey with 0.05% void coverage emerges with holes an order of magnitude larger. A second, quieter failure appears at scale: filtering a regional grid tile-by-tile without overlap stitches in visible seams along every chunk boundary. This page resolves both within the surface smoothing pipeline, the post-interpolation regularization stage of the parent workflow, with a NaN-safe weighted convolution, chunked rasterio I/O, and a deterministic validation gate suitable for automated pipelines.

Why the NaN Bleed Happens

scipy.ndimage.gaussian_filter is a separable linear convolution: each output cell is the weighted sum of its neighbours inside the kernel radius. Multiplication and addition involving IEEE-754 NaN are themselves NaN, so if even one void falls within the radius int(ceil(3 × σ)) of an output cell, that cell becomes NaN. The effect compounds in two dimensions and propagates outward by the full kernel radius on every pass. The snippet below reproduces it on a trivial grid:

import numpy as np
from scipy.ndimage import gaussian_filter

dem = np.full((5, 5), 10.0, dtype=np.float64)
dem[2, 2] = np.nan  # a single acoustic void

bled = gaussian_filter(dem, sigma=1.0, mode="nearest")
print(int(np.isnan(bled).sum()))  # -> 25: the whole array is now NaN

A single void destroys the entire array. On real multibeam grids the voids come from beam dropouts, water-column rejection, and the unfilled cells left behind by DEM interpolation — so the smoothing stage must treat NaN as a mask to be respected, not a value to be averaged.

Step-by-Step Fix

1. Apply a weighted NaN-safe convolution

Convolve a zero-filled copy of the data and a binary validity mask with the same Gaussian kernel, then divide. The division recovers the weighted mean of only the valid neighbours, and any output cell whose neighbourhood was entirely void stays NaN. This is the Pascucci/normalized-convolution approach and it preserves sharp bathymetric breaks such as dredged channels and pipeline corridors:

import logging
import numpy as np
from scipy.ndimage import gaussian_filter

logger = logging.getLogger(__name__)


def apply_nan_safe_gaussian(dem_array: np.ndarray, sigma: float) -> np.ndarray:
    """Gaussian-smooth a marine DEM without bleeding NaN voids.

    Convolves zero-filled data and a binary validity mask independently with
    the same kernel, then divides to recover the weighted mean. Cells with no
    valid neighbour inside the kernel footprint remain NaN.

    Args:
        dem_array: 2-D float array of depths, NaN at acoustic voids.
        sigma: Gaussian standard deviation in pixels (1.5-3.0 for bathymetry).

    Returns:
        Smoothed float32 array with original void positions preserved.

    Raises:
        ValueError: if sigma is not strictly positive.
    """
    if sigma <= 0:
        raise ValueError("sigma must be > 0; use 1.5-3.0 for bathymetric grids")

    valid_mask = ~np.isnan(dem_array)
    data_filled = np.where(valid_mask, dem_array, 0.0)
    weight_mask = valid_mask.astype(np.float64)

    # Convolve data and weights with the identical kernel; mode='nearest'
    # avoids the slope discontinuity that mode='constant' injects at edges.
    smoothed_data = gaussian_filter(data_filled, sigma=sigma, mode="nearest")
    smoothed_weights = gaussian_filter(weight_mask, sigma=sigma, mode="nearest")

    with np.errstate(invalid="ignore"):
        result = np.where(smoothed_weights > 0,
                          smoothed_data / smoothed_weights,
                          np.nan)

    logger.info("smoothed %d cells, %d remain void",
                int(valid_mask.sum()), int(np.isnan(result).sum()))
    return result.astype(np.float32)

Always pass mode='reflect' or mode='nearest' explicitly. The default mode='constant' pads with zeros, injecting a depth cliff at the grid edge that corrupts nearshore wave-transformation and tidal-prism calculations. For sigma tuning across sonar frequencies and the median-filter alternative for spike removal, see the parent surface smoothing reference.

2. Process region-scale grids in overlapping chunks

A 100k × 100k 1 m grid will not fit in RAM, and naive tiling produces seams because each tile’s edge cells are smoothed without their true neighbours. Iterate rasterio block windows, expand each window by a kernel-radius overlap before reading, smooth the expanded block, then trim the overlap back off before writing. The kernel never reads past a boundary without context, so no seam appears:

import logging
import numpy as np
import rasterio
from rasterio.windows import Window

logger = logging.getLogger(__name__)


def chunked_gaussian_filter(input_path: str, output_path: str,
                            sigma: float) -> None:
    """NaN-safe Gaussian filter over a region-scale DEM via windowed I/O.

    Each block is read with a kernel-radius overlap on every side so the
    convolution sees real neighbours across tile boundaries; the overlap is
    trimmed before the block is written back.

    Raises:
        ValueError: if the source carries no CRS (drift would be undetectable).
    """
    overlap = int(np.ceil(3 * sigma)) * 2

    with rasterio.open(input_path) as src:
        if src.crs is None:
            raise ValueError(f"{input_path} has no CRS; refusing to filter")

        meta = src.meta.copy()
        meta.update(dtype="float32", compress="zstd", nodata=np.nan,
                    tiled=True, blockxsize=256, blockysize=256)

        with rasterio.open(output_path, "w", **meta) as dst:
            for _, window in src.block_windows(1):
                col_off = max(0, window.col_off - overlap)
                row_off = max(0, window.row_off - overlap)
                width = min(src.width - col_off, window.width + 2 * overlap)
                height = min(src.height - row_off, window.height + 2 * overlap)

                expanded = Window(col_off, row_off, width, height)
                data = src.read(1, window=expanded).astype(np.float64)

                smoothed = apply_nan_safe_gaussian(data, sigma)

                # Trim the overlap back to the original (un-expanded) window.
                trim_x = window.col_off - col_off
                trim_y = window.row_off - row_off
                trimmed = smoothed[trim_y:trim_y + window.height,
                                   trim_x:trim_x + window.width]
                dst.write(trimmed.astype(np.float32), 1, window=window)

    logger.info("wrote %s with %d-px overlap", output_path, overlap)

The windowed read/write cycle inherits the source transform and CRS straight from src.meta, so georeferencing survives intact — never strip geospatial metadata on intermediate writes, and confirm horizontal and vertical referencing against the CRS alignment workflow before filtering. Convert CF-compliant NetCDF bathymetry to a Float32 GeoTIFF with rioxarray first; the trade-offs are covered in NetCDF vs GeoTIFF for marine data. Authoritative vertical-referencing requirements are in the IHO S-44 Standards for Hydrographic Surveys.

Naive Tiling vs Overlapping Windowed Filtering Left: two adjacent DEM tiles smoothed independently. Each tile's edge cells lack their true neighbours, so a discontinuity appears along the shared boundary as a seam. Right: the same write window is read with a kernel-radius overlap (a halo) on every side, the expanded block is Gaussian-smoothed so the kernel always sees real neighbours, then the halo is trimmed off before the original window is written, producing a seamless mosaic. Naive tiling edge cells smoothed without neighbours tile A tile B seam at boundary slope discontinuity Overlap → smooth → trim kernel always sees real neighbours expanded read window (halo = kernel radius) write window trim halo before write no seam continuous mosaic read +halo smooth trim+write

Verification and Acceptance Test

Gate the output deterministically: compute slope statistics and void coverage, and fail the build if smoothing flattened real geomorphology or if the void footprint grew. A NaN-safe filter should never increase void coverage beyond the original survey gaps.

import numpy as np
import rasterio


def validate_filtered_dem(dem_path: str,
                          max_slope_deg: float = 45.0) -> dict:
    """Compute slope and void statistics for a CI/CD gate.

    Returns a dict with slope stats plus boolean pass flags the pipeline
    can assert on directly.
    """
    with rasterio.open(dem_path) as src:
        data = src.read(1).astype(np.float64)
        res_y, res_x = src.res  # pixel size in CRS units

    valid = ~np.isnan(data)
    gy, gx = np.gradient(np.where(valid, data, np.nan), res_y, res_x)
    slope_deg = np.degrees(np.arctan(np.hypot(gx, gy)))
    valid_slope = slope_deg[valid]

    stats = {
        "mean_slope_deg": float(np.nanmean(valid_slope)),
        "max_slope_deg": float(np.nanmax(valid_slope)),
        "slope_violations": int(np.sum(valid_slope > max_slope_deg)),
        "nan_coverage_pct": float(100.0 * (~valid).sum() / valid.size),
    }
    return {
        "stats": stats,
        "passes_slope_gate": stats["slope_violations"] < 10,
        "passes_void_gate": stats["nan_coverage_pct"] < 0.1,
    }

Wire it into the pipeline as a hard assertion so a regression cannot ship:

python -c "from validate import validate_filtered_dem as v; \
import sys; r = v('smoothed.tif'); \
sys.exit(0 if r['passes_slope_gate'] and r['passes_void_gate'] else 1)"

Back this with pytest regression tests over known benchmark tiles, and difference the smoothed surface against the raw point cloud to confirm noise was suppressed while channel thalwegs, reef crests, and anthropogenic structures survived. Spike removal that precedes this stage is handled by multibeam point-cloud filtering and bathymetric artifact and noise removal.

Edge Cases and Gotchas

  • mode='constant' silently injects edge cliffs. The scipy default pads with cval=0.0, so a 12 m-deep survey gains a phantom step to 0 m at the grid border that flows straight into tidal-prism and wave models. Always set mode='nearest' or 'reflect', and re-check after any vertical-datum shift from the tidal datum transformation stage, which changes the sign and magnitude of edge values.
  • Anisotropic pixels skew the effective sigma. When src.res differs in x and y (common after reprojection), a scalar sigma smooths more aggressively along the finer axis. Pass a per-axis tuple sigma=(sy, sx) scaled by the pixel sizes, or resample to a square grid first.
  • NetCDF3 vs NetCDF4 nodata handling. Legacy NetCDF3 files often encode voids as a sentinel _FillValue (e.g. -9999) rather than NaN; convolving them as real depths poisons the filter far worse than NaN would. Mask the fill value to NaN on read before calling apply_nan_safe_gaussian, and prefer float32 over float64 outputs to halve the COG footprint without measurable accuracy loss at survey precision.

Up: Surface Smoothing Algorithms in Python