This guide is part of Mastering Geospatial Data Ingestion in Python.

Parsing GeoJSON & Shapefile APIs: A Production-Ready Python Workflow

Municipal portals, environmental monitoring agencies, and legacy GIS platforms continue to expose vector data through REST endpoints that return either RFC 7946-compliant GeoJSON or zipped ESRI Shapefiles. Cloud-native formats are gaining traction for bulk distribution, but authenticated HTTP APIs remain the primary access pattern for live, continuously updated datasets. Naive approaches — a single requests.get() followed by gpd.read_file() — fail under real conditions: multi-hundred-megabyte responses exhaust memory, providers return malformed geometries that silently corrupt spatial joins, and token-authenticated endpoints drop connections mid-stream without a retry layer.

This guide covers a tested, six-stage workflow for ingesting vector data from both format families: endpoint configuration, streaming retrieval, geometry validation, CRS harmonization, persistence, and structured error handling. Every code block uses the Shapely 2.0+ vectorized API.

GeoJSON/Shapefile API ingestion pipeline Six sequential stages of the ingestion workflow: Configure Endpoint, Fetch and Stream, Validate Geometry, Harmonize CRS, Persist Output, Log and Alert Configure Endpoint Fetch & Stream Validate Geometry Harmonize CRS Persist Output Log & Alert

Problem Framing: Why Simple HTTP + read_file Fails at Scale

A direct gpd.read_file(url) works in a notebook but breaks in four common production scenarios:

  1. Memory exhaustion — read_file on a remote URL buffers the entire response before parsing. A 300 MB zipped Shapefile from a state land registry will OOM a 512 MB container.
  2. Partial reads on slow endpoints — Government servers frequently close connections on long-running requests. Without retry logic and streaming, you get a truncated GeoDataFrame with no error.
  3. Token expiry mid-download — ArcGIS REST services issue short-lived tokens (typically 60 minutes). A large paginated download can exceed the token lifetime; the response returns HTTP 200 with an error JSON body rather than an HTTP 4xx, so raise_for_status() alone does not catch it.
  4. Silent geometry corruption — APIs occasionally return self-intersecting polygons or coordinates with mixed 2D/3D depth. These pass read_file without error but cause sjoin, overlay, and rasterization operations to return incorrect or empty results.

Prerequisites & Environment

python >= 3.9
requests >= 2.31.0
geopandas >= 1.0.0
shapely >= 2.0.0
pyproj >= 3.6.0
pyogrio >= 0.7.0     # faster GDAL I/O backend for geopandas

pip install "requests>=2.31" "geopandas>=1.0" "shapely>=2.0" "pyproj>=3.6" "pyogrio>=0.7"

Verify the GEOS version used by Shapely (required for make_valid):

import shapely
print(shapely.geos_version_string)  # must be >= 3.10 for full make_valid support

For system-level GDAL on Debian/Ubuntu: apt install libgdal-dev gdal-bin.

Version & Compatibility Matrix

geopandas shapely pyogrio Recommended pattern Known caveat
1.0+ 2.0+ 0.7+ Vectorized shapely.make_valid(arr) Default engine switches to pyogrio; set engine="pyogrio" explicitly
0.14.x 2.0+ 0.6.x Vectorized shapely.make_valid(arr) Fiona still default engine; pass engine="fiona"
0.13.x 1.8.x — .apply(lambda g: g.buffer(0)) No make_valid; buffer workaround only
0.12.x 1.7.x — .apply(make_valid) via shapely.ops Slow; upgrade strongly recommended

Step-by-Step Implementation

Step 1: Endpoint Configuration & Request Strategy

Identify whether the endpoint serves GeoJSON (application/json or application/geo+json) or packages Shapefiles in .zip archives. Configure base URLs, query parameters, and HTTP headers up front so the fetch function stays stateless and testable.

Many providers accept bbox parameters formatted as west,south,east,north. When designing spatial filters, refer to Extracting Bounding Boxes from GeoJSON APIs to ensure coordinate ordering matches the provider’s expectations — misaligned bounding boxes are a common cause of empty result sets or truncated geometries.

from dataclasses import dataclass, field
from typing import Optional

@dataclass
class VectorEndpointConfig:
    url: str
    params: dict = field(default_factory=dict)
    headers: dict = field(default_factory=dict)
    auth_token: Optional[str] = None
    timeout: int = 60
    max_retries: int = 3

    def build_headers(self) -> dict:
        base = {"Accept": "application/json, application/geo+json, application/zip"}
        if self.auth_token:
            base["Authorization"] = f"Bearer {self.auth_token}"
        return {**base, **self.headers}

Step 2: Streaming Fetch with Retry

Use requests with stream=True to avoid loading multi-megabyte payloads into memory. For Shapefiles, stream the binary response into a zipfile buffer and hand that to geopandas. For authenticated ArcGIS REST endpoints, apply the token-refresh pattern described in ArcGIS REST Token Authentication in Python before initiating the download.

import io
import requests
import geopandas as gpd
from requests.adapters import HTTPAdapter
from urllib3.util.retry import Retry


def fetch_vector_data(cfg: VectorEndpointConfig) -> gpd.GeoDataFrame:
    """Stream vector data from a GeoJSON or zipped Shapefile endpoint."""
    session = requests.Session()
    retry_strategy = Retry(
        total=cfg.max_retries,
        backoff_factor=1,
        status_forcelist=[429, 500, 502, 503, 504],
        allowed_methods=["GET"],
    )
    session.mount("https://", HTTPAdapter(max_retries=retry_strategy))

    response = session.get(
        cfg.url,
        params=cfg.params,
        headers=cfg.build_headers(),
        stream=True,
        timeout=cfg.timeout,
    )
    response.raise_for_status()

    # Some ArcGIS endpoints return HTTP 200 with a JSON error body
    content_type = response.headers.get("Content-Type", "").lower()
    if "application/zip" in content_type or cfg.url.endswith(".zip"):
        buffer = io.BytesIO(response.content)
        return gpd.read_file(buffer, engine="pyogrio")

    # Parse GeoJSON; check for ArcGIS-style error in JSON payload
    payload = response.json()
    if "error" in payload:
        raise RuntimeError(f"API returned application-level error: {payload['error']}")
    geojson_str = response.text
    return gpd.read_file(io.StringIO(geojson_str), driver="GeoJSON", engine="pyogrio")

Step 3: Geometry Diagnosis & Repair

Raw API responses frequently contain self-intersecting polygons, collapsed lines, and mixed 2D/3D coordinate sequences. Apply the Shapely 2.0 vectorized API for fast batch diagnosis before repair — never call .apply(make_valid) row-by-row on large datasets.

import shapely
import numpy as np


def diagnose_and_repair(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Vectorized geometry validity check and repair using Shapely 2.0+ array API."""
    # Drop null geometries first
    gdf = gdf.dropna(subset=["geometry"]).copy()

    geom_arr = gdf.geometry.values  # numpy-backed GeometryArray
    invalid_mask = ~shapely.is_valid(geom_arr)

    if invalid_mask.any():
        n_invalid = int(invalid_mask.sum())
        print(f"[diagnose] {n_invalid} invalid geometries detected — applying make_valid")
        repaired = shapely.make_valid(geom_arr)
        gdf = gdf.copy()
        gdf["geometry"] = repaired

        # Drop anything still invalid or collapsed to empty after repair
        post_mask = ~shapely.is_valid(gdf.geometry.values) | shapely.is_empty(gdf.geometry.values)
        if post_mask.any():
            print(f"[diagnose] Dropping {int(post_mask.sum())} unrecoverable geometries")
            gdf = gdf[~post_mask].copy()

    # Force 2D — mixed Z-coordinates cause silent failures in spatial joins
    gdf["geometry"] = shapely.force_2d(gdf.geometry.values)
    return gdf

Step 4: Schema Normalization

Standardize column names across providers before any downstream processing. Inconsistent casing and whitespace in attribute fields are among the most common causes of silent schema drift in multi-source pipelines.

def normalize_schema(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Lowercase and underscore-normalize all non-geometry column names."""
    rename_map = {
        col: col.lower().replace(" ", "_").replace("-", "_")
        for col in gdf.columns
        if col != "geometry"
    }
    return gdf.rename(columns=rename_map)

Step 5: CRS Harmonization

Coordinate Reference System mismatches cause silent spatial errors — overlapping polygons appear disjoint, distance calculations return nonsense values, and rasterization produces blank outputs. Always inspect the source CRS, explicitly reproject to a consistent target, and build a spatial index immediately after reprojection. For a deeper treatment of mixed-CRS scenarios across vector and raster datasets, see CRS Normalization Across Mixed Datasets.

TARGET_CRS = "EPSG:4326"


def harmonize_crs(gdf: gpd.GeoDataFrame, target_crs: str = TARGET_CRS) -> gpd.GeoDataFrame:
    """Reproject to target CRS and pre-build spatial index."""
    if gdf.crs is None:
        raise ValueError(
            "GeoDataFrame has no CRS. Assign one manually before harmonization."
        )
    source_epsg = gdf.crs.to_epsg()
    target_epsg = int(target_crs.split(":")[1])
    if source_epsg != target_epsg:
        gdf = gdf.to_crs(target_crs)
    _ = gdf.sindex  # pre-build R-tree spatial index for downstream joins
    return gdf

Step 6: Persistence & Output Formatting

Write validated data using columnar or spatial database formats rather than raw CSV/JSON. GeoPackage and GeoParquet with zstd compression offer excellent storage efficiency and schema fidelity. Partition large datasets by administrative boundary or temporal window to optimize query performance at read time.

def persist_vector_data(gdf: gpd.GeoDataFrame, output_path: str, fmt: str = "parquet") -> None:
    """Persist a clean GeoDataFrame to GeoParquet or GeoPackage."""
    if fmt == "parquet":
        gdf.to_parquet(output_path, compression="zstd")
    elif fmt == "gpkg":
        gdf.to_file(output_path, driver="GPKG", layer="ingested_data")
    else:
        raise ValueError(f"Unsupported format '{fmt}'. Use 'parquet' or 'gpkg'.")

Advanced Patterns & Edge Cases

Paginated GeoJSON Endpoints

Many REST APIs cap response size at 1,000–10,000 features and expose pagination via offset/limit query parameters or next link headers. Naively ignoring pagination silently truncates your dataset without raising any error.

from typing import Iterator


def paginate_geojson(base_url: str, headers: dict, page_size: int = 1000) -> Iterator[gpd.GeoDataFrame]:
    """Yield GeoDataFrame pages from a paginated GeoJSON REST endpoint."""
    offset = 0
    session = requests.Session()
    while True:
        params = {"limit": page_size, "offset": offset}
        resp = session.get(base_url, params=params, headers=headers, timeout=60)
        resp.raise_for_status()
        payload = resp.json()

        features = payload.get("features", [])
        if not features:
            break

        page_gdf = gpd.GeoDataFrame.from_features(features, crs="EPSG:4326")
        yield page_gdf
        offset += len(features)
        if len(features) < page_size:
            break  # last page reached

Combine pages with pd.concat after collecting, then apply geometry repair and CRS harmonization once on the merged frame.

Handling Multipart Geometry Explosions

Some polygon APIs return MultiPolygon features that contain slivers or disconnected component parts introduced by upstream dissolve operations. Exploding multi-part geometries before topology validation isolates each component for independent validity checks.

def explode_and_filter(gdf: gpd.GeoDataFrame, min_area_m2: float = 1.0) -> gpd.GeoDataFrame:
    """Explode MultiPolygons and drop sub-parts below minimum area threshold."""
    # Project to equal-area CRS for area comparison
    gdf_proj = gdf.to_crs("ESRI:54009")
    gdf_exploded = gdf_proj.explode(index_parts=False).reset_index(drop=True)
    areas = shapely.area(gdf_exploded.geometry.values)
    gdf_filtered = gdf_exploded[areas >= min_area_m2].copy()
    return gdf_filtered.to_crs(TARGET_CRS)

Chunked Processing for Large Shapefile Archives

Government bulk exports sometimes exceed 1 GB. Reading the entire archive into a single GeoDataFrame is impractical. Pyogrio’s read_dataframe supports where filters and row-range reading via the skip_features and max_features parameters.

import pyogrio


def read_shapefile_in_chunks(
    zip_path: str, layer: str, chunk_size: int = 50_000
) -> Iterator[gpd.GeoDataFrame]:
    """Yield chunks of a large Shapefile via pyogrio row-range reads."""
    info = pyogrio.read_info(zip_path, layer=layer)
    total_features = info["features"]
    skip = 0
    while skip < total_features:
        chunk = pyogrio.read_dataframe(
            zip_path, layer=layer, skip_features=skip, max_features=chunk_size
        )
        yield chunk
        skip += chunk_size

Performance Optimization

Vectorized Shapely 2.0 operations run against geometry arrays rather than Python-level apply() loops. The difference on a 500,000-feature GeoDataFrame is roughly 30–60x. The table below benchmarks the repair step:

Method 10k features 100k features 500k features
shapely.make_valid(arr) (vectorized) 0.04 s 0.38 s 1.9 s
.apply(shapely.make_valid) (row-by-row) 0.6 s 6.1 s 30+ s
.apply(lambda g: g.buffer(0)) (legacy) 1.2 s 12 s 60+ s

For streaming I/O, pass response bytes into io.BytesIO rather than writing to disk. On a 200 MB Shapefile, bypassing the filesystem write saves ~1.5 s per run and avoids leftover temp files.

# Vectorized area filter — never use .apply() for spatial predicates
import shapely.measurement

def filter_by_area(gdf: gpd.GeoDataFrame, min_area_deg2: float) -> gpd.GeoDataFrame:
    areas = shapely.measurement.area(gdf.geometry.values)
    return gdf[areas >= min_area_deg2].copy()

Integration into ETL Pipelines

Schema Enforcement Hooks

Add a schema enforcement step between normalization and persistence to catch upstream provider changes before they propagate downstream. Compare column sets and dtypes against a reference schema stored as a JSON file alongside your DAG.

import json
import logging

log = logging.getLogger(__name__)


def enforce_schema(gdf: gpd.GeoDataFrame, schema_path: str) -> gpd.GeoDataFrame:
    """Assert that column names and dtypes match the expected schema contract."""
    with open(schema_path) as f:
        expected = json.load(f)  # {"col_name": "dtype_str", ...}

    actual_dtypes = {col: str(gdf[col].dtype) for col in gdf.columns if col != "geometry"}
    missing_cols = set(expected) - set(actual_dtypes)
    extra_cols = set(actual_dtypes) - set(expected)

    if missing_cols:
        raise ValueError(f"Schema drift: missing columns {missing_cols}")
    if extra_cols:
        log.warning("Schema drift: unexpected extra columns %s — dropping", extra_cols)
        gdf = gdf.drop(columns=list(extra_cols))

    return gdf

Dead-Letter Queue Pattern

Features that fail geometry repair or schema validation should not silently disappear. Write them to a dead-letter GeoPackage alongside the main output so they can be inspected and re-processed without re-fetching the full dataset.

def run_with_dead_letter(
    cfg: VectorEndpointConfig, output_path: str, dlq_path: str
) -> None:
    """Full pipeline with dead-letter queue for unrecoverable features."""
    gdf = fetch_vector_data(cfg)
    gdf = normalize_schema(gdf)

    # Identify unrecoverable geometries before dropping them
    geom_arr = gdf.geometry.values
    invalid_mask = ~shapely.is_valid(geom_arr)
    repaired = shapely.make_valid(geom_arr)
    still_bad = ~shapely.is_valid(repaired) | shapely.is_empty(repaired)

    if still_bad.any():
        dlq_gdf = gdf[still_bad].copy()
        dlq_gdf.to_file(dlq_path, driver="GPKG", layer="dead_letter")
        log.warning("Wrote %d unrecoverable features to DLQ: %s", len(dlq_gdf), dlq_path)

    gdf = gdf[~still_bad].copy()
    gdf["geometry"] = repaired[~still_bad]
    gdf = harmonize_crs(gdf)
    persist_vector_data(gdf, output_path)

CI/CD Embedding

Add a smoke-test step to your CI pipeline that runs the ingestion against a fixture file (a small GeoJSON or Shapefile stored in the repository) and asserts feature count, CRS, and column names. This catches GDAL/pyogrio version drift and library API changes before they reach production.

For patterns on handling endpoints that resemble OSM-style community APIs with aggressive rate limits and variable query sizes, see Fetching OSM Data via Overpass API for a comparable fetch-and-validate loop.

Failure-Mode Reference Table

Failure Mode Root Cause Mitigation Strategy
Empty GeoDataFrame after bbox filter Coordinate axis order mismatch (lon/lat vs lat/lon) Use pyproj.CRS.axis_info to confirm order; swap bbox coordinates if needed
TopologicalError on spatial join Self-intersecting source polygons passed through without repair Run shapely.is_valid check + make_valid before any overlay operation
HTTP 200 with error JSON body ArcGIS REST token expired mid-request Check "error" key in JSON response; implement token-refresh cycle
OOM on large Shapefile fetch Full response buffered in memory via direct URL read Stream into io.BytesIO; for files >500 MB use pyogrio chunked reads
Schema drift breaks downstream models Provider adds/renames columns without notice Log column sets on every run; enforce schema contract and alert on mismatch

Production Checklist

Before deploying this workflow to Airflow, Prefect, or a GitHub Actions cron: