What makes a GIS pipeline reproducible rather than just automated?

Reproducibility means anyone can rerun it from the pinned source data and the same software versions and get byte-comparable outputs. That requires deterministic commands, recorded GDAL/PROJ versions, no manual steps, and explicit input/output checksums — automation alone does not guarantee any of these.

When should I move from QGIS clicks to GDAL and SQL?

Move to the command line as soon as the same operation must run on many files, repeat on data updates, or be auditable. One-off exploratory work stays in QGIS; anything that ships or recurs belongs in a scripted GDAL/ogr2ogr or PostGIS pipeline.

How do I pin GDAL and PROJ versions so results don't drift?

Run the pipeline in a pinned container (e.g. a tagged osgeo/gdal image) or a locked conda environment, and record gdalinfo --version and projinfo output in the run log. PROJ datum-grid changes can shift coordinates by metres between versions.

How do I validate vector and raster outputs in a pipeline?

Add automated gates: check feature counts, ST_IsValid geometry, extent and CRS for vectors; check pixel dimensions, NoData, band count and statistics for rasters. Fail the run on mismatch rather than trusting that an output file means success.

Reproducible GIS Pipelines with GDAL and PostGIS

Short answer

A reproducible GIS pipeline is one that another engineer can clone, point at the same pinned inputs, run on a clean machine, and obtain the same outputs — without opening a GUI or guessing parameters. That requires four things the typical desktop workflow lacks: deterministic, scripted operations (GDAL/OGR, SQL); pinned software versions (GDAL, PROJ, PostGIS); explicit validation gates; and a recorded provenance trail. Automation gets you speed; reproducibility gets you trust, and they are not the same thing.

Automation is not reproducibility

A bash script that loops over files is automation. It becomes reproducible only when the result is stable across people, machines and time. Two things break that stability most often:

PROJ version drift. Datum transformations depend on grid shift files bundled with PROJ. Upgrading PROJ can change a transformation pipeline, moving coordinates by a metre or more. The same ogr2ogr -t_srs EPSG:2154 command can therefore produce different coordinates on two machines.
GDAL default changes. Resampling defaults, creation options and even output driver behaviour evolve between releases. A command that omits -r, -tr or -co inherits whatever the installed GDAL decides.

The fix is to remove ambiguity. Pin the environment, set every parameter explicitly, and record what ran.

Pin the environment

Use a tagged container so the toolchain is frozen:

docker run --rm -v "$PWD":/data ghcr.io/osgeo/gdal:ubuntu-small-3.9.2 \
  gdalinfo --version

or a locked conda environment (conda env export > environment.yml). At the start of every run, log the versions so the output is self-describing:

gdalinfo --version          # GDAL 3.9.2, released 2024/09/04
projinfo EPSG:2154          # the exact transformation pipeline in use

Commit environment.yml (or the image tag) alongside the code. A result without a recorded toolchain is not reproducible, only repeatable on your laptop.

Make every command deterministic

Compare a vague conversion with an explicit one. The reprojection below leaves nothing to defaults — target CRS, resampling, resolution, NoData, tiling and compression are all stated:

gdalwarp \
  -s_srs EPSG:4326 -t_srs EPSG:2154 \
  -r bilinear -tr 25 25 -tap \
  -dstnodata -9999 \
  -co TILED=YES -co COMPRESS=DEFLATE -co BIGTIFF=IF_SAFER \
  -overwrite \
  input.tif output_l93.tif

Use -tap (target-aligned pixels) so re-running on tiles produces a consistent grid. For vectors, name the output schema explicitly rather than relying on shapefile field truncation:

ogr2ogr -f GPKG output.gpkg input.shp \
  -t_srs EPSG:2154 \
  -nln boreholes -nlt PROMOTE_TO_MULTI \
  -lco GEOMETRY_NAME=geom -lco FID=fid

Prefer GeoPackage over shapefile: no 10-character field limit, no 2 GB ceiling, real types, and a single file.

Make the pipeline rerunnable with a build tool

The cleanest way to enforce "rerun from a clean source" is a dependency graph. A Makefile rebuilds only what changed and documents the DAG in one place:

output_l93.tif: input.tif
	gdalwarp -s_srs EPSG:4326 -t_srs EPSG:2154 -r bilinear -tr 25 25 -tap \
	  -overwrite $< $@

slope.tif: output_l93.tif
	gdaldem slope $< $@ -compute_edges

.PHONY: clean
clean:
	rm -f output_l93.tif slope.tif

make clean && make is the reproducibility test: if it does not rebuild the whole chain from source, the pipeline has a hidden manual step. For heavier workflows, the same principle applies with Snakemake, Luigi or a dvc stage file.

Database-side processing with PostGIS

When several products query the same spatial truth, push the logic into PostGIS so the database — not a scattering of scripts — is the source of record. Load with ogr2ogr straight into the DB:

ogr2ogr -f PostgreSQL "PG:dbname=geo host=localhost" parcels.gpkg \
  -nln parcels -lco GEOMETRY_NAME=geom -lco FID=id \
  -t_srs EPSG:2154

Then make transformations explicit and indexed:

-- Set the declared SRID on raw data that arrived without one
UPDATE raw_points SET geom = ST_SetSRID(geom, 4326)
WHERE ST_SRID(geom) = 0;

-- Reproject into the analysis CRS (this changes coordinates)
ALTER TABLE raw_points ADD COLUMN geom_l93 geometry(Point, 2154);
UPDATE raw_points SET geom_l93 = ST_Transform(geom, 2154);

CREATE INDEX idx_raw_points_geom_l93 ON raw_points USING GIST (geom_l93);

The single most common PostGIS mistake is using ST_SetSRID (which only relabels) where ST_Transform (which recomputes coordinates) is needed. ST_SetSRID is for data that is already in the target CRS but lost its SRID tag; ST_Transform is for actually moving coordinates between systems.

Validation gates

Treat outputs as guilty until proven valid. Bake checks into the pipeline so a bad run fails loudly:

Vectors: feature count within expected bounds; SELECT count(*) FROM t WHERE NOT ST_IsValid(geom) returns zero (repair with ST_MakeValid); CRS equals the target; extent inside the area of use.
Rasters: gdalinfo -stats confirms band count, pixel dimensions, NoData and a sane min/max; -checksum gives a stable hash for regression tests.
End to end: store an MD5 of each canonical output and diff it on the next run. A changed checksum that you did not intend is a regression.

A successful exit code and a file on disk prove only that the program did something, not that it did the right thing.

Worked example: nightly DEM refresh

A recurring elevation product, driven by make:

Pull new source tiles into src/ (pinned by a manifest with URLs and checksums).
gdalbuildvrt mosaic.vrt src/*.tif to virtually mosaic without copying.
gdalwarp the VRT to the project CRS and grid with -tap, fixed resolution and DEFLATE.
gdaldem hillshade and gdaldem slope from the warped DEM.
Run the validation script (extent, NoData, checksum diff). Abort and alert on any mismatch.
Load tiles into PostGIS as out-of-db rasters and refresh the materialised views the web product reads from.

Because every step is in the Makefile with pinned inputs and a pinned GDAL image, the same DEM comes out tomorrow as today — unless the source data genuinely changed, which is the only difference you want.

Bathyl perspective

We reach for data engineering when GIS stops being a deliverable and becomes infrastructure. The aim is not sophistication for its own sake; it is a workflow that survives new data, new analysts and new products, where any result can be traced back to its source command and software version.

Reproducible GIS Pipelines