Converting Shapefiles to GeoJSON for Epi Pipelines

This guide is part of Spatial Data Types & Formats, and it solves one narrow operational problem: turning a legacy ESRI Shapefile into RFC 7946 GeoJSON that a surveillance dashboard or web map can consume without silently corrupting geometry, drifting coordinates, or leaking protected health information (PHI).

Problem Context & Constraints

Legacy ESRI Shapefiles persist in public health data exchanges because agency GIS workflows, state health-department portals, and federal boundary releases were standardized on the format decades ago. But modern syndromic surveillance dashboards, real-time case feeds, and interoperable web mapping stacks require strict GeoJSON. The naive export path — ogr2ogr in.shp out.geojson or a one-line gpd.read_file(...).to_file(...) — fails in this exact scenario for four compounding reasons:

• The Shapefile carries no opinion about its target CRS. A blind export writes whatever the source .prj declares. RFC 7946 mandates WGS 84 (EPSG:4326) with longitude-latitude ordering; a UTM, State Plane, or custom health-district grid that survives export will render in the Gulf of Guinea, not the county it describes.
• Legacy geometries are routinely invalid. Self-intersecting polygons, unclosed rings, and reversed winding orders pass quietly through a default export and then break spatial joins, weights-matrix construction, and hotspot routines downstream.
• Float noise becomes a payload and a privacy problem. Source coordinates often carry 10–15 decimals; serialized verbatim they bloat the file and encode a household to sub-millimeter precision.
• The attribute table ships PHI as-is. Patient identifiers, exact addresses, and free-text clinical notes ride along in the .dbf unless they are explicitly dropped.

A defensible converter therefore applies each safeguard in a fixed, auditable order rather than trusting the driver’s defaults. The same ordering discipline governs the broader Spatial Epidemiology Fundamentals & Data Standards layer this page sits under.

The stage order is load-bearing: CRS alignment must precede precision truncation (see the Coordinate Reference Systems for Public Health constraints), and geometry repair must precede truncation so that rounding never re-introduces a self-intersection on an already-valid ring.

Prerequisites

Pin the GIS stack so the conversion is reproducible across CI runners and analyst laptops:

# requirements.txt (pinned)
geopandas==1.0.1
shapely==2.0.6     # provides make_valid + force_2d
pyproj==3.6.1      # backs the datum transform
numpy==2.1.1

Input requirements for the converter below:

• A complete Shapefile sidecar set — .shp, .shx, .dbf, and crucially .prj. A missing .prj means the source CRS is undeclared and the run must abort rather than guess.
• A known list of PHI/PII column names for the source schema (for example patient_id, street_address, dob, mrn). This list is the de-identification contract and is logged with the run.
• Target CRS fixed to EPSG:4326 for RFC 7946 conformance; precision fixed by policy (six decimals here, but configurable to align with Precision Standards in Epi-Mapping).

Step-by-Step Solution

The implementation enforces datum-verified CRS alignment, geometry repair, post-transform precision truncation, and attribute sanitization, and it emits a manifest for the audit trail. It is designed for headless execution and CI/CD integration.

# Python 3.11 · geopandas==1.0.1 · shapely==2.0.6 · pyproj==3.6.1
import logging
import json
import re
import hashlib
from datetime import datetime, timezone
from typing import Optional, List
import geopandas as gpd
from shapely import force_2d
from shapely.validation import make_valid
from shapely.ops import transform

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(message)s",
)


def _round_coords(geom, precision: int = 6):
    """Recursively round every coordinate in a Shapely geometry."""
    def _round_xy(x, y, z=None):
        rounded = (round(x, precision), round(y, precision))
        return rounded if z is None else (*rounded, round(z, precision))
    return transform(_round_xy, geom)


def _sanitize_columns(gdf: gpd.GeoDataFrame, pii_columns: List[str]) -> gpd.GeoDataFrame:
    """Drop PHI/PII columns, then standardize remaining field names to snake_case."""
    drop_cols = [c for c in pii_columns if c in gdf.columns]
    if drop_cols:
        logging.info("Dropping PHI/PII columns: %s", drop_cols)
        gdf = gdf.drop(columns=drop_cols)
    gdf.columns = [re.sub(r"[^a-z0-9_]+", "_", c.lower()).strip("_") for c in gdf.columns]
    return gdf


def convert_shapefile_to_geojson(
    shp_path: str,
    out_path: str,
    pii_columns: List[str],
    target_crs: str = "EPSG:4326",
    precision: int = 6,
    drop_null_geoms: bool = True,
) -> dict:
    """Deterministic Shapefile -> RFC 7946 GeoJSON converter for epi pipelines.

    Order is fixed and audited: CRS align -> repair -> 2D -> truncate -> sanitize.
    Returns an audit manifest describing the transformations applied.
    """
    manifest = {"source": shp_path, "target_crs": target_crs, "precision": precision}

    logging.info("Loading shapefile: %s", shp_path)
    gdf = gpd.read_file(shp_path, encoding="utf-8")

    # Deterministic ordering: sort by a stable key so reruns are byte-identical.
    sort_key = "FID" if "FID" in gdf.columns else gdf.columns[0]
    gdf = gdf.sort_values(by=sort_key, kind="mergesort").reset_index(drop=True)

    # 1. CRS verification & datum-aware transformation.
    if gdf.crs is None:
        raise ValueError("Source .prj is missing — CRS undeclared. Abort, do not guess.")
    manifest["source_crs"] = gdf.crs.to_string()
    logging.info("Source CRS: %s | Target CRS: %s", manifest["source_crs"], target_crs)
    if str(gdf.crs) != target_crs:
        gdf = gdf.to_crs(target_crs)  # pyproj selects the datum transform pipeline
        logging.info("CRS transformation applied via pyproj.")
    minx, miny, maxx, maxy = gdf.total_bounds
    if minx < -180 or maxx > 180 or miny < -90 or maxy > 90:
        raise ValueError(f"Post-transform bounds out of WGS84 range: {gdf.total_bounds}")

    # 2. Geometry validation & repair (BEFORE truncation).
    invalid_mask = ~gdf.geometry.is_valid
    manifest["repaired_geometries"] = int(invalid_mask.sum())
    if invalid_mask.any():
        logging.warning("Repairing %d invalid geometries.", invalid_mask.sum())
        gdf.loc[invalid_mask, "geometry"] = gdf.loc[invalid_mask, "geometry"].apply(make_valid)
    if drop_null_geoms:
        null_mask = gdf.geometry.is_empty | gdf.geometry.isna()
        manifest["dropped_null_geometries"] = int(null_mask.sum())
        if null_mask.any():
            logging.info("Dropping %d null/empty geometries.", null_mask.sum())
            gdf = gdf[~null_mask].copy()

    # 3. Force 2D, then truncate precision (post-transform, post-repair).
    gdf["geometry"] = gdf.geometry.apply(force_2d)
    logging.info("Truncating coordinates to %d decimal places.", precision)
    gdf["geometry"] = gdf.geometry.apply(lambda g: _round_coords(g, precision))

    # 4. Attribute sanitization (the de-identification contract).
    manifest["dropped_pii_columns"] = [c for c in pii_columns if c in gdf.columns]
    gdf = _sanitize_columns(gdf, pii_columns)

    # 5. Serialization — minified RFC 7946, bbox stripped.
    geojson_dict = json.loads(gdf.to_json())
    geojson_dict.pop("bbox", None)
    payload = json.dumps(geojson_dict, separators=(",", ":"), ensure_ascii=False)
    with open(out_path, "w", encoding="utf-8") as f:
        f.write(payload)

    manifest["feature_count"] = len(geojson_dict["features"])
    manifest["output_sha256"] = hashlib.sha256(payload.encode("utf-8")).hexdigest()
    manifest["completed_utc"] = datetime.now(timezone.utc).isoformat()
    with open(out_path + ".manifest.json", "w", encoding="utf-8") as f:
        json.dump(manifest, f, indent=2)

    logging.info("Conversion complete. Features: %d | sha256: %s",
                 manifest["feature_count"], manifest["output_sha256"][:12])
    return manifest

The force_2d call strips Z/M dimensions that legacy survey shapefiles carry; the encoding="utf-8" read prevents mojibake in non-English jurisdictional names; and the SHA-256 of the serialized payload is the reproducibility anchor an auditor replays against.

Validation & Edge Cases

Three failure modes recur in real agency data. Each has a deterministic diagnostic.

1. Undeclared source CRS. A .prj-less Shapefile reads with gdf.crs is None, and the converter aborts rather than emitting silently wrong coordinates:

ValueError: Source .prj is missing — CRS undeclared. Abort, do not guess.

2. Coordinates outside WGS84 bounds after transform. A mislabeled .prj (for example a State Plane file tagged as EPSG:4326) survives to_crs() and then trips the bounds gate, surfacing the error before serialization rather than on a broken map:

ValueError: Post-transform bounds out of WGS84 range: [ 1320544.2  442118.9  1402887.5  501233.0 ]

3. Duplicate features after repair. make_valid() can split a self-touching polygon into a MultiPolygon, and overlapping administrative boundaries can produce coincident features. When aggregation is required, dissolve on the stable identifier after the repair stage:

# Collapse repaired multi-parts back to one feature per administrative unit.
gdf = gdf.explode(index_parts=False).dissolve(by="geoid").reset_index()
assert gdf["geoid"].is_unique, "Duplicate GEOIDs survived dissolve — inspect input overlaps."

A fourth quieter case: date fields. Shapefiles store dates as strings or integers; convert them to ISO 8601 before serialization so temporal joins in downstream syndromic feeds stay consistent.

Compliance Notes

For regulatory defensibility, the manifest emitted alongside the GeoJSON must persist, at minimum: source_crs and target_crs (proves the declared datum transform), precision (proves the coordinate resolution policy that governs de-identification headroom), dropped_pii_columns (the executed de-identification contract), repaired_geometries and dropped_null_geometries (proves no silent data loss), and output_sha256 with completed_utc (proves byte-level reproducibility). Sorting by a stable key before processing guarantees that an auditor replaying the run obtains an identical hash.

Note the boundary of what truncation buys you: six decimals (~0.11 m) is a payload and float-noise control, not a privacy threshold — it still encodes a residence. Coordinate-level suppression and aggregation are governed separately by the Compliance Mapping Frameworks gates, which decide whether a geometry is releasable at all before this converter ever runs.

Related Topics

• Spatial Data Types & Formats — the parent guide on format selection, validation, and conversion standards.
• Coordinate Reference Systems for Public Health — datum and projection enforcement that the CRS-alignment stage depends on.
• Precision Standards in Epi-Mapping — deterministic coordinate rounding and tolerance policy behind the truncation stage.