How to Make a GIS Workflow Reproducible in Python

On this page

A reviewer asks a simple question: "can you regenerate the figures from last year's report?" You still have the script. You still have the data — probably; there is a folder called final_v3 and one called final_v3_actual. You run it, the libraries have moved on, the supplier's file behind the same URL may or may not be the one you used, and the numbers come out different. Nothing was done wrong, exactly. The workflow was simply never built to be repeated. Reproducibility is the discipline of making sure that "run it again" produces the same answer — this year, next year, and on someone else's machine.

Problem statement

You cannot reliably reproduce your own results. The reasons compound:

The environment moved — a library upgrade, or a setting as small as PROJ's network access, changed a result and nothing recorded what produced the original.
The data moved — an input file was replaced at the same URL or path, so the version used back then no longer exists anywhere.
The parameters are gone — the threshold used for that run was a literal in a script that has since been edited.
The code is ambiguous — "the script" exists in three copies with slightly different logic and no version history.
Randomness — sampling, clustering, or a tie-break in a join returns a different answer each run.

A reproducible run records — and can restore — the code, environment, inputs and parameters behind each output, so the same command gives the same numbers a year later. Everything below runs on a real, small pipeline: HM Land Registry's 21,253 INSPIRE parcels for Rutland, filtered by area and sampled with a seed.

Quick answer

Pin five things and record all five with every run: the environment (locked dependency versions), the code (a commit hash), the inputs (immutable, checksummed files), the parameters (a config file), and any randomness (a fixed seed). Write them into a manifest stored beside the output.

Miss any one layer and "run it again" stops being a reliable question.

def file_hash(path: Path, chunk=1 << 20) -> str:
    h = hashlib.sha256()
    with Path(path).open("rb") as fh:
        for block in iter(lambda: fh.read(chunk), b""):
            h.update(block)
    return h.hexdigest()

def git_commit() -> str:
    try:
        out = subprocess.run(["git", "rev-parse", "--short", "HEAD"],
                             capture_output=True, text=True, check=True).stdout.strip()
        dirty = subprocess.run(["git", "status", "--porcelain", "--untracked-files=no"],
                               capture_output=True, text=True).stdout.strip()
        return out + ("-dirty" if dirty else "")
    except Exception:
        return "unknown"

def environment_report() -> dict:
    import pandas as pd, pyogrio, pyproj, shapely
    return {
        "python": sys.version.split()[0], "platform": platform.platform(),
        "geopandas": gpd.__version__, "pandas": pd.__version__, "shapely": shapely.__version__,
        "geos": shapely.geos_version_string, "pyproj": pyproj.__version__, "proj": pyproj.proj_version_str,
        "pyogrio": pyogrio.__version__, "gdal": pyogrio.__gdal_version_string__,
        "proj_network": pyproj.network.is_network_enabled(),
    }

def write_manifest(run_dir: Path, cfg: dict, inputs, outputs, stats: dict):
    manifest = {
        "run_at": datetime.now(timezone.utc).isoformat(timespec="seconds"),
        "code": {"commit": git_commit(), "entrypoint": "pipeline.py"},
        "environment": environment_report(),
        "config": json.loads(json.dumps(cfg, default=str)),
        "inputs": [{"path": str(p), "sha256": file_hash(p), "bytes": Path(p).stat().st_size} for p in inputs],
        "outputs": [{"path": str(p), "sha256": file_hash(p)} for p in outputs],
        "stats": stats,
    }
    (run_dir / "manifest.json").write_text(json.dumps(manifest, indent=2))
    return manifest

run_dir = Path("runs/first")
run_dir.mkdir(parents=True, exist_ok=True)
cfg = {"source": SOURCE, "output_dir": run_dir, "target_crs": "EPSG:27700", "min_area_m2": 50, "sample_n": 500, "seed": 42}
stats = run_pipeline(cfg)
manifest = write_manifest(run_dir, cfg, [SOURCE], [stats["output"]], stats)
print(json.dumps({k: manifest[k] for k in ("code", "inputs", "stats")}, indent=1)[:900])

{
 "code": {
  "commit": "3238b4f",
  "entrypoint": "pipeline.py"
 },
 "inputs": [
  {
   "path": "data/raw/inspire-rutland/2023-04-02/Rutland.zip",
   "sha256": "b33fece30cea1590a969cc9d6662a9264e4ddca81db3d0e17870901b2e022e6f",
   "bytes": 5911638
  }
 ],
 "stats": {
  "features_in": 21253,
  "features_out": 19470,
  "sample": 500,
  "sample_area_m2": 7919884.9,
  "output": "runs/first/parcels_sample.parquet"
 }
}

Pin, freeze the inputs, version, externalise, seed

Pin the environment, do not just record it

Recording versions tells you what happened; pinning them lets you get it back. A lock file captures the whole resolved dependency tree, not just the packages you named.

# uv — fast, lock-file first
uv lock                       # writes uv.lock with every transitive pin
uv sync --frozen              # recreate the exact environment anywhere

# pip-tools, or uv's drop-in for it
uv pip compile --generate-hashes requirements.in -o requirements.txt

# conda / mamba
conda env export --no-builds > environment.yml

For GIS specifically, GDAL is the version that matters most and the one people forget, because it arrives inside another package rather than as a direct dependency. It is the layer that actually reads your files. Ask the package that GeoPandas reads through — not gdalinfo on the command line, which can be a different GDAL entirely, as example 4 shows:

print(json.dumps(environment_report(), indent=1))

{
 "python": "3.14.4",
 "platform": "Linux-7.0.0-31-generic-x86_64-with-glibc2.43",
 "geopandas": "1.1.4",
 "pandas": "3.0.6",
 "shapely": "2.1.2",
 "geos": "3.13.1",
 "pyproj": "3.8.0",
 "proj": "9.8.1",
 "pyogrio": "0.13.0",
 "gdal": "3.12.4",
 "proj_network": false
}

PROJ deserves the same attention, and so does one setting in it. With network access on, PROJ downloads transformation grids it does not have; with it off, it falls back to a less accurate method. Here is the same British National Grid point transformed both ways:

from pyproj import Geod, Transformer

x, y = parcels_centroid                                       # the middle of Rutland, EPSG:27700
pyproj.network.set_network_enabled(False)
offline = Transformer.from_crs(27700, 4326, always_xy=True).transform(x, y)
pyproj.network.set_network_enabled(True)
online = Transformer.from_crs(27700, 4326, always_xy=True).transform(x, y)
pyproj.network.set_network_enabled(False)

_, _, shift = Geod(ellps="WGS84").inv(*offline, *online)
print(f"network off: {offline[0]:.7f}, {offline[1]:.7f}")
print(f"network on:  {online[0]:.7f}, {online[1]:.7f}   -> {shift:.2f} m apart")

network off: -0.6534383, 52.6617516
network on:  -0.6534269, 52.6617700   -> 2.19 m apart

Same code, same library versions, two answers about two metres apart: with the network on, PROJ fetched the OSTN15 grid from its CDN. Whether that happens depends on PROJ_NETWORK, on a config file, or on whether someone installed the grid package — none of which appears in a requirements file. It is why the environment report above records proj_network, and why a pipeline that must reproduce should either ship the grid or pin the setting.

Treat input data as immutable

The most common cause of an irreproducible GIS result is an input that changed under a name that did not. Adopt one rule: downloaded data is never modified and never overwritten. New version, new folder, dated.

data/
  raw/
    inspire-rutland/2023-04-02/Rutland.zip   ← never edited, never replaced
  interim/                                   ← regenerable, safe to delete
  processed/

The date in that folder name is not the download date. The setup below takes it from the server's Last-Modified header and saves the headers beside the file:

print(SOURCE)
print(json.loads(Path(SOURCE).with_name("source.json").read_text()))
print("sha256", file_hash(SOURCE)[:16], "…")

data/raw/inspire-rutland/2023-04-02/Rutland.zip
{'url': 'https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip', 'last_modified': 'Sun, 02 Apr 2023 05:28:03 GMT', 'etag': '"da0e4e809e74d225342c663639457097"', 'fetched': '2026-09-30'}
sha256 b33fece30cea1590 …

The URL — …/download/Rutland.zip — carries no version at all. The only records of which release you got are the header, the ETag and your own hash. Rutland's file says it was last modified in April 2023; a run next year may get the same bytes or new ones, and only the checksum in the manifest will tell you which.

The setup’s download, step by step: the release date comes from the server, not the calendar.

Put the code under version control and record the commit

A commit hash is the most compact possible description of "the code that ran". Recording it costs one subprocess call and makes git checkout <hash> a complete restoration of the logic. Flagging a dirty working tree matters as much as the hash:

print("clean:", git_commit())
Path("pipeline.py").write_text(Path("pipeline.py").read_text().replace('"min_area_m2"', '"min_area_m2"  '))
print("after an edit:", git_commit())
subprocess.run(["git", "checkout", "--", "pipeline.py"], check=True)
print("restored:", git_commit())

clean: 3238b4f
after an edit: 3238b4f-dirty
restored: 3238b4f

A run recorded as …-dirty is honest about being unreproducible, which is the second-best outcome and infinitely better than a hash that misrepresents what actually executed. The check ignores untracked files: a new output file in the folder is not a change to the code.

Externalise every parameter

A threshold living in a script is a parameter you will overwrite the next time you change it. In a config file it is a versioned, diffable artefact — and a copy of the resolved config goes into the run folder.

def snapshot_config(cfg, run_dir):
    resolved = json.loads(json.dumps(cfg, default=str))
    resolved["source"] = str(Path(resolved["source"]).resolve().relative_to(Path.cwd()))
    (run_dir / "config_used.json").write_text(json.dumps(resolved, indent=2))
    return resolved

print(snapshot_config(cfg, run_dir))

{'source': 'data/raw/inspire-rutland/2023-04-02/Rutland.zip', 'output_dir': 'runs/first', 'target_crs': 'EPSG:27700', 'min_area_m2': 50, 'sample_n': 500, 'seed': 42}

Resolved matters: after merging defaults, environment overlays, and command-line flags. That is what ran, and it is often not what any single file on disk says.

Fix every source of randomness

Randomness in GIS work is sneakier than in general data science because much of it is not obviously random.

kept = parcels[parcels.area >= 50]
unseeded = [set(kept.sample(n=500)["INSPIREID"]) for _ in range(2)]
seeded = [set(kept.sample(n=500, random_state=42)["INSPIREID"]) for _ in range(2)]
print("two unseeded samples identical:", unseeded[0] == unseeded[1])
print("two seeded samples identical:  ", seeded[0] == seeded[1])

print("seeded sample area:", round(float(kept.sample(n=500, random_state=42).area.sum()), 1), "m²")

two unseeded samples identical: False
two seeded samples identical:   True
seeded sample area: 7919884.4 m²

random.seed() and np.random.seed() do not reach pandas: sample() needs its own random_state, and so does anything in scikit-learn. Beyond explicit sampling, watch for these:

Row and file order. A join that produces ties resolves them by order, so unordered input gives different output. Path.glob() returns entries in the order the filesystem stores them — always sorted().
Set iteration. Sets have no order you should rely on; building a layer list from a set can shuffle results between interpreters.
Parallel workers. as_completed() and imap_unordered() return in completion order. Sort by a stable key before writing — see parallel batch processing.
Floating-point summation order. Summing in a different order can change the last digits of a total, though on this page's parcel areas it did not (example 3).

Write a manifest with every run

The manifest is the artefact that lets a future colleague — usually you — trust the output.

The manifest records output hashes, which raises a question worth testing: does writing the same data twice give the same bytes?

sample = gpd.read_parquet(stats["output"])
for suffix in ("gpkg", "parquet", "geojson", "fgb", "shp"):
    digests = []
    for i in range(2):
        target = Path(f"data/interim/same.{suffix}")                # the same path both times
        sample.to_parquet(target) if suffix == "parquet" else sample.to_file(target, layer="parcels")
        digests.append(file_hash(target)[:12])                     # for a shapefile: the .shp itself
        if suffix == "gpkg":
            shutil.copyfile(target, f"data/interim/same_{i}.gpkg")
        time.sleep(1.1)
    print(f"{suffix:8} {digests[0]} {digests[1]}  identical: {digests[0] == digests[1]}")

import sqlite3
stamps = [sqlite3.connect(f"data/interim/same_{i}.gpkg").execute("SELECT last_change FROM gpkg_contents").fetchone()[0]
          for i in range(2)]
print("gpkg_contents.last_change differs between the two:", stamps[0] != stamps[1])

def content_hash(frame: gpd.GeoDataFrame, key: str) -> str:
    frame = frame.sort_values(key)
    h = hashlib.sha256(frame.drop(columns="geometry").to_csv(index=False).encode())
    h.update(b"".join(frame.geometry.normalize().to_wkb()))
    return h.hexdigest()

print("content hash, both GeoPackages:",
      {content_hash(gpd.read_file(f"data/interim/same_{i}.gpkg"), "INSPIREID")[:12] for i in range(2)})

gpkg     310f91cc317f 407fc063cd8e  identical: False
parquet  aa595d2f8d25 aa595d2f8d25  identical: True
geojson  2f9ec5390155 2f9ec5390155  identical: True
fgb      4e7225d5877e 4e7225d5877e  identical: True
shp      6c8649575c6c 6c8649575c6c  identical: True
gpkg_contents.last_change differs between the two: True
content hash, both GeoPackages: {'f632156d32c7'}

Identical data written to the same path twice: four formats gave the same bytes, and the GeoPackage did not, because it stores a last-change timestamp in its gpkg_contents table. A byte hash of a GeoPackage output changes on every run, so a drift check built on it cries wolf. The shapefile's .shp was stable, but its .dbf header carries the date it was written, so the same test a day apart would fail there too. Hash outputs by content — sorted attributes plus normalised geometry, which gave one hash for both GeoPackages — or compare the stats, and keep byte hashes for inputs, which you never rewrite.

Make reruns easy and safe

Reproducibility that requires archaeology does not get used. One command should recreate the environment and rerun a specific configuration.

.PHONY: setup run verify
setup:
	uv sync --frozen

run:
	uv run python run_pipeline.py --config config/prod.yml

verify:
	uv run python tools/verify_manifest.py runs/first/manifest.json

The verify step asks whether the recorded inputs still exist, unchanged:

def verify_manifest(path):
    manifest = json.loads(Path(path).read_text())
    problems = []
    for entry in manifest["inputs"]:
        source = Path(entry["path"])
        if not source.exists():
            problems.append(f"missing input: {source}")
        elif file_hash(source) != entry["sha256"]:
            problems.append(f"input changed since the run: {source}")
    if manifest["code"]["commit"].endswith("-dirty"):
        problems.append("code had uncommitted changes when this ran")
    return problems or ["manifest verified: inputs unchanged"]

print(verify_manifest(run_dir / "manifest.json"))

# the same check against a copy that someone "fixed" in place
copy = Path("data/interim/Rutland.zip")
shutil.copyfile(SOURCE, copy)                  # copyfile: not the read-only mode
write_manifest(Path("data/interim"), cfg, [copy], [], {})
with copy.open("ab") as fh:
    fh.write(b"\n")
print(verify_manifest("data/interim/manifest.json"))

['manifest verified: inputs unchanged']
['input changed since the run: data/interim/Rutland.zip']

One appended byte, and the check says so. Without the manifest the next run would just have produced different numbers.

A self-describing run on Rutland's parcels

The data used on this page

HM Land Registry's INSPIRE Index Polygons for Rutland (© Crown copyright and database rights; OGL), downloaded into a folder named after the server's Last-Modified date. The setup also writes the small pipeline the page runs, pipeline.py, and commits it to a git repository in the working folder so the commit and dirty checks have something real to report.

import hashlib
import json
import os
import platform
import re
import shutil
import subprocess
import sys
import time
from datetime import datetime, timezone
from email.utils import parsedate_to_datetime
from pathlib import Path

import geopandas as gpd
import numpy as np
import pyproj
import requests

URL = "https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip"
found = sorted(Path("data/raw/inspire-rutland").glob("*/Rutland.zip"))
if found:
    SOURCE = found[-1]
else:
    response = requests.get(URL, timeout=300)
    released = parsedate_to_datetime(response.headers["Last-Modified"]).date().isoformat()
    SOURCE = Path("data/raw/inspire-rutland") / released / "Rutland.zip"
    SOURCE.parent.mkdir(parents=True)
    SOURCE.write_bytes(response.content)
    SOURCE.with_name("source.json").write_text(json.dumps(
        {"url": URL, "last_modified": response.headers["Last-Modified"], "etag": response.headers.get("ETag"),
         "fetched": datetime.now(timezone.utc).date().isoformat()}))
    SOURCE.chmod(0o444)                                    # read-only: raw data is never edited
Path("data/interim").mkdir(parents=True, exist_ok=True)

Path("pipeline.py").write_text('''from pathlib import Path
import geopandas as gpd

def run_pipeline(cfg):
    parcels = gpd.read_file(f"/vsizip/{cfg['source']}/Land_Registry_Cadastral_Parcels.gml",
                            columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO").to_crs(cfg["target_crs"])
    parcels["area_m2"] = parcels.area.round(1)
    kept = parcels[parcels["area_m2"] >= cfg["min_area_m2"]]
    sample = kept.sample(n=cfg["sample_n"], random_state=cfg["seed"]).sort_values("INSPIREID")
    out = Path(cfg["output_dir"]) / "parcels_sample.parquet"
    sample.reset_index(drop=True).to_parquet(out)
    return {"features_in": len(parcels), "features_out": len(kept), "sample": len(sample),
            "sample_area_m2": round(float(sample["area_m2"].sum()), 1), "output": str(out)}
''')
if not Path(".git").exists():
    for command in (["git", "init", "-q"], ["git", "add", "pipeline.py"],
                    ["git", "-c", "user.name=guide", "-c", "user.email=guide@example.org", "commit", "-qm", "pipeline"]):
        subprocess.run(command, check=True, env={**os.environ, "GIT_AUTHOR_DATE": "2026-09-30T12:00:00Z",
                                                  "GIT_COMMITTER_DATE": "2026-09-30T12:00:00Z"})

from pipeline import run_pipeline

parcels = gpd.read_file(f"/vsizip/{SOURCE}/Land_Registry_Cadastral_Parcels.gml", columns=["INSPIREID"], DOWNLOAD_SCHEMA="NO")
parcels_centroid = tuple(parcels.union_all().centroid.coords[0])
print(len(parcels), "parcels in", SOURCE.parent.name, "| CRS", parcels.crs.to_epsg())

21253 parcels in 2023-04-02 | CRS 27700

Example 1: A reproducibility preamble for any script

A few lines at the top of the entry point that pin randomness and record the environment.

import os, random

SEED = int(os.environ.get("SEED", 42))
random.seed(SEED)
np.random.seed(SEED)
rng = np.random.default_rng(SEED)                  # new code: pass a Generator around instead
print(random.random(), np.random.random(), rng.random())

def child_hash(**env):
    return subprocess.run([sys.executable, "-c", "print(hash('rutland'))"], capture_output=True, text=True,
                          env={**os.environ, **env}).stdout.strip()

print("two fresh interpreters, no PYTHONHASHSEED, same hash:", child_hash() == child_hash())
print("two fresh interpreters, PYTHONHASHSEED=42, same hash:", child_hash(PYTHONHASHSEED="42") == child_hash(PYTHONHASHSEED="42"))

0.6394267984578837 0.3745401188473625 0.7739560485559633
two fresh interpreters, no PYTHONHASHSEED, same hash: False
two fresh interpreters, PYTHONHASHSEED=42, same hash: True

String hashing is seeded when the interpreter starts, so hash() of the same string differs between runs unless PYTHONHASHSEED is set in the environment that launches Python — setting it inside a running script does nothing. Set it in the Makefile or the container, not in the code.

Example 2: A dated, self-describing run folder

Everything about a run in one directory: the data, the config, the log, the manifest.

def make_run_dir(base: Path, name: str, stamp: str | None = None) -> Path:
    stamp = stamp or datetime.now(timezone.utc).strftime("%Y-%m-%dT%H%M%SZ")
    run_dir = base / f"{stamp}_{name}"
    run_dir.mkdir(parents=True, exist_ok=False)     # never reuse a run folder
    (run_dir / "logs").mkdir()
    return run_dir

first = make_run_dir(Path("runs"), "parcels")
print("created", re.sub(r"\d", "#", first.name))
try:
    make_run_dir(Path("runs"), "parcels", stamp=first.name.split("_")[0])   # a second run in the same second
except FileExistsError as exc:
    print("second run refused:", type(exc).__name__)

created ####-##-##T######Z_parcels
second run refused: FileExistsError

exist_ok=False is deliberate: a collision means two runs are about to share a folder, and failing is the correct response. A stamp to the second is not unique on its own; add a short random suffix or the process id if runs can start together.

Example 3: Deterministic ordering everywhere it matters

from concurrent.futures import ThreadPoolExecutor, as_completed
import math

areas = parcels.area.to_numpy()
print("sum in file order:     ", repr(float(np.sum(areas))))
print("sum in sorted order:   ", repr(float(np.sum(np.sort(areas)))))
print("sum, plain Python loop:", repr(sum(areas.tolist())))
print("math.fsum (exact):     ", repr(math.fsum(areas)))

ordered = parcels.sort_values("INSPIREID").reset_index(drop=True)
chunks = [ordered.iloc[i:i + 3000] for i in range(0, len(ordered), 3000)]
with ThreadPoolExecutor(4) as pool:
    futures = [pool.submit(lambda c: (int(c["INSPIREID"].iloc[0]), round(float(c.area.sum()), 1)), c) for c in chunks]
    results = sorted(f.result() for f in as_completed(futures))    # completion order, then a stable sort
print(results[:2])

sum in file order:      428985424.6522182
sum in sorted order:    428985424.6522182
sum, plain Python loop: 428985424.6522182
math.fsum (exact):      428985424.6522182
[(31117368, 48268116.5), (31515320, 14762811.4)]

Here all four sums agree to the last digit: 21,253 areas of similar size, summed in double precision, leave nothing for the order to change. Summation order does move the last digits when values span many orders of magnitude, so a regression test should still compare with a tolerance rather than == — but it was not a cause of drift on this data, and the parallel results only needed the sort.

Example 4: Pinning the environment with Docker

When the environment must survive an operating-system upgrade, containerise it. The GDAL images carry no pip, so the build brings in uv to install the hashed requirements:

Path("docker").mkdir(exist_ok=True)
Path("docker/requirements.in").write_text("geopandas==1.1.4\npyogrio==0.13.0\nshapely==2.1.2\npyproj==3.8.0\n")
subprocess.run(["uv", "pip", "compile", "--quiet", "--python-version", "3.12", "--generate-hashes",
                "docker/requirements.in", "-o", "docker/requirements.txt"], check=True)
Path("docker/Dockerfile").write_text("""FROM ghcr.io/osgeo/gdal:ubuntu-small-3.12.2
COPY --from=ghcr.io/astral-sh/uv:0.9.5 /uv /usr/local/bin/uv
WORKDIR /app
COPY requirements.txt .
RUN uv pip install --system --break-system-packages --require-hashes -r requirements.txt
""")
build = subprocess.run(["docker", "build", "-q", "-t", "parcels-pipeline:2026-09-30", "docker"],
                       capture_output=True, text=True)
print("build:", "ok" if build.returncode == 0 else build.stderr[-300:])
check = subprocess.run(["docker", "run", "--rm", "parcels-pipeline:2026-09-30", "python3", "-c",
                        "import pyogrio, osgeo.gdal as g; print('GeoPandas reads with GDAL', "
                        "pyogrio.__gdal_version_string__, '| the image has GDAL', g.__version__)"],
                       capture_output=True, text=True)
print(check.stdout.strip())

build: ok
GeoPandas reads with GDAL 3.12.4 | the image has GDAL 3.12.2

The container answers a question most people get wrong. The image is "GDAL 3.12.2", but the pyogrio wheel brings its own GDAL, and that is the one GeoPandas reads files with. Pinning the base image pins the command-line tools; pinning pyogrio is what pins GeoPandas' GDAL. Tagging the image with a date rather than latest is what makes it an archive rather than a moving target, and --require-hashes ensures even a re-published package version cannot change the build. A plain pip install in that image fails with exit code 127: there is no pip to run.

What example 4’s container reported: pin the wheel to pin the GDAL GeoPandas uses.

Example 5: A regression test that catches silent drift

The strongest guarantee is an automated one: a small fixture through the real pipeline, with the numbers asserted.

Path("tests").mkdir(exist_ok=True)
Path("tests/test_pipeline_regression.py").write_text('''import math
import geopandas as gpd
from pipeline import run_pipeline

def test_known_input_gives_known_output(tmp_path):
    cfg = dict(source="SOURCE", output_dir=tmp_path, target_crs="EPSG:27700",
               min_area_m2=50, sample_n=500, seed=42)
    stats = run_pipeline(cfg)
    assert stats["features_out"] == FEATURES_OUT
    result = gpd.read_parquet(stats["output"])
    assert result.crs.to_epsg() == 27700
    assert math.isclose(result["area_m2"].sum(), SAMPLE_AREA, rel_tol=1e-9)
'''.replace("SOURCE", str(SOURCE)).replace("FEATURES_OUT", str(stats["features_out"]))
   .replace("SAMPLE_AREA", str(stats["sample_area_m2"])))

test = subprocess.run([sys.executable, "-m", "pytest", "-q", "-p", "no:cacheprovider", "tests"],
                      capture_output=True, text=True)
print(re.sub(r" in [\d.]+s", "", test.stdout.strip().splitlines()[-1]))

1 passed

The expected numbers are the ones this run produced — 500 sampled parcels out of the ones over 50 m², with their total area. When a library upgrade changes a result, this test fails on the upgrade rather than six months later in a report. Loosen the tolerance rather than deleting the assertion when a change is legitimate — and record why in the commit message.

Reproducibility is five properties

Reproducibility is not one property but five, and a workflow is only as reproducible as its weakest layer. Pinning every library while reading an input that gets replaced behind the same URL buys nothing. Checksumming the data while the analysis lives in an uncommitted notebook buys nothing either. The five layers — environment, code, data, parameters, randomness — have to be addressed together, which is why the manifest that records all five in one file is the single highest-value habit in this article.

Every "it gave different numbers this time" traces back to something that was not pinned.

GIS has a reproducibility hazard that general data analysis does not: the transformation stack underneath your code is itself versioned and configurable. PROJ ships better datum-shift grids and fetches them when allowed to; GDAL fixes edge cases in drivers; a wheel can carry a GDAL other than the one on the machine. These are improvements, and they change results — two metres, on this page, from one setting. That is exactly why recording the GDAL and PROJ versions, and the network setting, matters more than recording your own library's: when a number moves, the first question is which layer moved underneath it, and without a manifest that question takes days instead of minutes.

There is a distinction worth being precise about. Repeatability is getting the same answer on the same machine tomorrow. Reproducibility is a colleague getting the same answer on a different machine. Replicability is a different team reaching the same conclusion with their own data and code. Manifests and seeds give you the first. Lock files, immutable inputs, and committed code give you the second. The third is a scientific question, not an engineering one — but you cannot even attempt it without the first two.

None of this needs to be adopted at once. In practice the order of return is clear: write a manifest (an afternoon, and it makes every later question answerable), then stop overwriting input data (a folder convention), then pin the environment with a lock file, then add a seed and deterministic ordering, and finally add a regression test on a small fixture. Each step is useful alone, and the first one is useful immediately — including for runs that happened before you adopted the rest.

Hashes, data you do not control and notebooks

Hashing outputs

The manifest section’s test: only the GeoPackage changed its bytes between two identical writes.

An output hash only helps if writing the same data gives the same bytes. It does not for a GeoPackage, and a shapefile's .dbf records the date it was written. Hash those by content, or record counts and totals and compare those.

Hashing very large inputs

SHA-256 over a 50 GB raster takes minutes. Hash a header plus a sample of blocks, or record size, modification time, and any internal metadata the format provides. Note in the manifest which method you used, so nobody later mistakes a partial fingerprint for a full one.

Data you do not control

For a WFS, a live database or a URL with no version in it, there is nothing stable to point at. Record what you can: the query, the timestamp of the fetch, the server's Last-Modified header and ETag, and a hash of the fetched result. Then save that fetched snapshot as an immutable input — it is now the reproducible artefact.

Floating-point differences across machines

Identical code and versions can still differ in the last bits across CPU architectures or BLAS builds. Compare with tolerances (pytest.approx, numpy.isclose) rather than exact equality, and reserve byte-for-byte comparison for a single controlled environment such as a container.

The page’s steps as a pre-report check.

Notebooks

Notebooks record output but not execution order, and a notebook that ran cells out of order is not reproducible even by its author. Keep the pipeline in modules under version control and use notebooks to call it, or enforce top-to-bottom execution in CI with jupyter nbconvert --execute.

"-dirty" runs

Sometimes you must run with uncommitted changes. Record it honestly rather than suppressing it, and treat a -dirty manifest as provisional: fine for an exploratory run, not acceptable for a number that goes into a report.

Storage cost of immutable data

Never overwriting inputs means storage grows. That is usually the right trade — disk is cheaper than an irreproducible result — but manage it deliberately: compress older snapshots, move them to cheaper storage, keep the checksums forever even when the bytes are archived offline.

Internal links

FAQ

What does it take to make a GIS workflow reproducible?

Pin and record five things: the environment (a lock file, including the GDAL and PROJ versions and PROJ's network setting), the code (a git commit), the inputs (immutable files with checksums), the parameters (a config file saved with the run), and randomness (a fixed seed plus deterministic ordering). Missing any one of them is enough to break "run it again".

Why do my results change after upgrading GeoPandas or GDAL?

Because the layer underneath your code changed. GDAL fixes driver edge cases and PROJ ships improved datum-transformation grids — and downloads them when network access is on, which moved a British National Grid point about two metres in the example on this page. These are genuine improvements that change output. Recording those versions in a manifest turns a baffling difference into a one-line explanation.

What should a run manifest contain?

The UTC timestamp, the git commit (flagged if the tree was dirty), the Python and library versions including GDAL and PROJ, the checksum and size of every input, the fully resolved config, the checksum of every output (by content for a GeoPackage, which changes its bytes on every write), and summary statistics such as feature counts and runtime. Write it next to the output, not in a separate log.

How do I keep input data immutable?

Never write into the folder you downloaded to, and never overwrite a raw file. Save each acquisition under a dated folder, treat those folders as read-only, and checksum the files you read so the manifest records exactly which version a run used. Regenerable intermediates go somewhere clearly disposable.

Is `random.seed()` enough to make a run deterministic?

No. NumPy and pandas have their own generators (np.random.seed, random_state=), and much GIS non-determinism is not "random" at all — unsorted glob() results, set iteration order, and parallel workers returning out of order. Sort file lists, sort rows before writing, and impose a stable order on parallel results.

Do I need Docker for reproducible GIS work?

Not to start. A lock file, immutable inputs, a committed script, and a manifest cover most needs and cost far less. Reach for a container when the environment must survive an OS upgrade, when GDAL is difficult to install consistently across machines, or when someone else must reproduce your work on a system you cannot see.

How do I know reproducibility has broken before it matters?

Add a regression test: run a small committed fixture through the real pipeline in CI and assert on feature counts and a summed measurement. When a dependency upgrade changes a result, the test fails on the upgrade instead of six months later in a report — and you get to decide deliberately whether the new number is the better one.