Spatial Epi Public Health GIS

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.

Shapefile to GeoJSON Conversion Pipeline Six stages connected left to right: Read legacy shapefile (.shp / .prj), Verify datum and align CRS to EPSG:4326, Validate and repair geometry (make_valid), Truncate coordinates to 6 decimals, Sanitize attributes (drop PHI, snake_case), and Serialize compliant GeoJSON. Read Shapefile .shp / .prj legacy input Align CRS verify datum, EPSG:4326 Repair Geometry make_valid, fix rings Truncate Coords 6 decimals ~0.11 m Sanitize Attrs drop PHI, snake_case Serialize compliant GeoJSON Conversion pipeline — a legacy shapefile becomes a privacy-safe, CRS-aligned GeoJSON

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.

Why the Conversion Order Cannot Be Swapped Two ordering rules shown as correct-versus-broken pairs. Rule one: aligning CRS to EPSG:4326 before truncating coordinates keeps precision expressed in degrees, whereas truncating projected meters first then reprojecting snaps vertices off their true position. Rule two: repairing geometry before truncating keeps rings valid, whereas truncating an unrepaired ring lets rounding re-introduce a self-intersection. Rule 1 — align CRS before you truncate Align CRS → EPSG:4326 Truncate 6 dp in degrees ✓ precision is in the target CRS ~0.11 m, vertices stay on position Truncate round meters Reproject → degrees ✗ rounding error rides the datum transform vertices snap off their true location Rule 2 — repair geometry before you truncate make_valid ring is valid Truncate stays valid ✓ no self-intersection re-introduced downstream spatial joins hold Truncating an unrepaired ring lets rounding collapse near-coincident vertices into a fresh self-intersection.

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.

The round trip is where the losses become visible, and they are not symmetric: converting back does not restore what the first hop discarded.

The Round Trip Is Not Lossless in Either Direction Three boxes in a row: the original shapefile, the GeoJSON produced from it, and a shapefile regenerated from that GeoJSON. The first hop discards the .prj coordinate reference in favour of CRS84, coerces date fields to strings, and normalises multipart geometry. The second hop truncates field names to ten characters, drops records with null geometry, and re-guesses attribute types from their values. Geometry coordinates survive both hops; almost nothing else does. What survives a shapefile round trip geometry coordinates survive both hops — the metadata that made them interpretable does not counties.shp .prj + .dbf types multipart, nulls counties.geojson CRS84 by definition JSON scalar types only counties_rt.shp not the file you started with, despite the name − .prj → CRS84 − dates → strings − multipart normalised − names cut to 10 chars − null geometry dropped − types re-guessed Convert once, in one direction, and keep the original — a round trip is a rewrite, not a copy

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.