Converting shapefiles to GeoParquet without losing CRS metadata

You convert a county parcel shapefile to GeoParquet with a single gdf.to_parquet() call, the file writes without error, and everything looks fine — until a downstream spatial join returns zero matches or, worse, silently misaligns parcels against a zoning overlay by a few hundred meters. You inspect the output and find the culprit: pyarrow.parquet.read_schema() shows the geo metadata crs field is null, or gpd.read_parquet() loads a GeoDataFrame whose .crs is None. The geometry survived; the projection identity did not. Every coordinate is now an anonymous pair of numbers, and any consumer that assumes WGS84 will place a State Plane parcel in the wrong hemisphere of the coordinate space. This is a focused companion to the broader geospatial format conversion topic within automated feed ingestion and GIS data parsing, and it fixes one narrow failure: a shapefile-to-GeoParquet conversion that drops the CRS because the source .prj was missing, blank, or unrecognisable to PROJ.

Diagnosis: why the CRS comes out null jump to heading

GeoParquet stores its CRS inside a JSON blob in the file’s Arrow key-value metadata under the geo key, as PROJJSON. geopandas populates that blob from the GeoDataFrame’s .crs attribute at write time. If .crs is None when you write, the geo.columns.<geom>.crs field is emitted as null — no error, because a null CRS is a valid GeoParquet document. The loss almost always happens at read time, before you ever call to_parquet, for one of three reasons:

1. A missing or empty .prj sidecar. A shapefile’s CRS lives entirely in the companion .prj file. If the county packaged the .shp, .shx, and .dbf but omitted the .prj — extremely common with FTP dumps and zip exports — gpd.read_file() returns a GeoDataFrame with .crs = None and no warning.
2. A .prj PROJ cannot resolve. The file exists but contains an ESRI WKT variant, a truncated string, or a vendor projection that your PROJ build does not recognise, so OGR reads it as an unknown CRS and drops it rather than guessing.
3. A CRS lost in an earlier hop. The shapefile was itself produced by an upstream conversion that already discarded the projection, so the .prj misdeclares or omits the true CRS.

Reproduce it deterministically before changing anything. Inspect the sidecar and the loaded CRS directly:

from pathlib import Path
import geopandas as gpd

shp = Path("travis_county_parcels.shp")
prj = shp.with_suffix(".prj")

print("prj exists:", prj.exists())
print("prj bytes :", prj.stat().st_size if prj.exists() else 0)
print("prj text  :", prj.read_text().strip()[:120] if prj.exists() else "<missing>")

gdf = gpd.read_file(shp, engine="pyogrio")
print("loaded crs:", gdf.crs)          # None means the CRS is already gone
print("x range   :", gdf.total_bounds) # coordinate magnitude hints at the real CRS

A loaded crs: None confirms the diagnosis. The total_bounds are a strong clue to what the CRS should be: values like [3050000, 10050000, 3150000, 10120000] are feet in a State Plane zone, [-98.0, 30.1, -97.4, 30.6] are WGS84 degrees, and [600000, 3350000, 640000, 3380000] are UTM meters. Never write the file until the CRS is restored — a null CRS in GeoParquet is a silent, downstream-only failure.

Step-by-step implementation jump to heading

Each step isolates one concern; the final orchestration ties them together into a single guarded conversion.

Step 1 — Detect and inspect the source CRS jump to heading

Load the shapefile and branch on whether a CRS came through. Read the raw .prj bytes yourself rather than trusting that gpd.read_file surfaced the problem, because an empty or malformed .prj produces the same None as a missing one.

from pathlib import Path
from typing import Optional
import geopandas as gpd
from pyproj import CRS
from pyproj.exceptions import CRSError

def load_with_crs_status(shp_path: str) -> tuple[gpd.GeoDataFrame, Optional[CRS]]:
    """Return the layer and its resolved CRS, or (layer, None) if the CRS is gone."""
    gdf = gpd.read_file(shp_path, engine="pyogrio")
    if gdf.crs is not None:
        return gdf, CRS.from_user_input(gdf.crs)

    prj = Path(shp_path).with_suffix(".prj")
    if prj.exists() and prj.stat().st_size > 0:
        try:
            # The .prj exists but geopandas dropped it — try to parse it directly.
            return gdf, CRS.from_wkt(prj.read_text().strip())
        except CRSError:
            pass  # ESRI WKT variant PROJ can't read; fall through to inference.
    return gdf, None

Step 2 — Repair the CRS by inference when it is missing jump to heading

When no authoritative CRS survived, infer it deterministically. Use pyproj’s CRS.from_authority against a candidate list drawn from the county’s known projection, and validate the guess against the coordinate bounds so a wrong EPSG fails loudly instead of writing plausible-but-wrong data. Never guess blindly — pin candidates per jurisdiction.

from pyproj import CRS

# Candidate CRSs a given county is known to publish, most likely first.
COUNTY_CRS_CANDIDATES = {
    "travis_tx": [2277, 6578, 4326],   # TX Central State Plane (ft), then meters, then WGS84
    "king_wa":   [2926, 6597, 4326],   # WA North State Plane (ft)
}

def infer_epsg_from_bounds(gdf, candidates: list[int]) -> int:
    """Pick the candidate EPSG whose axis units match the coordinate magnitudes."""
    minx, miny, maxx, maxy = gdf.total_bounds
    looks_geographic = -180 <= minx <= 180 and -90 <= miny <= 90 and abs(maxx) <= 180
    for epsg in candidates:
        crs = CRS.from_epsg(epsg)
        is_geographic = crs.is_geographic
        # Match a degrees-looking extent to a geographic CRS, a large extent to projected.
        if looks_geographic == is_geographic:
            return epsg
    raise ValueError(
        f"No candidate in {candidates} matches bounds {gdf.total_bounds}; "
        "coordinate magnitudes fit neither geographic nor the projected candidates."
    )

Step 3 — Set the authoritative EPSG on the table jump to heading

Attach the resolved CRS with set_crs, never to_crs. set_crs labels the existing coordinates with their true identity; to_crs would transform them, which is exactly wrong when the coordinates are already correct and only the label was missing. This distinction is the single most common way a CRS repair goes wrong, and it is treated in depth in the CRS alignment strategies that govern projection identity across jurisdictions.

def apply_authoritative_crs(gdf, epsg: int):
    """Label coordinates with their true CRS (set), never reproject (to_crs)."""
    if gdf.crs is not None and gdf.crs.to_epsg() == epsg:
        return gdf  # already correct; do nothing
    # allow_override=True is required when a wrong/blank CRS is present.
    return gdf.set_crs(epsg=epsg, allow_override=True)

Step 4 — Write GeoParquet with embedded CRS and geo metadata jump to heading

With an authoritative CRS on the table, to_parquet embeds it as PROJJSON in the geo metadata automatically. Pin the GeoParquet schema version so the metadata layout is explicit and consumers know exactly what they are reading.

import geopandas as gpd

def write_geoparquet(gdf: gpd.GeoDataFrame, out_path: str) -> None:
    """Write GeoParquet only when the CRS is authoritative and embedded."""
    if gdf.crs is None:
        raise ValueError("Refusing to write GeoParquet with a null CRS.")
    # schema_version pins the geo metadata layout; write_covering_bbox aids pushdown.
    gdf.to_parquet(out_path, index=False, schema_version="1.0.0")

Step 5 — Round-trip verify the embedded CRS jump to heading

Do not trust that the write succeeded — prove it. Read the file back and both the GeoDataFrame CRS and the raw Arrow geo metadata must agree with the EPSG you set. This closes the loop the diagram shows, and a mismatch routes the layer back to inference rather than into the warehouse.

import json
import geopandas as gpd
import pyarrow.parquet as pq

def verify_geoparquet_crs(out_path: str, expected_epsg: int) -> None:
    """Assert the round-tripped CRS matches; raise so a bad file never ships."""
    back = gpd.read_parquet(out_path)
    if back.crs is None or back.crs.to_epsg() != expected_epsg:
        raise AssertionError(f"Round-trip CRS {back.crs} != EPSG:{expected_epsg}")

    # Confirm the CRS is physically present in the geo metadata, not just inferred.
    meta = pq.read_schema(out_path).metadata or {}
    geo = json.loads(meta.get(b"geo", b"{}").decode() or "{}")
    primary = geo.get("primary_column")
    crs_blob = geo.get("columns", {}).get(primary, {}).get("crs")
    if crs_blob is None:
        raise AssertionError("geo metadata present but crs field is null")

Verification and testing jump to heading

Confirm the repair fixed the identity rather than masking it:

• Non-null CRS on read. gpd.read_parquet(out).crs must be a concrete CRS whose .to_epsg() equals the EPSG you set — not None, and not a different code.
• CRS lives in the file, not the reader. Inspect the raw geo metadata with pyarrow.parquet.read_schema; the crs field must be non-null PROJJSON, proving the projection is embedded rather than reconstructed from a filename.
• Coordinates are unchanged. Because you used set_crs, the numeric coordinates before and after must be identical to the bit — assert gdf.total_bounds == back.total_bounds. Any drift means to_crs sneaked in and reprojected data that was already correct.
• A join lands correctly. Overlay a handful of known parcels against a zoning layer of a known CRS; the intersections must match ground truth. This is the test that catches an inferred-but-wrong EPSG that all the metadata checks pass.
• Round-trip through the source CRS. Transform the output to WGS84 and back to the source EPSG; a point should return within centimetres. Large drift signals the inferred zone was wrong.

Failure recovery jump to heading

When a shapefile defeats inference, keep the batch moving and stay reproducible:

• Quarantine, never default. A layer whose CRS cannot be resolved to a confident EPSG is written to a crs_unresolved/ partition with its bounds, .prj bytes, and candidate list, and the run continues. Never fall back to writing it as WGS84 — that manufactures a false projection identity that is far harder to detect later than a null one.
• Escalate ambiguous inference. When two candidates both match the bounds heuristic (feet versus survey-feet State Plane zones differ by ~2 ppm and can be indistinguishable at a glance), flag the layer for human confirmation against a landmark coordinate rather than picking one silently.
• Record the inference decision. Persist which EPSG was assigned, whether it came from the .prj or from bounds inference, and the candidate list, in the conversion manifest, so a later reviewer can audit why a projection was chosen — the same lineage discipline applied across the geospatial format conversion topic.
• Re-request the sidecar. If a source repeatedly ships shapefiles without a .prj, the durable fix is upstream: request the projection file from the county or pin its CRS in a per-jurisdiction lookup, rather than inferring on every run.

Frequently asked questions jump to heading

Related jump to heading

• Parent topic: Geospatial Format Conversion
• Section overview: Automated Feed Ingestion & GIS Data Parsing
• CRS Alignment Strategies — setting authoritative EPSG and distinguishing labeling from reprojection

Up: Geospatial Format Conversion