Analyze spatial data with GeoPandas: production guide

Create and activate an isolated environment, then install GeoPandas and the packages needed for your workflow. A coordinated conda-forge environment can be useful when GDAL and other native libraries need compatible binaries. Pin the tested versions in the project configuration.

# pip example
python -m venv .venv

# macOS/Linux
source .venv/bin/activate

# Windows PowerShell
# .venv\\Scripts\\Activate.ps1

python -m pip install --upgrade pip
python -m pip install geopandas pyarrow pytest

# Alternative conda-forge example
# conda create -n geo-prod -c conda-forge python geopandas pyarrow pytest
# conda activate geo-prod

Read a small or filtered sample, then inspect CRS, bounds, geometry types, nulls, empties, and identifiers before applying any spatial operation.

import geopandas as gpd

points = gpd.read_file(
    "data/observations.gpkg",
    layer="observations",
)
regions = gpd.read_file(
    "data/regions.gpkg",
    layer="neighborhoods",
)

print("points CRS:", points.crs)
print("regions CRS:", regions.crs)
print(points.geometry.geom_type.value_counts(dropna=False))
print("missing:", points.geometry.isna().sum())
print("empty:", points.geometry.is_empty.sum())
print("bounds:", points.total_bounds)
print(points.head())

Use to_crs when transforming known coordinates. Select the spatial predicate deliberately and check the resulting cardinality because a spatial join can produce multiple rows for one input feature.

# Transform points to the polygon layer's CRS for the join.
points = points.to_crs(regions.crs)

joined = gpd.sjoin(
    points,
    regions[["region_id", "name", "geometry"]],
    how="left",
    predicate="within",
)

summary = (
    joined.groupby("region_id", dropna=False)
    .size()
    .rename("observation_count")
    .reset_index()
)

result = regions.merge(summary, on="region_id", how="left")
result["observation_count"] = (
    result["observation_count"]
    .fillna(0)
    .astype("int64")
)

Run geometry checks and business-level assertions before publishing. Include known interior, boundary, and unmatched examples in the test fixture, and investigate unexpected one-to-many matches rather than silently deduplicating them.

invalid = ~result.geometry.is_valid
empty = result.geometry.is_empty
missing = result.geometry.isna()

if invalid.any():
    raise ValueError(f"Invalid geometries: {invalid.sum()}")
if empty.any() or missing.any():
    raise ValueError("Output contains empty or missing geometry")

assert result["observation_count"].ge(0).all()
assert result["region_id"].is_unique

# Optional diagnostics for a source layer:
# from shapely.validation import explain_validity
# print(points.loc[invalid_points].geometry.map(explain_validity).head())

For larger sources, select only required columns and filter during ingestion where the driver supports it. Bbox and mask filters must use the dataset CRS. Test Arrow transfer separately because compatibility and driver behavior affect the result.

import pyogrio

observations = pyogrio.read_dataframe(
    "data/observations.gpkg",
    layer="observations",
    columns=["observation_id", "category", "reported_at"],
    bbox=(-87.8, 41.6, -87.5, 42.1),
    use_arrow=True,
)

# For a production pipeline, record the filter, source CRS,
# row count, elapsed time, and package versions in the run manifest.

Write the result with an explicit schema and CRS, then store a manifest containing source identifiers, filters, predicate, units, package versions, and timestamps. Use versioned output paths or idempotent writes so reruns do not mix source vintages.

from pathlib import Path
import json
from datetime import datetime, timezone

output_dir = Path("outputs")
output_dir.mkdir(exist_ok=True)

result.to_parquet(
    output_dir / "neighborhood_observation_counts.parquet",
    index=False,
)
result.to_file(
    output_dir / "neighborhood_observation_counts.gpkg",
    layer="counts",
    driver="GPKG",
)

manifest = {
    "created_at": datetime.now(timezone.utc).isoformat(),
    "crs": result.crs.to_string() if result.crs else None,
    "predicate": "within",
    "geometry_type": sorted(result.geometry.geom_type.dropna().unique().tolist()),
    "row_count": len(result),
}
(output_dir / "neighborhood_observation_counts.manifest.json").write_text(
    json.dumps(manifest, indent=2),
    encoding="utf-8",
)