GDAL vs Rasterio for DEM and Orthomosaic Export

Both tools sit on top of the same C++ library, so this is not a question of which one is “more powerful” — a raster GDAL can write, rasterio can write identically, because rasterio is GDAL with a Python skin. The decision is about the interface you build your export stage around: the gdal_translate/gdaladdo command-line tools and the osgeo.gdal Python bindings on one side, versus rasterio’s pythonic dataset objects and windowed array I/O on the other. For a drone survey that turns point clouds into COGs, that interface choice determines how you script batch jobs, how you keep a 40-gigapixel mosaic inside a memory budget, how you unit-test the export, and how cleanly the export composes with the rest of a NumPy-based pipeline. This page runs the same DEM-and-orthomosaic COG export through both paths, side by side, and then lays out where each genuinely wins.

This is the decision layer of the DEM/DSM generation and raster export pipeline. If you have already picked rasterio and want the full authoring recipe, it lives in exporting Cloud-Optimized GeoTIFFs with Rasterio; the choice here is analogous to the API-versus-CLI trade-off examined for reprojection in PyProj vs GDAL for coordinate transforms.

Audience and prerequisites. This targets Python 3.10+ engineers who already know what a COG is (internal tiling, overviews, IFD ordering) and are deciding how to build the export stage of an automated pipeline. You should be comfortable reading both a shell invocation and a NumPy array loop, because the honest comparison requires seeing both. Both paths require a GDAL ≥ 3.4 with the COG driver; rasterio bundles its own in the wheel, while the CLI path uses a system or conda GDAL.

# Prerequisites

The two paths differ in what you install and how you invoke it. The CLI path needs the GDAL binaries on PATH (or the osgeo.gdal bindings for in-process calls); the rasterio path needs only the pip-installable wheel, which carries its own GDAL and never touches a system install.

Whichever path you pick, validate its output with the same cog_validate so the two are held to one standard. A subtle deployment note: mixing a system GDAL on PATH with rasterio’s bundled GDAL in one environment can produce two different GDAL versions in the same pipeline, which is exactly the kind of drift that makes a COG validate on one machine and fail on another — pin one source of GDAL per environment.

# Conceptual architecture

The key insight is that both paths funnel into the identical GDAL COG driver, so the output can be byte-identical; what differs is everything upstream of the driver. The CLI path shells out to gdal_translate, passing creation options as -co strings, and GDAL handles the read, resample, compress, and write internally as a black box. The rasterio path opens the source as a Python object, optionally reads and processes pixels window-by-window as NumPy arrays, and writes through the same driver — giving you a hook to inspect or transform pixels mid-flight that the CLI does not expose.

The sections below do the same export both ways — first the CLI, then rasterio — so you compare like for like, then a decision table and troubleshooting distill when each is the right call.

# Step 1: Export a DEM COG the GDAL CLI way

The CLI path is the shortest route from a finished raster to a COG, and it is the natural choice inside shell-driven batch pipelines, Makefiles, or containers where adding a Python dependency is friction. gdal_translate -of COG does the whole job — tiling, overviews, compression, IFD ordering — in one process, driven entirely by -co (creation option) flags.

# Float32 DEM -> COG: ZSTD + float predictor, average overviews.
gdal_translate -of COG \
  -co COMPRESS=ZSTD \
  -co PREDICTOR=3 \
  -co BLOCKSIZE=512 \
  -co OVERVIEW_RESAMPLING=AVERAGE \
  -co BIGTIFF=IF_SAFER \
  dem_source.tif dem_cog.tif

# Confirm the structure with the bundled validator.
python -m osgeo_utils.samples.validate_cloud_optimized_geotiff dem_cog.tif

Every knob is a string flag, which is the CLI’s strength and its weakness: it is trivial to template into a job runner, but there is no type checking, no autocompletion, and a typo like PREDICTER=3 is silently ignored rather than raising. For fully automated pipelines the same options are available in-process through the gdal.Translate binding, which avoids spawning a subprocess:

from osgeo import gdal

gdal.UseExceptions()   # make GDAL raise instead of returning None
gdal.Translate(
    "dem_cog.tif", "dem_source.tif",
    format="COG",
    creationOptions=["COMPRESS=ZSTD", "PREDICTOR=3", "BLOCKSIZE=512",
                     "OVERVIEW_RESAMPLING=AVERAGE", "BIGTIFF=IF_SAFER"],
)

gdal.UseExceptions() is essential here — without it the bindings return None on failure and a broken export slips through silently. This in-process form is the honest CLI-equivalent to compare against rasterio, because both run inside Python.

# Step 2: Export the same DEM COG the rasterio way

rasterio expresses the identical export as a copy through the COG driver, but as Python objects with real attributes rather than opaque strings. The payoff is not this trivial case — it is that the source is a first-class dataset you can inspect (src.dtypes, src.nodata, src.crs) and branch on before writing.

import rasterio
from rasterio.shutil import copy as rio_copy

def export_dem_cog(src_path: str, dst_path: str) -> None:
    with rasterio.open(src_path) as src:
        # Branch on real metadata: float DEM -> predictor 3, else 2.
        predictor = 3 if src.dtypes[0].startswith("float") else 2
    rio_copy(
        src_path, dst_path,
        driver="COG",
        compress="ZSTD",
        predictor=predictor,
        blocksize=512,
        overview_resampling="average",
        BIGTIFF="IF_SAFER",
    )

The decision that picks predictor from the dtype is the whole point: in the CLI path you would hardcode PREDICTOR=3 in the command and hope every input really is float, or write shell around gdalinfo output to branch. In rasterio the metadata is right there as typed Python. This is the same data-driven option selection the COG export guide builds on across compression and overview choices.

# Step 3: Where rasterio pulls ahead — windowed processing during export

The comparison stops being cosmetic when the export must transform pixels, not just re-encode them: clamping a DEM to a valid elevation range, applying a hillshade-friendly NoData, or masking an orthomosaic — on a raster too large to hold in RAM. The CLI treats the raster as a black box, so any pixel edit means a separate pre-processing pass and an intermediate file. rasterio reads and writes in windows, so you stream tile-sized NumPy arrays through your own function and never materialize the whole raster.

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

def export_clamped_dem(src_path: str, dst_path: str,
                       lo: float = -100.0, hi: float = 9000.0) -> None:
    """Clamp DEM elevations tile-by-tile while streaming to a tiled GeoTIFF."""
    with rasterio.open(src_path) as src:
        profile = src.profile.copy()
        profile.update(driver="GTiff", tiled=True, blockxsize=512,
                       blockysize=512, compress="ZSTD", predictor=3,
                       BIGTIFF="IF_SAFER")
        with rasterio.open(dst_path, "w", **profile) as dst:
            # Iterate the source's native tile grid; peak RAM = one block.
            for _, window in src.block_windows(1):
                data = src.read(1, window=window)
                nodata = src.nodata
                mask = data != nodata if nodata is not None else np.ones_like(data, bool)
                data[mask] = np.clip(data[mask], lo, hi)   # per-tile edit
                dst.write(data, 1, window=window)
    # Note: this writes a tiled GTiff; re-emit through the COG driver for
    # correct IFD ordering, e.g. rio_copy(dst_path, cog_path, driver="COG", ...).

This is the capability the CLI cannot match without an extra tool and an extra file on disk: arbitrary per-tile computation woven into the export, with peak memory pinned to one block regardless of the mosaic’s size. It is the same windowed discipline that keeps large exports inside a memory budget, exactly the failure mode diagnosed in decoding GDAL “cannot allocate memory” on export. Note the trade-off in the closing comment: a manual windowed write produces a tiled GeoTIFF, not a COG, so you still re-emit through the COG driver to fix IFD ordering — the CLI would have given you the COG in one step if you did not need the per-pixel edit.

# Comparison and parameter table

Head to head on the dimensions that actually decide an export stage. Neither column is universally better; the right pick depends on which row matters most to your pipeline.

# Verification and output inspection

Whichever path produced the file, validate it against the same standard so the comparison is honest. Assert the structural COG guarantees and that the codec/predictor you asked for actually landed — GDAL silently ignores an unsupported option on either path.

import rasterio
from rio_cogeo.cogeo import cog_validate

def verify_export(path: str, want_predictor: str = "3") -> None:
    with rasterio.open(path) as ds:
        assert ds.profile["tiled"] is True, "not tiled"
        assert len(ds.overviews(1)) >= 3, f"too few overviews: {ds.overviews(1)}"
        struct = ds.tags(ns="IMAGE_STRUCTURE")
        assert struct.get("PREDICTOR", want_predictor) == want_predictor, \
            f"predictor not applied: {struct.get('PREDICTOR')}"
    ok, errors, _ = cog_validate(path)
    assert ok, f"not a valid COG: {errors}"
    print("verified:", path)

Running this identical check on the CLI output and the rasterio output is the fastest way to prove they converged on the same file — if both pass with the same IMAGE_STRUCTURE tags and overview factors, the choice really was about the interface, not the result.

# Troubleshooting

The CLI and rasterio produce different-sized COGs from the same DEM. One path applied the float predictor and the other did not. Check the IMAGE_STRUCTURE tags on both — a common cause is PREDICTOR=3 mistyped or omitted in a -co flag, which GDAL ignores silently, while the rasterio kwarg would have raised. Align the options and the sizes converge.

gdal.Translate returns without error but the output is missing or wrong. You did not call gdal.UseExceptions(), so the binding signalled failure by returning None instead of raising. Enable exceptions at import time; this is the single most common reason a gdal.Translate pipeline “succeeds” yet produces no valid file.

A COG validates on the machine that used rasterio but fails where the CLI ran (or vice versa). The two environments ship different GDAL versions — rasterio’s bundled wheel versus a system/conda GDAL. Overview offset handling changed across GDAL releases, so pin one GDAL source per environment and confirm with rasterio.gdal_version() and gdalinfo --version.

The rasterio windowed write produces a tiled file that fails COG validation. A manual block_windows write yields a tiled GeoTIFF, not a COG — it has no correctly ordered overview IFDs. Re-emit the tiled result through driver="COG" (or gdal_translate -of COG) to add overviews and fix ordering, as noted in Step 3.

Batch throughput is fine on the CLI but the rasterio job is much slower. You are likely reading windows that do not align to the source’s native block grid, forcing repeated partial-block decompression. Iterate src.block_windows(1) (the native tiling) rather than an arbitrary fixed window, so each read maps to whole compressed blocks.

A multispectral uint16 ortho exports fine via DEFLATE but errors under COMPRESS=JPEG. JPEG in a COG is 8-bit only, on both paths. This is not a tool difference — neither the CLI nor rasterio can JPEG-compress 16-bit data. Use DEFLATE or ZSTD for anything beyond 8-bit, as the compression tuning guide details.

# Related

• Exporting Cloud-Optimized GeoTIFFs with Rasterio
• Generating DSM and DTM from Point Clouds with PDAL
• Troubleshooting DEM and Raster Export Failures
• PyProj vs GDAL for Coordinate Transforms

← DEM/DSM Generation & Raster Export Automation