GDAL & OGR Command-Line Automation Cookbook
Production-grade recipes for high-performance geospatial data engineering: multi-threaded gdalwarp reprojection, Cloud-Optimized GeoTIFF (COG) generation, VRT catalogs, and ogr2ogr spatial SQL.
The universal geospatial data abstraction library powering QGIS, PostGIS, Google Earth, and cloud spatial engines.
Streaming byte-range HTTP reads direct from AWS S3, Google Cloud Storage, and web servers without downloading.
Multi-threaded I/O chunking and multi-core warp filtering kernels maximizing server hardware utilization.
SpatiaLite embedded SQL engine executing complex spatial joins and spatial indexing directly in terminal pipelines.
Recipe 1: High-Performance Multi-Threaded Reprojection (gdalwarp)
Reprojecting large multi-gigabyte raster datasets from geographic coordinates (EPSG:4326) to projected metric frames (such as UTM Zone 33N, EPSG:32633) is often a computational bottleneck in automated pipelines. Out-of-the-box, gdalwarp runs single-threaded with a default 64 MB memory cache.
Resampling Algorithm Selection Guide
Selecting the proper resampling kernel is critical to preserve numerical integrity:
- Near (Nearest Neighbor): Mandatory for discrete categorical rasters (land cover classifications, soil classes, binary masks). Never creates new interpolated values.
- Bilinear: Fast linear interpolation. Suitable for smooth continuous imagery where slight blurring is acceptable.
- Cubic (Cubic Convolution): Evaluates a $4 \times 4$ pixel neighborhood. Industry standard for satellite surface reflectance and optical orthomosaics. Minimizes edge artifacts.
- Lanczos: High-order sinc windowing kernel evaluating a $6 \times 6$ pixel neighborhood. Optimal for sharp elevation models (DEMs) and sub-meter aerial photography.
- Average: Averages all contributing pixel values. Mandatory when downsampling high-resolution rasters to coarser cell sizes to avoid aliasing.
Production CLI Command
# High-throughput multi-threaded warp recipe
gdalwarp \
-s_srs EPSG:4326 \
-t_srs EPSG:32633 \
-r cubic \
-tr 10.0 10.0 \
-dstnodata -9999 \
-multi \
-wo NUM_THREADS=ALL_CPUS \
-wm 2048 \
--config GDAL_CACHEMAX 2048 \
-co COMPRESS=DEFLATE \
-co PREDICTOR=2 \
-co TILED=YES \
-co BLOCKXSIZE=512 \
-co BLOCKYSIZE=512 \
/data/raw_scene_4326.tif \
/data/warped_utm33n.tif
Parameter Breakdown:
-multi: Enables multi-threaded processing across independent chunks.-wo NUM_THREADS=ALL_CPUS: Spawns worker threads across all available CPU cores for the resampling convolution filter.-wm 2048: Allocates 2 gigabytes of RAM specifically for the internal warp buffer.--config GDAL_CACHEMAX 2048: Sets GDAL's global block cache to 2,048 MB, avoiding disk thrashing on large datasets.
Recipe 2: Generating Cloud-Optimized GeoTIFFs (COG)
Cloud-Optimized GeoTIFFs (COGs) are the foundational storage format of modern cloud-native geospatial stacks. A COG is an ordinary TIFF file organized with an internal tiled layout and hierarchical reduced-resolution overviews (pyramids). This allows client applications (such as QGIS, web map viewers, or Python rasterio) to fetch only the byte ranges corresponding to the visible map bounding box using standard HTTP GET range requests.
Single-Step Modern COG Generation
In modern GDAL (version 3.1 and later), the dedicated COG driver creates internal tiles, computes overviews, and applies lossless compression in a single command:
# Convert standard raster to Cloud-Optimized GeoTIFF with Zstandard compression
gdal_translate \
-of COG \
-co COMPRESS=ZSTD \
-co PREDICTOR=2 \
-co OVERVIEWS=AUTO \
-co RESAMPLING=CUBIC \
-co BLOCKSIZE=512 \
-co NUM_THREADS=ALL_CPUS \
/data/input_orthomosaic.tif \
/data/production_cog.tif
# Validate COG internal structure and byte offsets
python3 -m pip install -q validate_cloud_optimized_geotiff
python3 -m validate_cloud_optimized_geotiff /data/production_cog.tif
Using COMPRESS=ZSTD (Zstandard) yields superior decompression speeds and equal or higher compression ratios compared to traditional DEFLATE, significantly reducing egress bandwidth costs in cloud object stores (such as AWS S3 or Google Cloud Storage).
Recipe 3: Virtual Raster Datasets (VRT) for Seamless Mosaicking
When working with aerial survey projects comprising hundreds of individual flight orthophoto tiles, generating a physical merged GeoTIFF can consume hundreds of gigabytes of disk space and require hours of redundant processing. The Virtual Raster format (VRT) solves this by building an XML description file that references source files on disk.
Building a Spatial Mosaic Catalog
# Build a seamless spatial mosaic of all GeoTIFFs in a folder in seconds
gdalbuildvrt \
-resolution highest \
-r cubic \
-srcnodata -9999 \
-vrtnodata -9999 \
/data/orthos/regional_mosaic.vrt \
/data/orthos/*.tif
# Build a multi-spectral stacked raster from individual single-band files
gdalbuildvrt \
-separate \
/data/sentinel2_stack.vrt \
/data/B02_blue.tif \
/data/B03_green.tif \
/data/B04_red.tif \
/data/B08_nir.tif
The resulting .vrt file is lightweight (typically a few kilobytes) and opens instantly in QGIS or GDAL CLI tools as if it were a single monolithic multi-gigabyte raster image.
Recipe 4: High-Speed Vector Spatial SQL Filtering (ogr2ogr)
The ogr2ogr utility provides format conversion, reprojection, and spatial filtering across vector geodatabases (Shapefile, GeoJSON, GeoPackage, FlatGeobuf, PostGIS). Leveraging the embedded SQLite/SpatiaLite SQL engine allows filtering and aggregating attributes at the C++ level before writing to disk.
Spatial Filtering and Attribute Dissolve Recipe
# Convert Shapefile to high-performance GeoPackage with Spatialite SQL query
ogr2ogr \
-f "GPKG" \
-s_srs EPSG:4326 \
-t_srs EPSG:32633 \
-dialect SQLite \
-sql "SELECT
parcel_id,
owner_type,
ST_Area(GEOMETRY) as area_sqm,
GEOMETRY
FROM cadastral_parcels
WHERE zoning = 'AGRICULTURAL'
AND ST_Area(GEOMETRY) >= 10000" \
-nln "priority_solar_parcels" \
/data/output_parcels.gpkg \
/data/cadastral_parcels.shp
This single command reads an unprojected shapefile, filters parcels larger than 1 hectare classified as agricultural, computes their projected metric surface area, reprojects boundary coordinates into UTM Zone 33N, and stores the output in a single-file GeoPackage container with spatial indexes automatically built.
Recipe 5: Command-Line Band Math with gdal_calc.py
To execute complex spectral index equations without launching a desktop GUI, gdal_calc.py evaluates mathematical expressions on raster arrays:
# Calculate the Sealed Urban Index (SUI) directly in terminal:
# SUI = ((Green - SWIR1) * (SWIR2 - Red)) / ((Green + SWIR1) * (SWIR2 + Red)) + 0.5
gdal_calc.py \
-A /data/sentinel2_B03_green.tif \
-B /data/sentinel2_B04_red.tif \
-C /data/sentinel2_B11_swir1.tif \
-D /data/sentinel2_B12_swir2.tif \
--outfile=/data/sentinel2_SUI.tif \
--calc="(((A - C) * (D - B)) / (numpy.maximum(A + C, 1e-6) * numpy.maximum(D + B, 1e-6))) + 0.5" \
--type=Float32 \
--NoDataValue=-9999 \
--co="COMPRESS=DEFLATE" \
--co="TILED=YES"
Common Pitfalls and Field Troubleshooting in GDAL/OGR
Terminal geospatial pipelines present subtle configuration hazards:
1. The In-Place Overwrite Trap
Never set the output file path identical to the input file path in gdal_translate, gdalwarp, or ogr2ogr:
# CATASTROPHIC ERROR: Truncates and destroys the source file immediately!
gdalwarp -t_srs EPSG:32633 input.tif input.tif
GDAL opens the output file in write mode before reading input blocks. This truncates the source file to 0 bytes, permanently corrupting your dataset. Always output to a temporary scratch path or dedicated directory.
2. Degree vs. Meter Resolution Disparity
When reprojecting from geographic coordinates (EPSG:4326) to UTM, forgetting to specify the target resolution (-tr) causes GDAL to preserve the original angular pixel size in projected metric space:
# If source has dx = 0.0001 degrees, target will have dx = 0.0001 METERS (0.1 mm)!
# This creates a gigantic multi-terabyte raster that immediately exhausts disk storage.
Whenever converting from EPSG:4326 to a projected coordinate system, explicitly specify the metric grid resolution: -tr 10.0 10.0 (for 10m Sentinel-2) or -tr 30.0 30.0 (for 30m Landsat).
3. Axis Order Flips in GDAL 3 / PROJ 6+
Starting in GDAL 3.0, the EPSG geodetic authority rules are enforced strictly. In EPSG:4326, latitude comes first (Northing, Easting). If you pipe data into scripts expecting traditional GIS order (Easting, Northing), coordinates will be inverted. Resolve this across any command by setting:
export OGR_CT_FORCE_TRADITIONAL_GIS_ORDER=YES
Frequently Asked Questions on GDAL & OGR Processing
Why does GDAL 3.x swap latitude and longitude coordinate axes unexpectedly compared to GDAL 2.x?
Starting in GDAL 3.0 and PROJ 6, GDAL strictly enforces the authoritative coordinate axis order published by the EPSG registry. Geographic coordinate reference systems like EPSG:4326 define latitude (Northing) as the first axis and longitude (Easting) as the second axis. GDAL 2.x historically forced longitude first (X, Y). To maintain consistent (Easting, Northing / Lon, Lat) behavior across scripts, developers can pass the flag: --config OGR_CT_FORCE_TRADITIONAL_GIS_ORDER YES.
What configuration parameters maximize multi-core CPU throughput in gdalwarp?
By default, gdalwarp runs single-threaded with a modest 64 MB memory cache. To maximize hardware utilization across high-core workstations or cloud nodes, supply: -multi -wo NUM_THREADS=ALL_CPUS -wm 2048 --config GDAL_CACHEMAX 2048. The -multi flag enables multi-threaded chunk processing, -wo NUM_THREADS=ALL_CPUS parallelizes warp filtering kernels, and -wm allocates 2048 MB of RAM specifically for the warping buffer.
How does predictor compression (PREDICTOR=2) reduce Cloud-Optimized GeoTIFF file size?
Predictor compression applies a reversible mathematical filter before DEFLATE or LZW compression. Horizontal differencing (PREDICTOR=2) replaces each pixel value with the numerical difference between itself and its horizontal predecessor. Because adjacent raster pixels in continuous imagery share similar digital numbers, their difference values hover close to zero. This converts complex high-entropy data into long sequences of zeros, enabling the compression algorithm to reduce file sizes by an additional 20% to 40%.
Why is a Virtual Raster (VRT) superior to generating a physical merged mosaic of multiple tiles?
Merging hundreds of gigabytes of GeoTIFF tiles into a single physical raster requires hours of disk I/O and duplicates terabytes of storage. A Virtual Raster (VRT) is a lightweight XML file generated in seconds that points to the original tile files on disk. GDAL and desktop GIS software treat the VRT file as a seamless, continuous raster, reading only the specific bounding box pixels required during display or processing.
How does streaming remote rasters with /vsicurl/ eliminate local download bottlenecks?
The GDAL virtual file system handler /vsicurl/ wraps HTTP/HTTPS range requests. When opening a Cloud-Optimized GeoTIFF hosted on an Amazon S3 bucket, Google Cloud Storage, or HTTP web server, GDAL reads the file header (typically 50 KB) to parse metadata and internal tile byte offsets. When querying a local point or bounding box, GDAL issues targeted byte-range GET requests to download only the specific tiles needed, avoiding downloading gigabyte-scale source files.