Geospatial Data Ingestion & Processing Workflows

Every spatial analysis is only as trustworthy as the pipeline that fed it. This guide is the architectural map for building production-grade geospatial ETL in Python — the layer that sits between raw files on disk and the Spatial Analysis & Advanced Query Techniques you run downstream. It covers ingestion through Shapefile & GeoJSON Parsing, coordinate normalization via Coordinate Reference System Transformations, integrity enforcement with Topology Validation & Repair, attribute enrichment through Spatial Joins & Merging, and scalable storage in Cloud-Native Geospatial Formats. Where the input is a column of addresses rather than geometry at all, Geocoding & Address Pipelines covers the normalisation, rate limiting and caching that turn text into trustworthy coordinates. The primitives themselves — GeoDataFrames, geometry algebra, projection objects — are documented in Mastering Core Geospatial Python Libraries; here we focus on wiring them into a repeatable, memory-safe workflow for data scientists, GIS analysts, and urban planners.

Spatial data ingestion pipeline A five-stage pipeline: ingest raw formats, validate topology, transform to a target CRS, store as GeoParquet, and analyse — with each stage feeding the next. Ingest SHP / GeoJSON Validate topology Transform target CRS Store GeoParquet Analyse query From raw files to a queryable store Validate and reproject once, early — every downstream stage depends on clean, correctly-projected geometry
A robust pipeline validates and reprojects early, then stores in a columnar, cloud-native format ready for analysis.

Ecosystem Architecture & Dependency Management

The Python geospatial stack is a thin, well-designed set of wrappers over three C/C++ libraries: GDAL/OGR for format translation, PROJ for coordinate transformation, and GEOS for planar geometry. Almost every dependency conflict you will hit traces back to one of those three being present in two incompatible versions at once — for example a pyproj wheel bundling PROJ 9.x while a system GDAL links PROJ 8.x, so datum grids resolve differently depending on which library touches the geometry first. The single most reliable defense is environment isolation with pinned versions, installed from a channel that ships the binaries together.

For reproducible installs, prefer conda-forge (or mamba) where the C libraries are compiled as a coherent set:

mamba create -n geo-etl python=3.11 \
  "geopandas>=1.0" "shapely>=2.0" "pyproj>=3.6" \
  "rasterio>=1.3" "pyogrio>=0.7" "dask-geopandas>=0.4" \
  "pyarrow>=15" -c conda-forge

If you must use pip, install pyogrio and rasterio from wheels that vendor their own GDAL/PROJ, and never mix them with a apt/brew GDAL on the same PYTHONPATH. Once installed, run a verification script on first import in CI — it fails loudly instead of producing silently wrong coordinates months later:

"""Verify the geospatial toolchain is internally consistent. Run in CI."""
import geopandas as gpd
import shapely
import pyproj
import rasterio
from rasterio._version import gdal_version

print(f"geopandas : {gpd.__version__}")
print(f"shapely   : {shapely.__version__} (GEOS {shapely.geos_version_string})")
print(f"pyproj    : {pyproj.__version__} (PROJ {pyproj.proj_version_str})")
print(f"rasterio  : {rasterio.__version__} (GDAL {gdal_version()})")

# Fail fast on the classic mismatch: a CRS that PROJ cannot round-trip
crs = pyproj.CRS.from_epsg(25832)          # ETRS89 / UTM 32N
assert crs.to_epsg() == 25832, "PROJ cannot resolve EPSG:25832 — grid files missing"
assert shapely.geos_version >= (3, 10, 0), "GEOS too old for make_valid()"

Pin these versions in a lockfile (conda-lock, pip-tools, or uv) and rebuild the environment from the lock — not from loose ranges — so that a pyproj datum-grid update never changes analytical results between two runs of the same pipeline. Treat the toolchain like any other production dependency: version it, test it, and gate deployment on the verification script above.

Three version boundaries in this stack change behaviour rather than just fixing bugs, and a pipeline written on one side of them will not behave identically on the other. Knowing which side you are on is worth more than knowing the exact patch numbers.

The practical consequence for an ingestion pipeline is that "it worked last quarter" is not evidence the environment is sound. Record the resolved versions in the run's own output — a small metadata sidecar written next to every artefact — so that when a number changes six months later you can tell whether the data moved or the stack did.

import json
import geopandas as gpd
import shapely
import pyproj
import rasterio
from rasterio._version import gdal_version

def toolchain_fingerprint() -> dict:
    """Emit alongside every artefact so results stay attributable to a stack."""
    return {
        "geopandas": gpd.__version__,
        "shapely": shapely.__version__,
        "geos": shapely.geos_version_string,
        "pyproj": pyproj.__version__,
        "proj": pyproj.proj_version_str,
        "rasterio": rasterio.__version__,
        "gdal": gdal_version(),
        "proj_data_dir": pyproj.datadir.get_data_dir(),
    }

with open("stage/_toolchain.json", "w") as fh:
    json.dump(toolchain_fingerprint(), fh, indent=2)

Core Concepts & Data Model

The mental model for spatial ETL in Python is a single table where one column happens to hold geometry. That table is the GeoDataFrame — a pandas DataFrame with an active geometry column and an attached CRS — and understanding it end to end unlocks the entire workflow; the full object model is documented in GeoPandas DataFrames Explained. Each geometry cell is a Shapely object (Point, LineString, Polygon, or their multi-part variants), and the whole column is backed by a contiguous GEOS array in Shapely 2.0, which is what makes vectorized operations fast.

Three properties travel with the data and must never be treated as decoration:

import geopandas as gpd
from shapely.geometry import Point

# A GeoDataFrame is a table with one geometry column and one CRS
parcels = gpd.GeoDataFrame(
    {"parcel_id": ["A-101", "A-102"], "area_use": ["residential", "commercial"]},
    geometry=[Point(11.34, 44.49), Point(11.35, 44.50)],
    crs="EPSG:4326",                       # WGS84 lon/lat — degrees, not metres
)

print(parcels.geometry.name)               # -> 'geometry' (the active column)
print(parcels.crs.axis_info[0].unit_name)  # -> 'degree'
print(parcels.geom_type.unique())          # -> ['Point']

A geometry cell has four states, not two, and confusing them is behind a large share of pipeline defects. A geometry can be valid, invalid (a ring that crosses itself — still a geometry, and still has coordinates), empty (a real geometry object containing no coordinates, so is_empty is true and area is zero), or missing (None, no object at all). Each is detected by a different accessor, and most operations treat them differently: is_valid returns False for missing geometries as well as invalid ones, area returns zero for both empty and missing, and a spatial join drops empties and missings alike without comment. A validation gate that checks only is_valid therefore passes a frame that is half empty.

import geopandas as gpd
from shapely.geometry import Polygon
from shapely import wkt

audit = gpd.GeoDataFrame(
    {"parcel_id": ["A-101", "A-102", "A-103", "A-104"]},
    geometry=[
        Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),          # valid
        wkt.loads("POLYGON((0 0, 1 1, 1 0, 0 1, 0 0))"),    # invalid: bow-tie
        wkt.loads("POLYGON EMPTY"),                          # empty
        None,                                                # missing
    ],
    crs="EPSG:25832",
)

print(audit.geometry.is_valid.tolist())    # [True, False, True, False]
print(audit.geometry.is_empty.tolist())    # [False, False, True, False]
print(audit.geometry.isna().tolist())      # [False, False, False, True]

The attribute half of the frame carries its own trap: index alignment. Assigning a GeoSeries or a Series into a column aligns on the index, not on position. After an explode() the index repeats, after a filtered sjoin it has gaps, and an assignment from a differently-indexed result then fills the mismatched rows with NaN rather than raising. The habit that prevents it is to reset_index(drop=True) immediately after any operation that changes row cardinality, and to assign .values when you genuinely mean positional assignment.

One more structural detail worth internalising early: the CRS lives on the GeoSeries, not on the frame. A GeoDataFrame may hold several geometry columns — a parcel boundary, its centroid, a mailing point — each with its own CRS, and gdf.crs reports only the active one. set_geometry("centroid") changes which column every geometry operation applies to, and silently changes which CRS gdf.crs is talking about with it.

Rasters follow a parallel model. A Rasterio dataset is an N-band array plus a geotransform (an affine mapping array indices to world coordinates) and a CRS. The bridge between the two worlds — burning vectors into a raster grid or extracting pixel statistics per polygon — is where most real pipelines earn their keep, so keep both data models straight: vectors carry discrete geometry, rasters carry a regular grid, and both carry a CRS that must agree before they can be combined.

Vector and raster data models sharing one CRS Two parallel data models — a GeoDataFrame holding a vector table with one active geometry column, and a Rasterio dataset holding an N-band grid with an affine geotransform — both carry a coordinate reference system that must agree before the two can be combined. Two data models, one shared CRS GeoDataFrame — vector table parcel_id area_use geometry A-101 residential POLYGON… A-102 commercial POLYGON… one active geometry column + one CRS Rasterio dataset — raster grid N-band array + affine geotransform Shared CRS both must agree before combining
A GeoDataFrame stores discrete vector geometry in one active column; a Rasterio dataset stores a regular N-band grid. Both must carry the same CRS before you can burn vectors into the grid or extract pixel statistics per polygon.

Key Operations & Vectorized Workflows

Five operations carry the vast majority of ingestion pipelines. Each is vectorized — express it once against the whole column and let GEOS iterate in C, never in a Python for loop over rows.

1. Streamed reading. Large national or municipal datasets should never be loaded whole. Read in bounded chunks so peak memory stays flat regardless of file size:

import geopandas as gpd
import pyogrio

def ingest_stream(path: str, chunk_size: int = 50_000):
    """Yield GeoDataFrame chunks so memory stays bounded on multi-GB inputs."""
    info = pyogrio.read_info(path)
    for start in range(0, info["features"], chunk_size):
        yield gpd.read_file(
            path, engine="pyogrio",
            rows=slice(start, start + chunk_size),
        )

for batch in ingest_stream("parcels_national.gpkg"):
    print(f"{len(batch)} features | CRS: {batch.crs.to_epsg()}")

2. Filtering and column pruning. Drop unused attribute columns and clip to a bounding box before any heavy geometry work — narrow tables and fewer features cut both memory and compute:

bbox = (11.20, 44.40, 11.50, 44.60)                     # (minx, miny, maxx, maxy)
parcels = gpd.read_file("parcels_national.gpkg", bbox=bbox, engine="pyogrio")
parcels = parcels[["parcel_id", "area_use", "geometry"]]  # prune before processing

Push both filters into the reader rather than applying them afterwards. bbox= becomes a spatial-index lookup in formats that carry one — FlatGeobuf's packed Hilbert R-tree, GeoPackage's RTree table — so the features outside the box are never parsed at all, and columns= stops the driver from decoding attribute values you are about to discard. Selecting columns after the read costs the full parse and saves only memory; selecting them in the read saves both. The gap widens sharply on GeoParquet, where column pruning is a property of the file layout rather than a filter applied to parsed rows:

parcels = gpd.read_parquet(
    "stage/parcels.parquet",
    columns=["parcel_id", "area_use", "geometry"],   # only these column chunks read
)

3. Vectorized geometry math. Buffers, centroids, simplification, and area all apply to the entire column at once. Compute metric quantities only in a projected CRS (see the next section):

parcels_utm = parcels.to_crs("EPSG:25832")              # metric CRS for area/buffer
parcels_utm["area_m2"] = parcels_utm.area               # square metres, vectorized
parcels_utm["setback"] = parcels_utm.buffer(-3.0)       # 3 m inward setback

4. Spatial predicates and joins. Enrich features by their location using an R-tree-backed spatial join — the workhorse of attribute integration, covered in depth in Spatial Joins & Merging:

zones = gpd.read_file("zoning_districts.gpkg").to_crs(parcels_utm.crs)
enriched = gpd.sjoin(parcels_utm, zones[["zone_code", "geometry"]],
                     how="left", predicate="within")

5. Dissolve and aggregate. Collapse many features into summary geometry grouped by an attribute — the spatial equivalent of a GROUP BY:

by_use = parcels_utm.dissolve(by="area_use", aggfunc={"area_m2": "sum"})

6. Clip and overlay. Where a join answers "which zone is this parcel in", an overlay answers "what is the geometry of the part that falls inside" — it cuts new features at the boundary rather than copying attributes across it. Use clip for the common case of trimming a layer to a study area, and overlay when you need the intersection, difference or union as geometry; the algebra and its performance profile are covered in Geometric Intersections & Overlays:

study_area = gpd.read_file("floodplain_boundary.gpkg").to_crs(parcels_utm.crs)
in_floodplain = parcels_utm.clip(study_area)                    # trimmed geometry
exposed = gpd.overlay(parcels_utm, study_area, how="intersection")
exposed["exposed_m2"] = exposed.area                            # measurable share

7. Assign a stable feature identity. Every ingestion pipeline that runs more than once needs to answer "is this the same parcel as last month". Source keys are unreliable — vendors renumber, merge and re-issue — so derive a deterministic identifier from the fields that genuinely identify the feature, and never from the geometry alone, which changes with every resurvey and every reprojection round trip:

import hashlib
import pandas as pd

def feature_key(row) -> str:
    payload = f"{row['municipality_code']}|{row['parcel_id']}"
    return hashlib.blake2b(payload.encode("utf-8"), digest_size=8).hexdigest()

parcels_utm["feature_key"] = parcels_utm.apply(feature_key, axis=1)
assert not parcels_utm["feature_key"].duplicated().any(), "key is not unique"

Chaining these operations — read, prune, transform, join, dissolve, overlay, key — expresses most ingestion jobs without a single explicit loop over rows. The one operation that resists vectorization is anything reaching into a geometry's internals; there, work against gdf.geometry.values with a Shapely 2 ufunc rather than .apply(), which pays a Python dispatch per row. When a dataset outgrows one machine, the same API scales out through Dask-GeoPandas, discussed under Production Patterns below.

CRS / Projection Considerations

Coordinate reference systems are the single largest source of silently wrong results in spatial pipelines, so this stage deserves disproportionate care; the full transformation reference lives in Coordinate Reference System Transformations and the underlying PROJ mechanics in Coordinate Systems with PyProj. Three rules prevent nearly all CRS defects.

Normalize once, immediately after ingestion. Assign a CRS if one is missing, then reproject every dataset to a single project-wide target before any analysis. Never let two frames with different CRSs meet in a join.

import geopandas as gpd

surveys = gpd.read_file("field_survey.shp")

# set_crs only labels; to_crs actually transforms coordinates — never confuse them
if surveys.crs is None:
    surveys = surveys.set_crs("EPSG:4326")   # label existing lon/lat data
surveys_metric = surveys.to_crs("EPSG:25832")  # reproject to ETRS89 / UTM 32N

Choose a projected CRS for anything metric. Areas, lengths, buffers, and distances are meaningless in geographic degrees. Use a locally appropriate projected system — a UTM zone, a national grid such as ETRS89/UTM or British National Grid (EPSG:27700) — and specifically avoid Web Mercator (EPSG:3857) for measurement: its scale distortion reaches tens of percent away from the equator, so a "500 m" buffer in EPSG:3857 is not 500 m on the ground.

Decision tree for choosing a target CRS A decision tree rooted in the question of what the coordinates are for. One branch, measuring on the ground — area, length, buffer, distance — leads to a projected CRS such as a local UTM zone or national grid, and away from Web Mercator, whose scale error grows with latitude. The other branch, storing or drawing, leads to EPSG:4326 for lon/lat interchange and to EPSG:3857 for tile rendering only. A closing band says to confirm the CRS unit name reads metre before measuring. Pick the CRS from what the numbers must mean What is the CRS for? Measuring on the ground area · length · buffer · distance Storing or drawing interchange · tiles · web display Projected CRS local UTM zone or national grid EPSG:25832 · 27700 Never Web Mercator EPSG:3857 scale error grows with latitude by tens of percent EPSG:4326 lon/lat storage and interchange EPSG:3857 tile rendering only, never math Confirm before measuring: the active CRS unit name must read metre, not degree.
The target CRS follows from the job: metric work needs a local projected system, interchange stays in lon/lat, and Web Mercator is a rendering grid only.

Mind axis order and deprecated syntax. PROJ 6+ honours the authority-defined axis order, which for EPSG:4326 is latitude-then-longitude. GeoPandas and Shapely store coordinates x/y (lon/lat), so when you hand-build a transformer, pass always_xy=True to keep everything lon/lat. And never use the removed +init=epsg: PROJ string form — it is rejected outright by PROJ 6+.

from pyproj import Transformer

# always_xy=True -> input & output are (lon, lat) / (x, y), never (lat, lon)
tf = Transformer.from_crs("EPSG:4326", "EPSG:25832", always_xy=True)
easting, northing = tf.transform(11.34, 44.49)   # (lon, lat) in -> (E, N) out
print(round(easting), round(northing))           # ~ 685000 4928000

# WRONG — raises CRSError on PROJ 6+:
# CRS("+init=epsg:4326")

Verify the transform survived by checking units after reprojection: gdf.crs.axis_info[0].unit_name should read metre, not degree, before you compute a single area.

Architecturally, the target CRS is a property of the project, not of any one script. Put it in configuration, import it everywhere, and never let a module hard-code a second code — the failure mode of a pipeline with two target CRSs in it is not an exception but a join that returns fewer rows than it should, months after the second code was introduced. Where the extent is not known in advance, derive the target from the data itself rather than guessing a zone, using estimate_utm_crs() as described in Choosing a UTM Zone Automatically in Python:

import geopandas as gpd

sensors = gpd.read_file("river_gauges.geojson")      # RFC 7946: lon/lat, EPSG:4326
target = sensors.estimate_utm_crs()                  # zone derived from the extent
print(target.to_authority())                         # ('EPSG', '32632') for 6-12°E
sensors_metric = sensors.to_crs(target)

That helper is only correct while the data fits inside one zone. Continental extents need an equal-area CRS instead, and a layer that spans several zones and is reprojected into one of them will carry growing distortion at its edges — the reason the choice belongs to the pipeline's configuration rather than to whichever function happens to need metres first.

Production Patterns & Performance

Pipelines that pass on a laptop sample fail on the full archive for predictable reasons: memory, I/O, and index reuse. Address all three deliberately.

Bound memory with windowed raster reads. Never read a full scene into RAM. Rasterio's windowing streams one tile at a time, which is mandatory for cloud-optimized GeoTIFFs and large DEMs:

import rasterio
from rasterio.windows import Window

with rasterio.open("elevation_dem.tif") as dem:
    for row in range(0, dem.height, 1024):
        for col in range(0, dem.width, 1024):
            win = Window(col, row, 1024, 1024)
            tile = dem.read(1, window=win, masked=True)   # one 1024² tile at a time
            # ... process tile, write result, release memory before the next window

Build the spatial index once, and rebuild it after edits. Accessing .sindex triggers R-tree construction and caches it. Any operation that mutates geometry (repair, reproject, explode) invalidates the cached tree — force a rebuild before the next spatial query, or joins silently use stale bounds.

import geopandas as gpd
from shapely.validation import make_valid

boundaries = gpd.read_file("admin_boundaries.gpkg")
invalid = ~boundaries.is_valid
print(f"invalid geometries: {invalid.sum()}")
if invalid.any():
    boundaries.loc[invalid, "geometry"] = boundaries.loc[invalid, "geometry"].apply(make_valid)
_ = boundaries.sindex        # rebuild R-tree after modifying geometry

Full validation strategy — self-intersection detection, ring orientation, sliver removal — is covered in Topology Validation & Repair.

Store in a columnar, cloud-native format. GeoParquet and FlatGeobuf give columnar compression, predicate pushdown, and partitioned reads — a decisive upgrade over shipping Shapefiles between stages. The trade-offs across GeoParquet, FlatGeobuf, and cloud-optimized rasters are detailed in Cloud-Native Geospatial Formats:

parcels_metric.to_parquet("stage/parcels.parquet")       # columnar, CRS preserved

Scale out when a single machine is not enough. Terabyte archives need partitioned execution. Dask-GeoPandas mirrors the GeoPandas API over lazy partitions, adding a spatial shuffle so joins stay local — the full scaling playbook is in Scaling with Dask-GeoPandas:

import dask_geopandas as dgpd

ddf = dgpd.read_parquet("census/*.parquet")
ddf = ddf.spatial_shuffle(by="geometry")     # co-locate nearby features per partition
result = ddf.to_crs("EPSG:25832").buffer(250).compute()
result.to_parquet("stage/census_buffers.parquet")

Push heavy queries into a database. When the workload is repeated aggregation or joins across large tables, hand it to an engine built for it: PostGIS for a transactional spatial database, or DuckDB spatial for embedded analytics that read GeoParquet directly. The Python pipeline becomes the loader and orchestrator rather than the compute engine.

Handle complex geometry explicitly. Real archives contain multipart features, Z/M dimensions, and non-planar networks that break naive cleaning. Strip dimensions and explode multiparts before planar analysis:

import geopandas as gpd
from shapely import force_2d

network = gpd.read_file("transport_network.gpkg")
network["geometry"] = force_2d(network.geometry.values)      # drop Z/M for 2D work
network = network.explode(index_parts=False).reset_index(drop=True)
assert network.geom_type.isin(["Point", "LineString", "Polygon"]).all()

Cross the vector–raster boundary deliberately. Most production pipelines eventually have to combine the two data models, and there are only two directions to do it in. Burning vectors into the grid (rasterio.features.rasterize) turns polygons into a mask or a categorical band, which is what you want when the raster is the analytical unit — a land-cover reclassification, a cost surface, a training mask. Sampling the grid under vectors goes the other way and produces one row of statistics per feature, which is what you want when the feature is the analytical unit. The rule that decides the CRS is asymmetric: reproject the vectors to the raster, never the raster to the vectors, because warping resamples every pixel and degrades the data while reprojecting geometry is lossless.

import geopandas as gpd
import rasterio
from rasterio.features import rasterize

with rasterio.open("elevation_dem.tif") as dem:
    floodplain = gpd.read_file("floodplain_boundary.gpkg").to_crs(dem.crs)
    mask = rasterize(
        ((geom, 1) for geom in floodplain.geometry),
        out_shape=(dem.height, dem.width),
        transform=dem.transform,
        fill=0,
        all_touched=False,       # True inflates area at the polygon edge
        dtype="uint8",
    )

all_touched is the argument to think about rather than accept: the default includes only pixels whose centre falls inside the polygon, while True includes any pixel the boundary clips, which systematically overstates area for small features. The per-feature direction, including partial-pixel weighting, is covered in Zonal Statistics & Raster Sampling.

Make every stage idempotent. A spatial pipeline reruns constantly — a supplier reissues a delivery, a node is preempted, someone fixes a bug and replays a month. Any stage that appends rather than replaces will double its output on the second run, and spatial data is unusually good at hiding that: duplicated polygons overlap perfectly, so the map looks right while every area sum is twice what it should be. Write each stage's output to a temporary path and move it into place atomically on success, key partitions by their inputs so a rerun overwrites the same partition, and check for exact duplicates on the stable feature key before the data leaves the stage.

import shutil
import tempfile
from pathlib import Path
import geopandas as gpd

def write_partition(gdf: gpd.GeoDataFrame, destination: str) -> None:
    """Replace-on-success: a failed run leaves the previous artefact intact."""
    assert not gdf["feature_key"].duplicated().any(), "duplicate keys in partition"
    destination = Path(destination)
    destination.parent.mkdir(parents=True, exist_ok=True)
    with tempfile.TemporaryDirectory(dir=destination.parent) as tmp:
        staged = Path(tmp) / destination.name
        gdf.to_parquet(staged)
        shutil.move(staged, destination)     # atomic within one filesystem

Partition on the column you filter by. A folder of GeoParquet parts is only fast if a reader can skip most of it. Hive-style partitioning by an administrative code — region=BY/parcels.parquet — lets both Dask and DuckDB prune whole directories from the query plan before opening a file, which beats any spatial index for the "one region at a time" access pattern that dominates municipal and national workloads. Within a partition, sorting by a spatial key keeps each row group's bounding box tight so the per-file statistics are worth consulting.

Measure the data, not just the run. A pipeline that succeeds is not the same as a pipeline that is right. Emit a small set of numbers per stage and store them next to the artefact: feature count, the share of invalid geometries before repair, the count of null and empty geometries pruned, the total bounds, and the set of column names. The value is in the deltas — a validity rate that drops from 99.8% to 91% between two deliveries is a supplier regression, and a total_bounds that shifts by a few hundred kilometres is a CRS mislabelling. Neither raises an exception, and neither is visible without a record to compare against.

import geopandas as gpd

def stage_metrics(gdf: gpd.GeoDataFrame, stage: str) -> dict:
    minx, miny, maxx, maxy = gdf.total_bounds
    return {
        "stage": stage,
        "features": len(gdf),
        "invalid": int((~gdf.geometry.is_valid).sum()),
        "empty": int(gdf.geometry.is_empty.sum()),
        "missing": int(gdf.geometry.isna().sum()),
        "epsg": gdf.crs.to_epsg() if gdf.crs else None,
        "bounds": [round(v, 1) for v in (minx, miny, maxx, maxy)],
        "columns": sorted(c for c in gdf.columns if c != "geometry"),
    }

Schema drift deserves its own alarm rather than a downstream KeyError. Comparing the recorded columns list against the previous run catches a supplier who renamed zone_code to zoning_code, which is otherwise discovered by a join that suddenly matches nothing.

Know where each tool stops. The ladder is fairly predictable on ordinary hardware. Plain GeoPandas is comfortable while the live geometry fits in perhaps a third of RAM — a few million simple features on a 32 GB machine, far fewer for vertex-dense polygons. Beyond that, streaming batches keep memory flat at the cost of losing whole-dataset operations, and Dask-GeoPandas buys out-of-core execution at the cost of a shuffle whenever features must meet their neighbours. PostGIS wins when many clients query the same indexed data concurrently, and DuckDB wins for single-process analytical scans over columnar files. Moving up the ladder before you have to costs more in complexity than it returns in throughput.

Common Mistakes

Frequently Asked Questions

How do I handle a pipeline with datasets in several different CRSs? Normalize to one project-wide target CRS immediately after ingestion. Assign a CRS with set_crs only if it is genuinely missing, then to_crs every frame to the same projected system before any join or measurement. Two frames with different CRSs must never meet in sjoin — the result is silent misalignment, not an error.

What is the most memory-efficient way to process large vector and raster files? For vectors, read in bounded chunks with pyogrio and store intermediates as partitioned GeoParquet; for rasters, use Rasterio windowed reads over cloud-optimized GeoTIFFs so only one tile is resident at a time. Both keep peak memory flat regardless of total file size.

When should I reach for Dask-GeoPandas, PostGIS, or DuckDB instead of plain GeoPandas? Stay on GeoPandas while data fits comfortably in RAM. Move to Dask-GeoPandas for out-of-core, multi-core execution of the same API; to PostGIS when you need a transactional, indexed spatial database shared across services; and to DuckDB spatial for fast embedded analytics that query GeoParquet directly without a server.

Should I keep data as Shapefiles between processing stages? No. Shapefiles truncate field names to 10 characters, cap at 2 GB, and split across sidecar files. Use GeoParquet or FlatGeobuf for intermediate storage — columnar compression, preserved CRS, and predicate pushdown make them faster and safer.

Why do my area and distance calculations look wrong? Almost always because the geometry is in a geographic CRS (degrees) or in Web Mercator. Reproject to a locally appropriate projected CRS — a UTM zone or national grid — and confirm gdf.crs.axis_info[0].unit_name reads metre before computing anything metric.

How do I make a spatial pipeline safe to rerun? Make every stage replace rather than append. Write each output to a temporary path and move it into place only on success, key partitions by their inputs so a replay overwrites the same files, and assert uniqueness on a stable feature key before the data leaves the stage. Duplication is the failure mode to design against, because overlapping duplicate polygons look correct on a map while doubling every count and area downstream.

What should I monitor to catch bad data before it reaches analysis? Four numbers per stage, recorded next to the artefact: feature count, the share of invalid geometries before repair, the count of null and empty geometries pruned, and the layer's total bounds — plus the sorted column list to catch schema drift. The absolute values matter less than the deltas between runs; a validity rate that falls sharply or bounds that jump by hundreds of kilometres are supplier and CRS problems respectively, and neither raises an exception on its own.

Where should the target CRS be defined in a multi-stage pipeline? In one configuration value that every stage imports, never hard-coded in the module that first needs metres. A pipeline containing two different target codes does not fail loudly — it produces joins that quietly match fewer features than they should. Where the extent is not known ahead of time, derive the code from the data with estimate_utm_crs() and record the answer in the run's metadata.

When do I use Shapely versus GeoPandas? Use Shapely for low-level, per-geometry logic and custom validation; use GeoPandas for tabular, vectorized operations across an entire column. In practice you write GeoPandas for the pipeline and drop to Shapely for the tricky single-geometry cases.

Related