Dask vs Ray for Bathymetric Processing at Scale

Choosing between Dask and Ray is one of the first architectural decisions in scaling a pipeline within the Bathymetric Processing & Terrain Modeling domain, and it is one of the easiest to get wrong. Both frameworks distribute Python across cores and nodes, both are pip-installable, and both will happily run a small survey on a laptop — which is exactly why teams pick one on convenience and then hit a wall when the survey grows from a single harbour to a full coastal charting campaign of tens of billions of soundings. The two engines are optimized for different shapes of work. Dask is an array and dataframe engine: it excels when the computation is a windowed operation over one large logical grid — out-of-core gridding, chunk-aligned digital elevation model (DEM) algebra, focal statistics — and it integrates natively with xarray so a NetCDF or Zarr cube behaves like a lazy NumPy array. Ray is a task and actor engine: it excels when the work is an embarrassingly parallel map over independent units — process 4,000 survey tiles, run a point-cloud classifier per tile, fan out a machine-learning inference over swaths — where each unit is a self-contained function call with its own cost profile. This page gives a reference comparison, two production implementations, the validation gates that keep either engine honest, the failure modes that bite in production, and the handoff to the compute and orchestration layers.

When the Choice Actually Matters #

For a survey small enough to fit in RAM, neither engine buys you anything a plain NumPy or scipy script does not — and the framework overhead makes them slower. The decision matters when the dataset exceeds a single machine’s memory, when wall-clock time on one core is unacceptable, or when the pipeline must run unattended on a schedule where a memory blow-up is a paged incident rather than a re-run. At that scale the question is no longer “which is faster” — both saturate a compute cluster — but “which execution model matches the data structure,” because a mismatch forces you to fight the framework on every operation.

The clean heuristic: if you can express the whole computation as algebra over one coordinate grid, reach for Dask; if you must express it as a loop over independent things, reach for Ray. A gridding step that resamples scattered soundings into a raster and then runs a focal slope calculation is grid algebra — Dask. A step that reads 4,000 .bag tiles, cleans each with a per-tile model, and writes 4,000 cleaned tiles is a map over independent units — Ray. Many real pipelines contain both, and the mature answer is often to use each where it fits rather than force one engine to do everything.

Head-to-Head Reference Comparison #

The table condenses the axes that actually decide a marine deployment. Read it as “for this axis, which engine’s default behaviour costs me less engineering.”

Two rows deserve emphasis. xarray integration is often the whole decision: if your source is a Zarr or NetCDF cube and your operation is “smooth this cube, compute a slope, write it back,” Dask lets you write ordinary array code and streams it out-of-core for free, whereas Ray forces you to hand-partition the cube and reassemble it. Memory model is the other: Ray’s plasma object store keeps large read-only NumPy arrays once and hands zero-copy views to every task that needs them, which is a real win when 200 tile-cleaning tasks all reference the same 2 GB sound-velocity lookup or model weights.

Dask: Out-of-Core Gridding and DEM Algebra #

The canonical Dask use in bathymetry is resampling scattered soundings onto a target grid and then running windowed DEM operations that never fit one machine’s RAM. The code below builds a chunk-aligned xarray.DataArray backed by dask.array, applies a memory-bounded focal smoothing with map_overlap so chunk seams do not leak edge artifacts, and writes to Zarr. It is typed, logs progress, and raises on any non-finite leakage rather than silently emitting NaNs into a downstream DEM. This is the same grid algebra that underlies surface smoothing algorithms and DEM interpolation, scaled out-of-core.

import logging
from typing import Final

import dask.array as da
import numpy as np
import xarray as xr
from scipy.ndimage import uniform_filter

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(name)s | %(message)s",
)
log = logging.getLogger("dask_gridding")

CHUNK: Final[int] = 2048          # cells per chunk edge; tune to worker RAM
SMOOTH_WIN: Final[int] = 5        # focal window in cells (odd)
OVERLAP: Final[int] = SMOOTH_WIN // 2


def _focal_mean(block: np.ndarray) -> np.ndarray:
    """Windowed mean over a single chunk plus its overlap halo.

    Operates on a plain NumPy block supplied by Dask's map_overlap. NaN
    soundings are treated as missing: we smooth only finite depths so the
    mean is not poisoned by no-data cells.
    """
    finite = np.isfinite(block)
    filled = np.where(finite, block, 0.0)
    weight = finite.astype(np.float64)
    num = uniform_filter(filled, size=SMOOTH_WIN, mode="nearest")
    den = uniform_filter(weight, size=SMOOTH_WIN, mode="nearest")
    with np.errstate(invalid="ignore", divide="ignore"):
        out = np.where(den > 0.0, num / den, np.nan)
    return out


def smooth_dem_out_of_core(
    grid: xr.DataArray,
    output_zarr: str,
) -> xr.DataArray:
    """Focal-smooth a chunk-aligned bathymetric DEM out-of-core and write Zarr.

    Args:
        grid: 2-D DataArray (dims 'y', 'x') backed by dask.array, EPSG-tagged.
        output_zarr: destination Zarr store path or URI.

    Returns:
        The lazy smoothed DataArray (already persisted to Zarr).

    Raises:
        ValueError: if the grid is not 2-D, not Dask-backed, or if smoothing
            introduces non-finite cells where the input was finite.
    """
    if grid.ndim != 2:
        raise ValueError(f"expected a 2-D DEM, got dims={grid.dims}")
    if not isinstance(grid.data, da.Array):
        raise ValueError("grid must be dask-backed; call .chunk() first")

    log.info("smoothing DEM shape=%s chunks=%s", grid.shape, grid.data.chunksize)
    smoothed = xr.apply_ufunc(
        lambda b: da.map_overlap(
            _focal_mean, b, depth=OVERLAP, boundary="none", dtype="float64"
        ),
        grid,
        dask="allowed",
        keep_attrs=True,
    )

    # Guard: smoothing must not create NaNs where input depths were valid.
    leaked = ((~np.isfinite(smoothed)) & np.isfinite(grid)).sum().compute()
    if int(leaked) > 0:
        raise ValueError(f"{int(leaked)} cells went non-finite during smoothing")

    smoothed.attrs.update(grid.attrs)
    smoothed.to_dataset(name="depth").to_zarr(output_zarr, mode="w")
    log.info("wrote smoothed DEM to %s", output_zarr)
    return smoothed


if __name__ == "__main__":
    # Example: a 40k x 40k DEM (~12 GB float64) that never fully loads.
    lazy = da.random.normal(-25.0, 3.0, size=(40_000, 40_000)).rechunk(CHUNK)
    dem = xr.DataArray(lazy, dims=("y", "x"), attrs={"crs": "EPSG:32618"})
    smooth_dem_out_of_core(dem, output_zarr="s3://survey/dem_smoothed.zarr")

Ray: Task-Parallel Tile Processing #

The canonical Ray use is a map over independent survey tiles where each tile’s cost varies — a dense reef tile takes far longer than an empty deep-water tile — and a central graph would only get in the way. The code declares a @ray.remote task that cleans one tile, submits all tiles, and drains results with ray.wait so memory stays bounded (never hold all futures’ results at once). A shared read-only reference (here a sound-velocity table) is placed in the object store once with ray.put and passed by reference, so all tasks share one zero-copy copy. Failures raise per tile and are collected rather than aborting the fleet.

import logging
from dataclasses import dataclass
from typing import Final

import numpy as np
import ray

logging.basicConfig(level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s")
log = logging.getLogger("ray_tiles")

MAX_IN_FLIGHT: Final[int] = 64     # cap concurrent tasks to bound object-store use


@dataclass(frozen=True)
class TileResult:
    tile_id: str
    kept: int
    dropped: int


@ray.remote(num_cpus=1, max_retries=2, retry_exceptions=True)
def clean_tile(tile_path: str, svp_ref: np.ndarray, k: float = 3.5) -> TileResult:
    """Clean one multibeam tile: robust vertical-outlier cull.

    Args:
        tile_path: object-storage URI of an (N, 3) XYZ tile.
        svp_ref: shared read-only sound-velocity table (zero-copy from plasma).
        k: robust threshold in scaled-MAD units.

    Raises:
        ValueError: if the tile is empty or malformed, so the driver can
            quarantine it rather than silently emit a hole in the mosaic.
    """
    xyz = np.load(tile_path)  # (N, 3) float64
    if xyz.ndim != 2 or xyz.shape[1] != 3 or xyz.shape[0] == 0:
        raise ValueError(f"malformed tile {tile_path} shape={xyz.shape}")

    z = xyz[:, 2]
    med = np.median(z)
    mad = 1.4826 * np.median(np.abs(z - med)) or 1e-6
    keep = np.abs(z - med) <= k * mad
    tile_id = tile_path.rsplit("/", 1)[-1]
    log.info("tile=%s kept=%d dropped=%d svp_n=%d", tile_id, int(keep.sum()),
             int((~keep).sum()), svp_ref.shape[0])
    return TileResult(tile_id, int(keep.sum()), int((~keep).sum()))


def run_tile_fleet(tile_paths: list[str], svp: np.ndarray) -> list[TileResult]:
    """Map clean_tile over all tiles with a bounded in-flight window.

    Raises:
        RuntimeError: if any tile task fails after retries, aggregating the
            failures so a partial mosaic is never silently accepted.
    """
    ray.init(ignore_reinit_error=True)
    svp_ref = ray.put(svp)                      # one shared copy in plasma
    pending = [clean_tile.remote(p, svp_ref) for p in tile_paths[:MAX_IN_FLIGHT]]
    queue = tile_paths[MAX_IN_FLIGHT:]
    done: list[TileResult] = []
    failures: list[str] = []

    while pending:
        ready, pending = ray.wait(pending, num_returns=1)
        try:
            done.append(ray.get(ready[0]))
        except Exception as exc:  # noqa: BLE001 - collect, do not abort fleet
            failures.append(str(exc))
            log.error("tile task failed: %s", exc)
        if queue:
            pending.append(clean_tile.remote(queue.pop(), svp_ref))

    if failures:
        raise RuntimeError(f"{len(failures)} tile(s) failed: {failures[:3]} ...")
    log.info("fleet complete: %d tiles cleaned", len(done))
    return done

Validation Gates and Quality Control #

Neither engine is trustworthy until three properties are proven, because distributed execution introduces failure modes that a single-process script cannot have.

1. Deterministic output regardless of partitioning. The gridded result must not depend on chunk size, worker count, or task ordering. Floating-point non-associativity means a naive parallel sum can differ in the last bits between two chunkings; for bathymetry that is usually within tolerance, but a logic bug (an overlap halo that is too small, a reduction that double-counts a seam) produces large, real differences. Run the pipeline at two chunkings and assert agreement within a tight tolerance.

import numpy as np
import xarray as xr


def assert_partition_invariant(
    build,               # callable(chunk: int) -> xr.DataArray
    chunk_a: int,
    chunk_b: int,
    atol: float = 1e-6,
) -> None:
    """Assert the pipeline yields the same DEM under two chunkings.

    Raises:
        AssertionError: if any cell differs by more than atol, which signals a
            seam/overlap bug rather than benign float reordering.
    """
    a = build(chunk_a).compute()
    b = build(chunk_b).compute()
    diff = np.abs(np.nan_to_num(a.values) - np.nan_to_num(b.values))
    worst = float(diff.max())
    if worst > atol:
        raise AssertionError(f"partitioning changed output: max|Δ|={worst:.3e} > {atol:.0e}")

2. Memory ceiling. Set an explicit per-worker limit and prove a single oversized unit cannot breach it. On Dask, cap worker memory and let chunk spill engage; on Ray, cap the object-store fraction and num_cpus so concurrency is bounded. The gate is operational: run the largest expected tile or chunk under the configured limit in CI and confirm no worker exceeds it.

3. Reproducibility. Pin library versions (Dask, distributed, Ray, NumPy, xarray) in the image and record them in the output manifest. A DEM produced under NumPy 1.26 and re-run under 2.0 can differ in reduction order; without a recorded environment you cannot explain a diff months later. The relative error tolerance for the invariance gate is defined as

$$\varepsilon_{\text{rel}} = \frac{\lVert \mathbf{D}_a - \mathbf{D}_b \rVert_\infty}{\max(1,\ \lVert \mathbf{D}_a \rVert_\infty)} \le 10^{-6}$$

where $\mathbf{D}_a$ and $\mathbf{D}_b$ are the DEMs produced under two partitionings. A breach above $10^{-6}$ is a logic bug to fix, not a tolerance to widen.

Common Failure Modes and Diagnosis #

Dask Graph Explosion #

Symptom: The scheduler pegs one core for minutes before any work starts; the dashboard shows millions of tasks; the driver’s RAM climbs before computation.

Diagnosis: Chunks are too small, so a modest array becomes a graph of millions of nodes and the central scheduler chokes building and optimizing it. This is the classic anti-pattern of chunks=(1, 1) or per-row chunking. Confirm with len(arr.__dask_graph__()); a healthy DEM job has thousands of tasks, not millions. Fix by increasing chunk edge length so each task does meaningful work (target 100 MB–1 GB per chunk) and by avoiding element-wise Python in favour of vectorized array ops.

Ray Object-Store Spilling #

Symptom: Throughput collapses partway through a tile fleet; logs show Spilled N MiB to disk; local disk fills.

Diagnosis: Task results (or ray.put objects) accumulate in plasma faster than the driver consumes them, so Ray spills to disk and then thrashes. Almost always caused by submitting all futures at once and holding every result. Confirm from the Ray dashboard object-store panel. Fix with a bounded in-flight window and ray.wait (as in the fleet code above), and return small summaries from tasks rather than large arrays — write big outputs to object storage inside the task and return only a manifest.

Serialization of Large NumPy Arrays #

Symptom: Tasks spend most of their wall-clock time in transfer, not compute; CPU utilization is low while network is saturated.

Diagnosis: A large read-only array (SVP table, model weights, a reference grid) is captured by closure and re-serialized into every task. On Dask this is repeated cloudpickle over the wire; on Ray it defeats the zero-copy store. Confirm by checking task input sizes. Fix on Dask by client.scatter(big_array) once and passing the future; on Ray by ray.put once and passing the ObjectRef, never the raw array.

Version Skew Between Driver and Workers #

Symptom: Cryptic deserialization errors (ModuleNotFoundError inside a worker, pickle protocol mismatch, KeyError on an attribute), or silently wrong results after a rolling image update.

Diagnosis: Driver and workers run different library versions — a partially rolled-out Kubernetes deployment, or a local driver against a remote cluster on a different image. Both Dask and Ray assume homogeneous environments. Confirm by logging dask.__version__ / ray.__version__ and numpy.__version__ from a task and comparing to the driver. Fix by pinning one image for the whole cluster and forbidding driver-side ad-hoc installs; this is precisely what the containerization and image-pinning practices in the deployment domain enforce.

Pipeline Integration and Downstream Handoff #

Whichever engine you pick, it is a stage, not the whole system, and it hands off in two directions. Upstream, the point-cloud cleaning and artifact-and-noise removal stages feed cleaned tiles or scattered soundings into the scaled gridding step here. Downstream, the smoothed DEM flows to charting, volume, and hydrodynamic products.

Operationally, the engine runs inside two infrastructure layers. Job submission, ret/limits, and OOM handling belong to Kubernetes Job Scheduling for Bathymetric Processing: a Dask cluster or a Ray cluster is provisioned as pods with explicit memory requests and limits, and the same Handling OOM Kills in Bathymetric Kubernetes Jobs guidance applies to both — a spilled Ray object store and an under-chunked Dask worker are the two commonest causes of a 137 exit. Scheduling, retries, and dependency ordering belong to the orchestrator: Orchestrating Marine Pipelines with Dask & Prefect shows how a Prefect flow spins up an ephemeral Dask cluster for the gridding task and tears it down after, and the same pattern wraps a Ray job as a flow step. Reproducibility depends on the image: pin Dask, Ray, NumPy, and GDAL/PROJ exactly, as covered in the containerization cluster, so a re-run months later produces the identical DEM.

A practical hybrid many teams converge on: Ray for the per-tile cleaning and classification fan-out, Dask for the final out-of-core mosaic-and-smooth, joined by writing intermediate tiles to object storage and re-opening them as a chunked xarray cube. Each engine does the shape it is built for, and neither is forced to emulate the other.

{
  "engine": "dask",
  "dask_version": "2024.7.1",
  "numpy_version": "2.0.1",
  "chunk_edge_cells": 2048,
  "worker_memory_limit_gb": 8,
  "partition_invariance_max_delta": 3.1e-07,
  "output_zarr": "s3://survey/dem_smoothed.zarr",
  "crs": "EPSG:32618",
  "written_at": "2026-07-13T09:20:00Z"
}

The manifest is what lets a downstream consumer — or a future you — reproduce and trust the DEM: the engine and versions explain any numerical diff, and partition_invariance_max_delta is the proof the gate ran and passed.

Related #

• DEM Interpolation Techniques for Seafloor Mapping — the gridding algorithms this scales out-of-core
• Surface Smoothing Algorithms in Python — focal operations distributed here via Dask map_overlap
• Kubernetes Job Scheduling for Bathymetric Processing — provisioning Dask and Ray clusters with memory limits
• Orchestrating Marine Pipelines with Dask & Prefect — scheduling and tearing down ephemeral clusters around these stages

Up: Bathymetric Processing & Terrain Modeling