From Drone Capture to Publishable Map: A Python + GDAL Pipeline

Drone survey output is not a map. It’s a folder of overlapping JPEG images with GPS metadata embedded in the EXIF. Getting from that folder to a georeferenced raster you can publish involves several steps that each look simple and collectively take longer than expected the first time.

This is the pipeline I use. It assumes you’ve already run photogrammetry (I use OpenDroneMap) and have a GeoTIFF output. The processing from that point forward is where most of the actual work lives.

What You Start With

After OpenDroneMap processes a typical flight, you get something like:

odm_orthophoto/
  odm_orthophoto.tif      # full-resolution GeoTIFF, often 2–8 GB
  odm_orthophoto.png      # web-resolution preview
odm_dem/
  dtm.tif                 # Digital Terrain Model
  dsm.tif                 # Digital Surface Model
odm_georeferencing/
  odm_georeferencing_model_geo.txt

The odm_orthophoto.tif is your main product. It’s a 4-band (RGBA) GeoTIFF in the CRS you specified during processing — usually WGS84 (EPSG:4326) or UTM for survey work.

Step 1: Check the Input

from osgeo import gdal
import json

gdal.UseExceptions()

def inspect_raster(path: str) -> dict:
    ds = gdal.Open(path)
    if ds is None:
        raise ValueError(f"Could not open {path}")

    gt = ds.GetGeoTransform()
    proj = ds.GetProjection()

    return {
        "driver": ds.GetDriver().ShortName,
        "size": (ds.RasterXSize, ds.RasterYSize),
        "bands": ds.RasterCount,
        "crs": proj[:80] + "..." if len(proj) > 80 else proj,
        "origin": (gt[0], gt[3]),
        "pixel_size": (gt[1], abs(gt[5])),
        "nodata": [ds.GetRasterBand(i+1).GetNoDataValue() for i in range(ds.RasterCount)],
    }

info = inspect_raster("odm_orthophoto.tif")
print(json.dumps(info, indent=2))

The GetGeoTransform() tuple deserves attention: gt[0] and gt[3] are the top-left corner coordinates. gt[1] is the pixel width in map units. gt[5] is the pixel height — it’s negative because rasters count rows from top to bottom but coordinate systems count up from the bottom. New people get confused by this every time.

Step 2: Reproject to Target CRS

If you’re delivering to a client expecting Texas State Plane or if you need to overlay with existing spatial data in a different CRS:

from osgeo import osr
import subprocess

def reproject_raster(src_path: str, dst_path: str, target_epsg: int) -> None:
    """Reproject a raster to a target CRS using gdalwarp."""
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(target_epsg)

    warp_options = gdal.WarpOptions(
        dstSRS=srs.ExportToWkt(),
        resampleAlg=gdal.GRA_Bilinear,
        multithread=True,
        warpMemoryLimit=4096,  # MB
        creationOptions=["COMPRESS=LZW", "TILED=YES", "BIGTIFF=YES"],
    )

    result = gdal.Warp(dst_path, src_path, options=warp_options)
    if result is None:
        raise RuntimeError(f"Warp failed for {src_path}")
    result = None  # flush and close

    print(f"Reprojected to EPSG:{target_epsg}: {dst_path}")

# Texas State Plane Central (NAD83) = EPSG:32139
reproject_raster(
    "odm_orthophoto.tif",
    "orthophoto_stateplane.tif",
    32139
)

BIGTIFF=YES is important once your file exceeds 4GB. Without it, GDAL will fail silently or produce a corrupt file. TILED=YES makes the file accessible randomly without reading the whole thing into memory — critical for files larger than your RAM.

Step 3: Generate a Cloud-Optimized GeoTIFF (COG)

A regular GeoTIFF requires the whole file to be downloaded before it can be displayed. A COG stores internal tiles and overviews in a specific order that allows HTTP range requests — a viewer can fetch just the tiles it needs for the current zoom level.

def convert_to_cog(src_path: str, cog_path: str) -> None:
    translate_options = gdal.TranslateOptions(
        format="GTiff",
        creationOptions=[
            "COMPRESS=LZW",
            "TILED=YES",
            "BLOCKXSIZE=512",
            "BLOCKYSIZE=512",
            "COPY_SRC_OVERVIEWS=YES",
            "BIGTIFF=YES",
        ],
    )

    # Build overviews on the source first
    src = gdal.Open(src_path, gdal.GA_Update)
    overview_levels = [2, 4, 8, 16, 32, 64]
    src.BuildOverviews("AVERAGE", overview_levels)
    src = None

    # Translate to COG
    gdal.Translate(cog_path, src_path, options=translate_options)
    print(f"COG written: {cog_path}")

convert_to_cog("orthophoto_stateplane.tif", "orthophoto_stateplane_cog.tif")

BuildOverviews is the slow step on large files — expect 10–20 minutes on a 4GB raster. Run it with GDAL_CACHEMAX set high:

export GDAL_CACHEMAX=4096
python process.py

Step 4: Extract the DEM to Contours

For survey deliverables, contour lines at 1-foot or 2-foot intervals are often required alongside the orthophoto:

def raster_to_contours(dem_path: str, shp_path: str, interval: float) -> None:
    dem_ds = gdal.Open(dem_path)
    band = dem_ds.GetRasterBand(1)

    drv = ogr.GetDriverByName("ESRI Shapefile")
    contour_ds = drv.CreateDataSource(shp_path)
    srs = osr.SpatialReference()
    srs.ImportFromWkt(dem_ds.GetProjection())
    contour_layer = contour_ds.CreateLayer("contours", srs=srs)

    field_defn = ogr.FieldDefn("Elevation", ogr.OFTReal)
    contour_layer.CreateField(field_defn)

    gdal.ContourGenerate(band, interval, 0, [], False, 0, contour_layer, -1, 0)

    contour_ds = None
    dem_ds = None
    print(f"Contours written to {shp_path}")

from osgeo import ogr
raster_to_contours("odm_dem/dtm.tif", "contours_1ft.shp", interval=0.3048)  # 1 foot in meters

On Error Handling in Production

The GDAL Python bindings are thin wrappers and their error behavior is inconsistent. gdal.UseExceptions() at the top of your script converts GDAL errors to Python exceptions rather than returning None silently. Without it, gdal.Open("nonexistent.tif") returns None and you get a NoneType has no attribute crash several lines later with no useful context.

Log the input file size, band count, and CRS at the start of every run. When something goes wrong at step 3 of a 40-minute pipeline, you want to know whether the input was valid.