Converting Parcel Shapefiles to GeoParquet for Batch Runs

Most parcel fabrics are still distributed as shapefiles, a format with a ten-character field-name limit, no proper date or boolean type, no null representation for numbers, and a read cost that dominates many batch pipelines. GeoParquet fixes all of that, but a careless conversion silently changes the data — truncated names become ambiguous, missing values become zeros, and a spatially unordered file loses the bounding-box read that made the conversion worthwhile. This guide converts a county fabric correctly and proves nothing changed. It is the vector half of cloud-native geospatial formats for compliance pipelines.

Prerequisites

Step-by-step

Step 1: Read with the encoding and CRS stated explicitly

Two things routinely go wrong at read time and both are silent. A shapefile’s attribute encoding is declared in an optional .cpg file that is frequently absent or wrong — the common case is a file declared as UTF-8 that is actually Latin-1, which mangles owner names and street names. And a missing .prj leaves the CRS unset, which every subsequent step then guesses at.

import geopandas as gpd

def read_shapefile(path, encoding="utf-8", expect_epsg=None):
    gdf = gpd.read_file(path, encoding=encoding)
    if gdf.crs is None:
        raise ValueError(f"{path} has no .prj — the CRS must be supplied explicitly")
    if expect_epsg and gdf.crs.to_epsg() != expect_epsg:
        raise ValueError(f"{path} is EPSG:{gdf.crs.to_epsg()}, expected {expect_epsg}")
    return gdf

Where the encoding is uncertain, reading as Latin-1 never raises and never loses bytes, so it is the safer fallback: any mojibake is visible in the output rather than being a decode error that stops the run.

What the shapefile format cannot carryEach format limitation against what it does to a parcel fabric and what the conversion has to restore.Effect on the fabricRestored byTen-character fieldnamesZONING_DIS and ZONING_D_1A rename map from the data dictionaryNo null for numbersMissing encoded as 0 or -9999Sentinel conversion, deliberatelyNo date typeDates as strings or numbersParsing to a real timestampEncoding declaredoptionallyMangled owner and street namesReading with a stated encoding
The conversion is often the step where the data stops being lossy — the shapefile is the lossy representation, not the Parquet.

Step 2: Restore field names and types

The DBF limit truncates field names to ten characters, which turns ZONING_DISTRICT_CODE and ZONING_DISTRICT_NAME into ZONING_DIS and ZONING_D_1. Map them back from the data dictionary before anything downstream references them, and fix the types the format could not carry.

import pandas as pd

RENAMES = {
    "ZONING_DIS": "zoning_district_code",
    "ZONING_D_1": "zoning_district_name",
    "LAST_UPDAT": "last_updated",
    "DWELL_UNIT": "dwelling_units",
}

def restore_schema(gdf, renames=RENAMES, date_cols=("last_updated",),
                   string_cols=("zoning_district_code",), null_sentinels=(-9999, 0)):
    gdf = gdf.rename(columns=renames)
    for c in date_cols:
        if c in gdf:
            gdf[c] = pd.to_datetime(gdf[c], errors="coerce", utc=True)
    for c in string_cols:
        if c in gdf:
            gdf[c] = gdf[c].astype("string").str.strip()
    return gdf

The null_sentinels parameter names the problem it exists for. Shapefiles cannot store a null number, so producers encode missing values as 0 or -9999, and a dwelling unit count of zero is indistinguishable from unknown. Converting these to genuine nulls is a decision that needs the data dictionary and needs recording — turning 0 into null is wrong if the field genuinely means zero units.

Step 3: Sort spatially before writing

This is the step that determines whether the converted file is fast. Parquet stores row groups with per-group statistics, and a bounding-box read skips any row group whose bounds do not intersect the query. That only helps if rows near each other in space are near each other in the file — which, for a shapefile in whatever order the producer happened to write it, they are not.

import numpy as np

def hilbert_sort(gdf):
    """Order rows so that spatial neighbours are file neighbours."""
    return gdf.iloc[gdf.geometry.hilbert_distance().argsort()].reset_index(drop=True)

The effect is large. On a county fabric, an unsorted GeoParquet read for a small area touches nearly every row group, while a Hilbert-sorted one touches a handful — the difference between reading the whole county and reading a tile, which is precisely what chunking county-scale runs by spatial tile depends on.

What the spatial sort buys on a small-area readRow groups touched by a bounding-box read of one tile, from an unsorted file and from a Hilbert-sorted one.Hilbert-sorted, with covering bbox3 row groups readHilbert-sorted, no covering bbox12 row groups readUnsorted — producer order47 row groups readIllustrative, for one tile of a county fabric at 25,000 rows per group. The verification step measures this against your own data.
A conversion that loses the sort passes every correctness check and makes every tiled run several times slower — a regression with no obvious cause.

Step 4: Write with a row group size matched to the read pattern

Row groups are the unit of skipping, so their size trades index granularity against per-group overhead. Very large groups skip nothing useful; very small ones make the metadata large and the reads chatty.

def write_geoparquet(gdf, path, row_group_size=25_000):
    gdf.to_parquet(
        path,
        index=False,
        compression="zstd",
        row_group_size=row_group_size,
        write_covering_bbox=True,     # per-row bbox columns: enables bbox pushdown
        schema_version="1.1.0",
    )

write_covering_bbox is what makes gpd.read_parquet(path, bbox=...) genuinely selective rather than a post-read filter, and zstd compresses parcel attribute data considerably better than snappy at similar read speed. A county fabric that arrives as a 900 MB shapefile set typically lands around 120–180 MB.

Step 5: Prove the conversion changed nothing

A conversion is only useful if it is faithful, and faithfulness is checkable rather than assumed.

def compare(src, dst, key="parcel_id", tol=1e-6):
    problems = []
    if len(src) != len(dst):
        problems.append(f"row count {len(src)} -> {len(dst)}")
    if src.crs != dst.crs:
        problems.append(f"CRS {src.crs} -> {dst.crs}")

    a = src.set_index(key).geometry.sort_index()
    b = dst.set_index(key).geometry.sort_index()
    if not a.index.equals(b.index):
        problems.append("parcel ids differ")
    else:
        area_delta = (b.area - a.area).abs().max()
        if area_delta > tol:
            problems.append(f"max area change {area_delta:.6f}")
        if not b.geom_equals_exact(a, tolerance=tol).all():
            problems.append("geometry differs beyond tolerance")
    return problems

assert not compare(original, gpd.read_parquet(out_path))

Comparing areas as well as geometry catches the case where a coordinate precision setting quietly rounded vertices. The geometry itself should be byte-identical in practice — GeoParquet stores WKB — so any difference at all is worth investigating rather than tolerating.

Step 6: Record the conversion in the lineage

The converted file is what every run reads, so it needs the same provenance treatment as the source.

manifest = {
    "source": {"path": str(shp_path), "sha256": sha256_of(shp_path),
               "encoding": encoding, "crs": str(original.crs)},
    "output": {"path": str(out_path), "sha256": sha256_of(out_path),
               "rows": len(gdf), "row_group_size": 25_000, "sorted": "hilbert"},
    "renames": RENAMES,
    "null_sentinels_converted": [-9999],
    "geoparquet_version": "1.1.0",
}
The conversion is a step, so it gets a recordSource hash, output hash, the rename map and the sentinel decisions are recorded together, because none of them can be reconstructed from the output alone.Hash the sourceplus the encoding actually usedRecord every renamethe part a future reader cannot recoverRecord the sentinel conversionsturning 0 into null is a decisionHash the output, and prove equalityrow count, CRS, geometry, area
A column called zoning_district_code gives no hint that it was ZONING_DIS, and reconciling against the original a year later needs that mapping.

The renames belong in the manifest specifically because they are the part a future reader cannot reconstruct. A column called zoning_district_code gives no hint that it was ZONING_DIS in the source, and reconciling a pipeline output against the original shapefile a year later requires that mapping.

Verification

Beyond the equality check, confirm the file reads back the way the pipeline will read it and that the spatial ordering is doing its job.

import time

full = gpd.read_parquet(out_path)
assert full.crs.to_epsg() == original.crs.to_epsg()
assert full["dwelling_units"].isna().sum() == expected_unknown_count

# A small-area read should be much cheaper than a full read.
t0 = time.perf_counter(); gpd.read_parquet(out_path); t_full = time.perf_counter() - t0
t0 = time.perf_counter(); tile = gpd.read_parquet(out_path, bbox=sample_bbox)
t_bbox = time.perf_counter() - t0
print(f"full {t_full:.2f}s / bbox {t_bbox:.2f}s for {len(tile)} of {len(full)} rows")
assert t_bbox < t_full / 4, "bbox read is not selective — check the spatial sort"

That last assertion is the one worth keeping in CI. A conversion that loses the Hilbert sort still passes every correctness check and quietly makes every tiled run several times slower, which is a regression nobody attributes to the conversion.

Common Pitfalls

  • Trusting the .cpg file. A wrong encoding declaration mangles text silently. Latin-1 as a fallback never raises and leaves the damage visible.
  • Leaving truncated field names in place. ZONING_D_1 is meaningless in six months and unmappable without the data dictionary.
  • Converting sentinel values to null without checking. A dwelling unit count of 0 may be genuine; -9999 never is.
  • Skipping the spatial sort. Everything still works and every tiled read gets slower, which is a performance regression with no obvious cause.
  • Omitting write_covering_bbox. Without it, a bbox= read filters after reading rather than skipping row groups.
  • Assuming multipart geometry survives unchanged. It does, but shapefiles do not distinguish polygon from multipolygon, so a mixed layer needs its geometry type checked rather than declared.

Frequently Asked Questions

Is GeoPackage a better target than GeoParquet?

For interactive editing and for tools that expect OGC formats, yes — GeoPackage is a proper database with indexes and it opens in every desktop GIS. For batch analytics it is not close: columnar storage, predicate pushdown and compression make GeoParquet several times faster for the read-mostly, column-subset access a compliance pipeline performs. Many projects keep both, generated from the same conversion.

How large should row groups be?

Twenty-five thousand rows is a reasonable default for parcels. Tune it against the actual read pattern: if tiles typically pull a few thousand parcels, smaller groups skip more; if the pipeline mostly reads whole counties, larger groups reduce overhead. The measurement in the verification step is the way to choose rather than the guidance.

Does the CRS survive a round-trip?

Yes — GeoParquet stores the CRS as PROJJSON, which is strictly more faithful than a shapefile’s .prj, since a WKT1 .prj cannot express everything a modern CRS definition carries. In practice this is one of the strongest reasons to convert: it is common for the shapefile to be the lossy representation, not the Parquet. Confirming what actually came through is the job of validating CRS metadata before a compliance run.

Should the conversion repair geometry?

No — keep conversion and repair as separate steps with separate records. A conversion that also repairs cannot be verified by the equality check above, and the whole value of that check is that it isolates format change from data change. Repair afterwards, per fixing invalid parcel polygons with make_valid.

Can the file be partitioned as well as sorted?

Yes, and it is worth it above a few million rows: writing a dataset partitioned by county or by a coarse H3 cell lets a reader skip whole directories before touching any Parquet metadata. Below that scale, Hilbert sorting plus row-group statistics does most of the work with none of the directory management.

What about the parcels’ attribute nulls in downstream code?

They become real nulls, which is the point, and it will surface arithmetic that was quietly treating missing as zero. That is a genuine improvement even though it looks like new breakage — a density calculation over parcels with unknown unit counts should not silently report zero units, and after conversion it will raise or produce NaN instead.

Part of: Cloud-native geospatial formats for compliance pipelines