Reproducible Spatial Workflows

4.8. Reproducible Spatial Workflows#

4.8.1. Learning objectives#

By the end of this section, you should be able to:

  • explain what can change between two runs of the same analysis and how to control each cause;

  • separate configuration from analysis;

  • preserve raw inputs and keep raw, intermediate, and final outputs apart;

  • choose an export format and explain its trade-offs;

  • export results with a provenance record and a methods note;

  • test sensitivity to thresholds before drawing a conclusion; and

  • write an analysis plan that another person can review and rerun.

Section goal

Turn the exploratory skills of S01–S07 into a result that another person can regenerate and check. A small synthetic layer shows each practice working: configuration, raw input, export, provenance, and a sensitivity sweep.

4.8.2. From exploration to an auditable result#

Interactive notebooks are useful for discovery, but a defensible result should be regenerable by another person [Sandve et al., 2013, University of Southern California et al., 2021]. Keep the source, decisions, transformations, checks, and outputs visible rather than relying on hidden notebook state.

Two runs of the same analysis can differ for several reasons, and each has its own remedy. Table 4.29 lists them.

Table 4.29 What can change between runs, and the practice that controls it.#

Cause

What goes wrong

Control

Hidden notebook state

Cells run out of order, or variables survive from deleted cells

Restart the kernel and run all cells from the top

Changing data

Open-data services update or remove features

Cache raw responses, record the query and date, store a checksum

Software versions

Libraries change CRS, predicate, or overlay behaviour

Record the Python and library versions

Machine-specific paths

A personal directory, or a source file that gets overwritten

Relative paths, parameters, and a read-only raw folder

Randomness

Samples that differ on every run

Fix the seed, or sample deterministically

Unstated choices

Thresholds and predicates buried in code

Name them as parameters at the top

4.8.3. A reusable workflow#

  1. Define the question and unit of analysis. Decide whether the unit is a facility, polygon, pixel, edge, or route.

  2. Record provenance. Save the provider, URL or query, retrieval date, license, data vintage, and software versions.

  3. Preserve raw inputs. Do not overwrite the source while cleaning or reprojecting.

  4. Inspect and validate. Check geometry types, CRS, missing values, identifiers, bounds, and expected ranges.

  5. Transform deliberately. Use a documented CRS, predicate, overlay, threshold, or raster mask.

  6. Validate intermediate results. Check row counts, area or length totals, and the number of unmatched records.

  7. Export data and methods. Store results with metadata and a short methods note.

  8. Review sensitivity. Test thresholds, boundary rules, source dates, and alternative assumptions that could change the conclusion.

4.8.4. Separate configuration from analysis#

Stable configuration makes examples easier to rerun. Put every choice that could change the answer, such as the study area, data date, CRS, thresholds, and output location, in one place at the top, and have later cells read it. A threshold buried in the middle of an analysis is easy to miss and hard to vary, whereas a named constant can be reviewed and swept, as the sensitivity section below does:

Hide code cell source

 
from pathlib import Path

STUDY_AREA = "Calgary, Alberta, Canada"
SOURCE_DATE = "2026-07-26"
ANALYSIS_CRS = "EPSG:32611"        # metric CRS used for every measurement
SERVICE_RADIUS_M = 500             # threshold, in metres
OUTPUT_DIR = Path("output/chapter_04")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

Keep paths relative to the project or pass them as parameters. Avoid embedding a personal home directory, writing over source files, or depending on variables created in an earlier interactive cell.

For downloaded data, cache raw responses when the license permits it and record the query parameters. Open-data services change over time, so a rerun may not return the same features even when the code is unchanged. A checksum, a short fingerprint computed from the bytes of the raw file, proves later that two runs used identical input.

4.8.5. Export spatial results with metadata#

GeoPackage [Open Geospatial Consortium, 2023] is a useful default for derived vector layers because it keeps geometry and attributes together. Table 4.30 compares it with the alternatives.

Table 4.30 Common export formats and their trade-offs for spatial data.#

Format

Strengths

Limitations

GeoPackage (.gpkg)

One file, several layers, stored CRS, long column names

Less common in web tools than GeoJSON

GeoJSON

Plain text and web friendly

Expects longitude and latitude, and is bulky for large layers

Shapefile

Read by almost every tool

Several files per layer, column names cut to 10 characters, one geometry type, 2 GB limit

GeoParquet

Compact and fast for large tables

Needs recent tools

CSV with WKT

Simplest to inspect

No CRS, no data types, and geometry stored as text

GeoTIFF

Standard for rasters, with CRS, transform, and no-data stored

Cells only, with no vector features

The cell below runs the whole pattern on a synthetic layer. A stand-in “download” is saved as the raw input and fingerprinted. The analysis then reads it, measures areas in a metric CRS, and exports a result that keeps the original CRS, because delivery often requires it. Raw and final outputs live in separate folders.

Hide code cell source

 
import hashlib

import geopandas as gpd
from shapely.geometry import box

RAW_DIR, FINAL_DIR = OUTPUT_DIR / "raw", OUTPUT_DIR / "final"
for folder in (RAW_DIR, FINAL_DIR):
    folder.mkdir(exist_ok=True)

# A stand-in for a downloaded response: three synthetic parks in longitude and latitude.
downloaded = gpd.GeoDataFrame(
    {
        "asset_cd": ["P-001", "P-002", "P-003"],
        "site_name": ["Central Green", "Riverside", "Bowness Park"],
    },
    geometry=[
        box(705500, 5659000, 705600, 5659100),    # 1 ha
        box(705700, 5659000, 705900, 5659100),    # 2 ha
        box(706000, 5659000, 706100, 5659300),    # 3 ha
    ],
    crs=ANALYSIS_CRS,
).to_crs("EPSG:4326")

raw_path = RAW_DIR / "parks_raw.geojson"
downloaded.to_file(raw_path, driver="GeoJSON")
raw_checksum = hashlib.sha256(raw_path.read_bytes()).hexdigest()

# Analysis: measure in the metric CRS, deliver in the source CRS.
raw = gpd.read_file(raw_path)
result = raw.assign(area_ha=raw.to_crs(ANALYSIS_CRS).area / 10_000)

assert len(result) == len(raw)
assert result.geometry.notna().all()
assert (result["area_ha"] >= 0).all()

result.to_file(FINAL_DIR / "parks.gpkg", layer="parks", driver="GPKG")

# Read the export back to confirm that it round-trips.
check = gpd.read_file(FINAL_DIR / "parks.gpkg", layer="parks")
assert check.crs == result.crs and len(check) == len(result)
assert check.geometry.geom_equals(result.geometry).all()

print(result.drop(columns="geometry").round(3).to_string(index=False))
print("raw input fingerprint:", raw_checksum[:12])

asset_cd     site_name  area_ha
   P-001 Central Green      1.0
   P-002     Riverside      2.0
   P-003  Bowness Park      3.0
raw input fingerprint: 6bdc2dda186b

The export round-trips: three rows in and three out, with the CRS and geometry unchanged, and the areas come back as exactly 1, 2, and 3 ha even though the stored coordinates are degrees, because the measuring was done in the metric CRS. The 12-character fingerprint is the start of the full SHA-256 hash that the next step stores.

The export is only half of the deliverable. A provenance record states where the data came from and how it was processed [Wilkinson et al., 2016], and a short methods note turns that record into text a reader can check. Software versions belong in the record because library behavior changes between releases; the versions printed below depend on the environment that ran the cell.

Hide code cell source

 
import importlib.metadata as package_metadata
import json
import platform

metadata = {
    "study_area": STUDY_AREA,
    "source": "synthetic park layer (stand-in for a downloaded response)",
    "source_date": SOURCE_DATE,
    "input_sha256": raw_checksum,
    "crs": str(result.crs),
    "analysis_crs": ANALYSIS_CRS,
    "rows": len(result),
    "geometry_types": result.geom_type.value_counts().to_dict(),
    "software": {
        "python": platform.python_version(),
        **{name: package_metadata.version(name) for name in ["geopandas", "shapely", "pyproj"]},
    },
}
(FINAL_DIR / "metadata.json").write_text(json.dumps(metadata, indent=2), encoding="utf-8")
print(json.dumps(metadata, indent=2))

{
  "study_area": "Calgary, Alberta, Canada",
  "source": "synthetic park layer (stand-in for a downloaded response)",
  "source_date": "2026-07-26",
  "input_sha256": "6bdc2dda186bce907f2d2f54478b82af786e85ce089769826324fb556ce7b62c",
  "crs": "EPSG:4326",
  "analysis_crs": "EPSG:32611",
  "rows": 3,
  "geometry_types": {
    "Polygon": 3
  },
  "software": {
    "python": "3.13.15",
    "geopandas": "1.2.0",
    "shapely": "2.2.0",
    "pyproj": "3.8.0"
  }
}

The record holds what a rerun must match: the input hash, CRS, row count, geometry types, and the versions of Python (3.13.15), GeoPandas (1.2.0), Shapely (2.2.0), and pyproj (3.8.0). If a later run gives different numbers, compare these fields first to see whether the input or the software changed.

Hide code cell source

 
methods_note = f'''# Methods note: park areas

- **Question:** How large is each park site, in hectares?
- **Unit of analysis:** one park site per row.
- **Source:** {metadata["source"]}, retrieved {metadata["source_date"]} (SHA-256 {raw_checksum[:12]}).
- **Processing:** areas measured in {ANALYSIS_CRS}; layer delivered in {metadata["crs"]}.
- **Checks:** {metadata["rows"]} rows in and out, no missing geometry, areas non-negative.
- **Limitation:** areas are planar and depend on the chosen CRS.
'''
(FINAL_DIR / "methods.md").write_text(methods_note, encoding="utf-8")
print(methods_note)

for path in sorted(OUTPUT_DIR.rglob("*")):
    if path.is_file():
        print(path.as_posix())

# Methods note: park areas

- **Question:** How large is each park site, in hectares?
- **Unit of analysis:** one park site per row.
- **Source:** synthetic park layer (stand-in for a downloaded response), retrieved 2026-07-26 (SHA-256 6bdc2dda186b).
- **Processing:** areas measured in EPSG:32611; layer delivered in EPSG:4326.
- **Checks:** 3 rows in and out, no missing geometry, areas non-negative.
- **Limitation:** areas are planar and depend on the chosen CRS.

output/chapter_04/.DS_Store
output/chapter_04/final/metadata.json
output/chapter_04/final/methods.md
output/chapter_04/final/parks.gpkg
output/chapter_04/metadata.json
output/chapter_04/parks.gpkg
output/chapter_04/raw/parks_raw.geojson

The listing shows raw input and final outputs in separate folders, with the metadata and the methods note beside the data they describe. The GeoPackage default also has a practical reason. The same layer written as a shapefile keeps its columns only if their names are 10 characters or fewer [ESRI, 1998, GDAL/OGR contributors, 2026]:

Hide code cell source

 
import tempfile

with tempfile.TemporaryDirectory() as folder:
    folder = Path(folder)
    wide = result.rename(columns={"site_name": "community_district_name"})
    wide.to_file(folder / "layer.shp")
    wide.to_file(folder / "layer.gpkg", layer="layer", driver="GPKG")

    shapefile_parts = sorted(p.suffix for p in folder.iterdir() if p.stem == "layer" and p.suffix != ".gpkg")
    print("shapefile parts   :", shapefile_parts)
    print("shapefile columns :", list(gpd.read_file(folder / "layer.shp").columns))
    print("GeoPackage columns:", list(gpd.read_file(folder / "layer.gpkg").columns))

shapefile parts   : ['.cpg', '.dbf', '.prj', '.shp', '.shx']
shapefile columns : ['asset_cd', 'community_', 'area_ha', 'geometry']
GeoPackage columns: ['asset_cd', 'community_district_name', 'area_ha', 'geometry']

/var/folders/fd/9ctr8dys3tv3jyfrxn6_8y5m0000gp/T/ipykernel_18008/1658088252.py:6: UserWarning: Column names longer than 10 characters will be truncated when saved to ESRI Shapefile.
  wide.to_file(folder / "layer.shp")
/Users/hdtg3/miniconda/lib/python3.13/site-packages/pyogrio/raw.py:733: RuntimeWarning: Normalized/laundered field name: 'community_district_name' to 'community_'
  ogr_write(

The shapefile needs five files and cuts community_district_name down to community_, whereas the GeoPackage is one file with the full name. GDAL warns about the truncation, but the warning is easy to miss.

The same export on the repository layer [Open Calgary: The City of Calgary Open Data Portal, 2026] follows. It reads a downloaded layer, writes it to a GeoPackage, and records the study area, date, CRS, row count, and geometry types:

Hide code cell source

 
import json
import geopandas as gpd

result = gpd.read_file("../data/Parks_Sites_20261019.geojson")
result.to_file(OUTPUT_DIR / "parks.gpkg", layer="parks", driver="GPKG")

metadata = {
    "study_area": STUDY_AREA,
    "source_date": SOURCE_DATE,
    "crs": str(result.crs),
    "rows": len(result),
    "geometry_types": result.geom_type.value_counts().to_dict(),
}
(OUTPUT_DIR / "metadata.json").write_text(
    json.dumps(metadata, indent=2),
    encoding="utf-8",
)

167

Keep raw, intermediate, and final outputs separate. If a result is reprojected for measurement, retain the original CRS for delivery when that is part of the data contract.

4.8.6. Make maps and tables reviewable#

Every result should state the study area, date, CRS, source, operation, threshold, and unit. A plausible map can still contain a degree-based distance error, an accidental one-to-many join, or a raster no-data value interpreted as zero.

Useful review outputs include:

  • a source and transformation summary;

  • input and output row counts;

  • geometry validity and missing-value counts;

  • unique identifiers before and after joins;

  • total area or length in metric units; and

  • a map of a representative sample and any unmatched features.

For graph-specific provenance, apply the same principles to Chapter 5.

4.8.7. Test sensitivity before concluding#

A threshold is a choice, so a conclusion is only as firm as its stability when the choice changes. The efficient pattern is to compute the underlying quantity once, here the distance from each district’s representative point to its nearest library, and then compare it with several thresholds. The toy example below has three libraries on a line and four districts, and SERVICE_RADIUS_M from the configuration is one of the tested radii.

Hide code cell source

 
import pandas as pd
from shapely.geometry import Point

libraries = gpd.GeoSeries([Point(0, 0), Point(1500, 0), Point(3000, 0)])
districts = gpd.GeoDataFrame(
    {"district": ["D1", "D2", "D3", "D4"]},
    geometry=[Point(200, 100), Point(1800, 0), Point(2600, 200), Point(4900, 0)],  # representative points
)
districts["nearest_m"] = districts.geometry.apply(lambda point: libraries.distance(point).min())
display(districts[["district", "nearest_m"]].round(0).set_index("district"))

rows = []
for radius in sorted({250, SERVICE_RADIUS_M, 1000, 2000}):
    covered = districts.loc[districts["nearest_m"] <= radius, "district"].tolist()
    rows.append({"radius_m": radius, "share_covered": f"{len(covered) / len(districts):.0%}",
                 "covered": ", ".join(covered)})
display(pd.DataFrame(rows).set_index("radius_m"))
nearest_m
district
D1 224.0
D2 300.0
D3 447.0
D4 1900.0
share_covered covered
radius_m
250 25% D1
500 75% D1, D2, D3
1000 75% D1, D2, D3
2000 100% D1, D2, D3, D4

At the 500 m default, three of four districts are covered, and that headline holds at 1,000 m. Full coverage appears only at 2,000 m because of the outlying district D4, and at 250 m only one district qualifies. A report that gave only the 500 m figure would hide that spread, so the sweep belongs in the methods note. The distances above also show how near the 500 m result is to a threshold: D3 is only 53 m inside the radius (447 m), so a 400 m radius would cover just half of the districts. The same approach applies to boundary rules (within compared with intersects), source dates, and resolution.

4.8.8. Capstone prompt#

Choose a local question such as “Which community districts contain a library?” Write a one-page analysis plan before coding. Include:

  1. the source layers and data types;

  2. the unit of analysis and expected geometry;

  3. the CRS and threshold units;

  4. the operation sequence;

  5. at least three validation checks;

  6. the output data and map; and

  7. one limitation and one sensitivity test.

Table 4.31 lists what each part must state, and a gap that each one commonly hides.

Table 4.31 What each part of the analysis plan must state.#

Plan item

It must state

A common gap

Source layers and data types

Provider, geometry type, and date of each layer

Not saying which layer is the left side of a join

Unit of analysis and expected geometry

What one output row describes, and its geometry type

Counting joined rows instead of features

CRS and threshold units

The storage CRS, the measurement CRS, and the units

A threshold in metres applied to degrees

Operation sequence

Each step in order, with its predicate or overlay

Leaving out repair and reprojection

Validation checks

Expected row counts, totals, and unmatched records

Checks that could never fail

Output data and map

Format, fields, and what the map shows

A map with no source, date, or units

Limitation and sensitivity test

One thing the data cannot show, and one parameter varied

Reporting a single threshold as the answer

4.8.9. Summary#

Reproducibility connects spatial reasoning to trustworthy data engineering. Preserve sources, make choices explicit, validate every stage, export useful representations with a provenance record and a methods note, test how conclusions change with their thresholds, and communicate what the data cannot establish.