Troubleshooting DEM and Raster Export Failures

The final stage of a drone photogrammetry pipeline — turning a classified point cloud into a DEM/DSM and writing an orthomosaic to a distributable GeoTIFF or Cloud-Optimized GeoTIFF (COG) — is where the most expensive failures surface, because they surface late. Feature matching converged, bundle adjustment closed, the dense cloud looks correct in a viewer, and then writers.gdal emits an all-nodata raster, gdalwarp is killed by the OOM reaper, or the ortho opens as a solid black rectangle in QGIS. This page is a diagnostic aggregator for that stage. It does not re-teach how to build a DSM or write a COG — the DSM and DTM generation with PDAL and exporting Cloud-Optimized GeoTIFFs with rasterio guides do that. Instead it gives you a repeatable methodology: reproduce the failure deterministically, inspect the raster’s own metadata, isolate which pipeline stage introduced the fault, and apply the narrowest fix that resolves it.

The failure modes cluster into a small, enumerable set. PDAL rasterization produces empty or all-nodata output when the requested grid does not intersect the point cloud, or every point was filtered out upstream. Outputs carry the wrong CRS — or none at all — when a georeferencing string was dropped between the reconstruction engine and the writer. Vertical values are internally consistent but wrong by 20–40 metres when ellipsoidal and orthometric height are mixed, or feet are written where metres are declared. Memory blows up on large gdalwarp or writers.gdal calls because the whole raster is allocated at once. And a GeoTIFF that opens black is almost always a nodata, band-order, or overview problem rather than a corrupt file. Each of these has a deterministic Python-observable signature, and the sections below walk the methodology that turns a vague “the export is broken” into a specific, fixable stage fault.

Audience and prerequisites. This guide targets Python 3.10+ on a 64-bit OS with at least 8 GB RAM and command-line GDAL available on PATH. You should be comfortable reading gdalinfo output, know the difference between a projected and geographic CRS, and understand that a raster band’s nodata value, mask/alpha band, and colour interpretation are independent pieces of metadata that any one stage can corrupt. The diagnostic commands are non-destructive: they read metadata and statistics, never rewrite pixels, so you can run them freely against a suspect output before deciding on a fix.

# Prerequisites

Install the inspection and remediation stack in an isolated environment and pin the versions below. The single most common cause of “it works on my machine” in raster export is a GDAL version skew between the command-line tools and the rasterio/pdal wheels that link their own GDAL, because they can resolve nodata, COG structure, and CRS strings differently.

Confirm the GDAL versions agree before trusting any diagnosis: gdalinfo --version on the shell and rasterio.__gdal_version__ in Python should report the same major/minor. A mismatch there is itself a root cause, not a footnote. When the fault is specifically a COG that fails validation, cross-reference the dedicated fixing COG validation errors in GDAL workflow, which covers the validate_cloud_optimized_geotiff failure codes in depth.

# Conceptual architecture

Diagnosis is a decision tree, not a checklist. You start from a single observed symptom, inspect the raster’s declared metadata to narrow the candidate stages, isolate the offending stage by re-running it in isolation against a minimal input, and only then apply a fix — re-verifying against the same metadata you inspected. The tree below is the mental model the rest of this page implements. It deliberately separates symptom (what you see) from stage (where the fault was introduced) from fix (the narrowest change), because the same symptom — a black raster — can originate in rasterization, in the writer’s nodata handling, or in the overview build.

The methodology below is deterministic: each numbered step is a self-contained routine you run in order, and each narrows the search space before the next. Resist the temptation to skip to a fix — a nodata tweak applied to a raster whose real fault is a non-intersecting PDAL bounds will “succeed” while producing the same empty output.

# Step 1: Reproduce the failure deterministically

Before inspecting anything, capture the failure so it happens the same way every time. Intermittent OOM kills and “sometimes black” tiles are almost always deterministic once you pin the exact input, the GDAL environment, and the command. The routine below records the environment that GDAL actually sees and runs the failing export under a captured stderr, so the error text — ERROR 1: Cannot allocate memory, no data, Attempt to create ... failed — is preserved verbatim rather than paraphrased from memory.

import os
import subprocess
from pathlib import Path


def reproduce_export(cmd: list[str], workdir: Path) -> dict:
    """Run a failing export command with a pinned, recorded GDAL environment."""
    env = os.environ.copy()
    # Pin the cache so an OOM is reproducible rather than machine-dependent.
    env.setdefault("GDAL_CACHEMAX", "256")          # MB, not a fraction
    env.setdefault("CPL_DEBUG", "ON")               # verbose GDAL internals
    proc = subprocess.run(
        cmd, cwd=workdir, env=env,
        capture_output=True, text=True,
    )
    report = {
        "returncode": proc.returncode,              # -9 = SIGKILL (OOM reaper)
        "gdal_cachemax": env["GDAL_CACHEMAX"],
        "stderr_tail": proc.stderr.strip().splitlines()[-5:],
        "stdout_tail": proc.stdout.strip().splitlines()[-3:],
    }
    return report

A returncode of -9 means the kernel OOM-killer sent SIGKILL; the process never got to print an error, which is why the log looks truncated. That signature routes you straight to the memory branch and the decoding GDAL “Cannot allocate memory” on export deep-dive. Any other non-zero code with a message in stderr_tail means the writer ran far enough to complain, and the message text is your first classifier in Step 2.

# Step 2: Inspect the raster’s declared metadata

The output raster carries the evidence of its own defect. Before touching pixels, read the CRS, nodata value, band count, colour interpretation, mask flags, and block layout — these five fields distinguish a CRS fault from a nodata fault from an untiled-read fault. The function below returns a structured profile using rasterio, avoiding a subprocess round-trip and giving you typed values you can assert on.

import rasterio
from rasterio.enums import MaskFlags


def inspect_raster(path: str) -> dict:
    """Extract the metadata fields that classify an export failure."""
    with rasterio.open(path) as ds:
        band_stats = []
        for i in range(1, ds.count + 1):
            # Cheap sampled stats; approx=True avoids a full read on huge files.
            stats = ds.statistics(i, approx=True, clear_cache=True)
            band_stats.append({
                "band": i,
                "colorinterp": ds.colorinterp[i - 1].name,
                "min": stats.min, "max": stats.max, "mean": stats.mean,
                "mask": [f.name for f in ds.mask_flag_enums[i - 1]],
            })
        return {
            "crs": ds.crs.to_string() if ds.crs else None,   # None = no CRS!
            "is_projected": bool(ds.crs and ds.crs.is_projected),
            "nodata": ds.nodata,
            "dtype": ds.dtypes[0],
            "count": ds.count,
            "blockxsize": ds.block_shapes[0][1],             # 1 => untiled
            "overviews": ds.overviews(1),
            "all_valid": MaskFlags.all_valid in ds.mask_flag_enums[0],
            "bands": band_stats,
        }

Read the profile like a differential diagnosis. crs is None immediately explains a raster that “lands in the ocean” when overlaid — the CRS branch. A band whose min == max (often both 0 or both the nodata value) explains a black or flat raster without you opening a viewer. A blockxsize equal to the full raster width means the file is striped, not tiled, which forces GDAL to read entire scanlines and is the structural cause behind both the memory branch and slow COG serving. And all_valid == True on a raster that clearly has borders means no mask exists, so the fill value is being drawn as real data — the black-ortho branch.

# Step 3: Isolate the offending stage

With the symptom classified, re-run only the suspect stage against a minimal input so the fault reproduces without the rest of the pipeline confounding it. For a PDAL empty-raster case, the decisive check is whether the requested grid even intersects the point cloud: compare pdal info bounds against the pipeline’s output extent before blaming the writer. For a CRS case, isolate whether the georeferencing was present on the input to the writer or lost by it.

import json
import subprocess


def pdal_bounds(las_path: str) -> dict:
    """Return the point cloud's native bounds and CRS via `pdal info`."""
    out = subprocess.run(
        ["pdal", "info", "--summary", las_path],
        capture_output=True, text=True, check=True,
    )
    summary = json.loads(out.stdout)["summary"]
    bounds = summary["bounds"]
    return {
        "minx": bounds["minx"], "maxx": bounds["maxx"],
        "miny": bounds["miny"], "maxy": bounds["maxy"],
        "srs": summary.get("srs", {}).get("horizontal", ""),
        "num_points": summary["num_points"],
    }


def grid_intersects(cloud: dict, req_minx: float, req_miny: float,
                    req_maxx: float, req_maxy: float) -> bool:
    """True only if the requested raster window overlaps the cloud extent."""
    return not (req_maxx < cloud["minx"] or req_minx > cloud["maxx"]
                or req_maxy < cloud["miny"] or req_miny > cloud["maxy"])

If grid_intersects returns False, the raster is empty because the requested origin_x/origin_y or a hard-coded bounds string in the writers.gdal stage sits in a different coordinate space than the cloud — frequently a metres-vs-degrees confusion, or a UTM extent requested against a cloud still in a local frame. That is a rasterization-stage fault, and the fix lives in the DSM and DTM generation with PDAL guide, alongside the related fixing PDAL pipeline JSON schema errors guide when the pipeline definition itself is malformed. If num_points is near zero, an upstream filters.range or filters.smrf classification removed every point, and no writer setting can rasterize an empty cloud.

# Step 4: Apply the narrowest fix and re-verify

Only now do you change anything, and you change exactly one thing. Assigning a missing CRS is a metadata-only rewrite that must never touch pixels; setting a nodata value or building a mask changes how existing pixels are interpreted, not their values; capping memory changes only the writer’s environment. The routine below handles the two most common metadata-only fixes — stamping a known CRS onto a raster that lost it, and declaring the correct nodata — using rasterio’s update mode so the pixel data is preserved byte-for-byte.

import rasterio
from rasterio.crs import CRS


def repair_metadata(path: str, epsg: int | None = None,
                    nodata: float | None = None) -> None:
    """Stamp a known CRS and/or nodata onto an output without touching pixels."""
    with rasterio.open(path, "r+") as ds:      # r+ = in-place metadata update
        if epsg is not None and ds.crs is None:
            # Only assign when absent; never silently overwrite a real CRS.
            ds.crs = CRS.from_epsg(epsg)
        if nodata is not None:
            ds.nodata = nodata                 # declares fill, does not fill
    # Re-inspect to confirm the fix landed and nothing else changed.
    profile = inspect_raster(path)
    assert profile["crs"] is not None, "CRS still missing after repair"
    if nodata is not None:
        assert profile["nodata"] == nodata, "nodata not applied"

Assigning a CRS is correct only when the pixel grid was already in that projection but the label was dropped — it re-labels, it does not reproject. If the grid is genuinely in the wrong projection you must gdalwarp/reproject the pixels, which resamples them. For a vertical-datum fault (heights consistent but off by a fixed offset), the fix is a 3D or compound target CRS so the geoid separation is applied, exactly as covered in the coordinate transformation workflows in pyproj guide. For the memory and black-tile branches, the fixes are involved enough to warrant their own pages: decoding GDAL “Cannot allocate memory” on export and diagnosing empty or black orthomosaic tiles.

# Symptom-signature reference

This table is the fast path: match the error text or visible symptom in the first column, and it names the likely stage, the single command that confirms it, and the remediation. Use it to jump straight to the right numbered step rather than working the whole tree.

# Verification and output inspection

A fix is not done until the output passes the same metadata gate you used to diagnose it. The block below is a single acceptance function that asserts the four properties every deliverable raster must have: a real CRS, a declared nodata (or a mask), non-degenerate band statistics (proving the pixels are not all fill), and a tiled, over-viewed structure suitable for distribution. Run it in CI so a regression re-introducing any of these fails the build rather than a user’s viewer.

import rasterio
from rasterio.enums import MaskFlags


def assert_deliverable(path: str) -> None:
    """Gate a DEM/ortho export against the four must-have properties."""
    with rasterio.open(path) as ds:
        # 1. Georeferencing exists and is projected (metre-based distances).
        assert ds.crs is not None, "output has no CRS"
        assert ds.crs.is_projected, f"CRS is not projected: {ds.crs}"

        # 2. Fill is declared, either as nodata or an explicit mask/alpha.
        has_mask = any(MaskFlags.per_dataset in ds.mask_flag_enums[b]
                       or MaskFlags.alpha in ds.mask_flag_enums[b]
                       for b in range(ds.count))
        assert ds.nodata is not None or has_mask, "no nodata and no mask"

        # 3. Band statistics are non-degenerate (not an all-black raster).
        s = ds.statistics(1, approx=True)
        assert s.min != s.max, f"band 1 is flat ({s.min}); likely all-fill"

        # 4. Tiled with overviews for memory-safe reads and fast serving.
        assert ds.block_shapes[0][1] != ds.width, "raster is untiled (striped)"
        assert ds.overviews(1), "no overviews built"
    print(f"deliverable OK: {path}")

The s.min != s.max assertion is the cheapest possible black-raster detector and belongs in every export pipeline — it catches an all-nodata DSM and an all-zero ortho without a human ever opening QGIS. Pair this gate with the COG-specific validate_cloud_optimized_geotiff check when the deliverable is a COG, and with a round-trip CRS assertion when a reprojection was part of the fix. Choosing whether GDAL or rasterio drives the final write is itself a decision with tradeoffs, laid out in GDAL vs rasterio for DEM and orthomosaic export.

# Troubleshooting

PDAL’s writers.gdal runs without error but the output is entirely nodata. Where do I start? Do not touch the writer yet. Run pdal info --summary on the input and compare its bounds against the extent your pipeline requests. In the vast majority of cases the requested grid does not intersect the cloud — a UTM origin requested against a cloud still in a local or geographic frame — so every cell is nodata by definition. If the bounds do intersect, check num_points: an over-aggressive filters.range or ground filter may have removed every point before rasterization.

gdalinfo shows no coordinate system but the pixels look correct. Can I just assign one? Only if you are certain the pixel grid is already in that projection and the label was merely dropped. Assigning a CRS in r+ mode re-labels the file without moving a pixel, which is correct for a lost tag. If the grid is actually in a different projection, assigning the wrong CRS silently misplaces the raster; you must gdalwarp -t_srs to resample the pixels into the target grid instead.

Elevation values are internally consistent but every height is about 20 to 40 metres off against my GCPs. You are mixing ellipsoidal and orthometric height. The reconstruction produced ellipsoidal heights and the writer stamped a 2D CRS, so no geoid separation was applied. Rewrite the vertical reference with a 3D or compound CRS that names the geoid model, and re-check a band mean against a surveyed control point. A constant offset is the signature of a datum problem, not a per-pixel error.

My DEM heights are exactly 3.28 times too large. That factor is the metres-to-feet ratio. The source declared US survey feet (or the sensor exported feet) while the writer wrote a metre unit, or vice versa. Confirm the declared UNIT in gdalinfo against the band’s magnitude, then either convert the values with an explicit scale or correct the declared unit — never both.

gdalwarp is killed with no error message and a non-zero exit. A killed process with exit signal -9 is the kernel OOM reaper, not a GDAL bug. GDAL tried to allocate a whole-raster buffer or an oversized cache. Cap GDAL_CACHEMAX to a fixed megabyte value, add -wm to bound warp memory, and write tiled output with BIGTIFF=YES. The full recipe is in the memory-specific deep-dive linked in the reference table above.

The COG I wrote opens fine locally but fails validate_cloud_optimized_geotiff. Opening locally only proves the pixels are readable; COG validity is about structure — tiling, internal overviews, and IFD ordering. A file written untiled, or with external .ovr overviews, is a valid GeoTIFF but an invalid COG. Rewrite it tiled with internal overviews using the COG driver, following the dedicated COG validation workflow rather than patching the existing file.

# Related

• Generating DSM and DTM from Point Clouds with PDAL
• Exporting Cloud-Optimized GeoTIFFs with Rasterio
• GDAL vs Rasterio for DEM and Orthomosaic Export
• Decoding GDAL “Cannot Allocate Memory” on Export
• Diagnosing Empty or Black Orthomosaic Tiles

← DEM/DSM Generation & Raster Export Automation