Shapefile & GeoJSON Parsing: Architecting Robust Vector Ingestion

Reliable vector ingestion is the entry point of Geospatial Data Ingestion & Processing Workflows: before any layer can be reprojected, joined, or validated, it has to be read correctly and turned into a typed, CRS-tagged model. This stage bridges the legacy ESRI Shapefile — a 1990s multi-file format that still dominates public data portals — with modern RFC 7946 GeoJSON, and hands a clean GeoDataFrame to downstream steps such as Coordinate Reference System Transformations and Spatial Joins & Merging. The focus throughout is spatial accuracy, memory-efficient parsing, and production-grade pipeline design rather than one-off scripts.

Architecture & Data Structures

A Shapefile is not one file but a family: the mandatory .shp (geometry), .shx (a fixed-length offset index into the .shp), and .dbf (a dBASE IV attribute table), plus the optional-but-critical .prj (WKT projection) and .cpg (attribute encoding). Miss the .prj and the layer arrives with no CRS; miss the .cpg and non-ASCII attribute values decode with the wrong codepage. GeoJSON collapses all of that into a single UTF-8 text document: a FeatureCollection whose features each carry a geometry and a flat properties object, with coordinates fixed by RFC 7946 to WGS84 longitude-latitude order.

Both formats resolve to the same in-memory object — a GeoDataFrame — so the reader is the only format-aware component in the pipeline. Everything downstream operates on the geometry column (a Shapely-backed GeoSeries) and ordinary pandas columns, regardless of where the data came from.

import geopandas as gpd

# One reader, either format — the result type is identical
parcels = gpd.read_file("parcels.shp", engine="pyogrio")
zoning = gpd.read_file("zoning.geojson", engine="pyogrio")

print(type(parcels))            # <class 'geopandas.geodataframe.GeoDataFrame'>
print(parcels.geometry.name)    # 'geometry'  — the active geometry column
print(parcels.crs)              # EPSG:25832  (None if the .prj was missing)
print(parcels.geom_type.value_counts())

The single active-geometry-column model is why this stage is worth isolating: once a layer is a well-formed GeoDataFrame with a known CRS, the rest of the workflow never needs to know it began life as a Shapefile.

Two structural facts about the .shp explain most of the surprises that follow. The first is that a shape type is a property of the file, not of the record: the header declares Point, PolyLine, Polygon or MultiPoint once, and every record must match it. There is consequently no distinct multipolygon type — a Polygon record simply carries a list of rings, and whether a given ring is a second outer part or a hole in the first is decided purely by its winding direction. ESRI's rule is the inverse of the GeoJSON one: outer rings run clockwise, holes counter-clockwise. A converter that copies coordinates without re-orienting them produces a file that opens fine and renders with its islands punched out as holes. The second fact is the 2 GB ceiling. Record offsets in the .shx are stored as 16-bit word counts in a 32-bit field, which caps each component file at roughly two gigabytes; a national parcel layer hits that wall in the .dbf long before the geometry runs out, and the write fails or truncates rather than raising a clean error.

Elevation and measure values ride along as separate shape variants — PointZ, PolygonZM and their siblings — and they are the reason a layer sometimes arrives with a geom_type you expected but coordinates you did not. Shapely 2.0 carries Z through faithfully, drops M silently (GEOS has no measure dimension), and reports the difference only through has_z. Two more sidecars may appear next to the mandatory files: .sbn/.sbx are ESRI's proprietary spatial index, which GDAL ignores entirely, and .qix is the quadtree index GDAL and MapServer do read, turning a bounding-box read into an index seek instead of a full scan. Neither carries data, so both are safe to delete — but deleting the .qix silently costs you every read-time spatial filter described later on this page.

On the GeoJSON side the equivalent structural details are three optional members that no parser validates for you. A bbox may sit on the collection, on a feature, or on a geometry, and is never recomputed on edit — a stale one is a live hazard for any consumer that trusts it for culling. A feature id sits beside properties, not inside it, so a round-trip through a GeoDataFrame loses it unless you ask for it explicitly. And foreign members — arbitrary extra keys at any level — are explicitly permitted and explicitly ignorable, which is how vendor metadata ends up in a file that every reader silently discards.

Environment Configuration & Dependency Resolution

The modern stack favours pyogrio over the older fiona for its multi-threaded GDAL/OGR bindings and Arrow-based memory mapping, which deliver roughly 5–10× read throughput on large vector datasets. pyogrio is the default engine in GeoPandas 1.0+, but pinning it explicitly keeps behaviour reproducible across machines. Install the whole stack from a single channel — mixing pip GDAL wheels with conda-forge builds is the root cause of most PROJ and encoding failures.

# Reproducible single-channel install
conda install -c conda-forge \
  "geopandas>=1.0" "pyogrio>=0.8" "shapely>=2.0" "pyproj>=3.6"

import geopandas as gpd, pyogrio, shapely, pyproj

for mod in (gpd, pyogrio, shapely, pyproj):
    print(mod.__name__, mod.__version__)

# Confirm GDAL sees the Shapefile and GeoJSON drivers as writable
drivers = pyogrio.list_drivers()
assert drivers.get("ESRI Shapefile") == "rw"
assert drivers.get("GeoJSON") == "rw"

Shapely 2.0 is a hard requirement, not a nicety: its vectorized, GEOS-backed geometry engine is what makes bulk validity checks and repairs fast enough to run inside an ingestion loop rather than as an offline batch. If shapely.__version__ reports 1.x, geometry operations silently fall back to slow, per-object Python loops.

Which engine is active changes behaviour, not just speed, so pin it rather than inferring it. GeoPandas 1.0 defaults to pyogrio and falls back to fiona only if pyogrio is absent; on 0.x the default is the reverse. Set it once at import time — gpd.options.io_engine = "pyogrio" — and every read_file in the process is deterministic regardless of what is installed. The two engines differ in ways that surface as data bugs rather than errors: fiona hands GDAL fields back one feature at a time as Python objects and preserves None for missing values, while pyogrio builds whole columns at once and represents a missing numeric as NaN, which promotes the column to float. Layer-open semantics differ too — fiona's read_file accepts a zip:// URI, pyogrio wants GDAL's /vsizip/ prefix — and only pyogrio exposes the record-batch iterator used below.

import geopandas as gpd
from osgeo import gdal   # only if the GDAL Python bindings are installed

gpd.options.io_engine = "pyogrio"           # deterministic, regardless of install order

# GDAL's own version gates several export options used later on this page
print(gdal.VersionInfo("RELEASE_NAME"))     # e.g. 3.8.4

Three capability gates are worth checking before a pipeline is written rather than after it fails in staging. The RFC7946=YES GeoJSON creation option needs GDAL 2.4 or newer; the GeoJSONSeq newline-delimited driver arrived in the same generation; and Arrow-accelerated reads (gpd.read_file(..., use_arrow=True)) require both pyogrio 0.8+ and a GDAL built with Arrow support, falling back silently to the ordinary path when either is missing. Mixing installation channels is what makes these gates unpredictable — a pip install gdal wheel next to a conda-forge libgdal gives you two GDAL libraries in one process, and the one that loaded first wins.

Vectorized Operations & Core Workflow

The canonical ingestion path reads the source, enforces a CRS, and filters columns early to keep memory flat. Passing columns=[...] to the reader pushes projection down into GDAL so unused attribute columns are never materialised — a large win on wide .dbf tables.

import geopandas as gpd

# High-throughput ingestion with explicit CRS enforcement.
# `columns` prunes attributes at the GDAL layer, before any Python objects exist.
parcels = gpd.read_file(
    "parcels.shp",
    engine="pyogrio",
    columns=["parcel_id", "owner", "area_m2", "geometry"],
)

# A .prj-less Shapefile reads back as crs=None. Backfill the KNOWN authority
# code — set_crs only labels, it never moves coordinates.
if parcels.crs is None:
    parcels = parcels.set_crs("EPSG:25832")   # ETRS89 / UTM 32N — a metric CRS

print(parcels.geom_type.value_counts())
print(f"{len(parcels):,} features, {parcels.memory_usage(deep=True).sum() / 1e6:.1f} MB")

For datasets larger than RAM, read in record batches instead of one call. pyogrio.read_dataframe streams fixed-size chunks so peak memory stays proportional to the batch, not the file.

import pyogrio
import geopandas as gpd

batches = []
for batch in pyogrio.read_dataframe(
    "national_parcels.shp",
    columns=["parcel_id", "area_m2", "geometry"],
    batch_size=50_000,
    return_fids=False,
    as_iterator=True,
):
    # process each 50k-row GeoDataFrame here (filter, repair, append)
    batches.append(batch[batch["area_m2"] > 0])

parcels = gpd.GeoDataFrame(
    __import__("pandas").concat(batches, ignore_index=True),
    crs="EPSG:25832",
)

Batching bounds memory, but the bigger win is not reading the rows at all. GDAL evaluates three independent filters below the Python layer, and each one discards features before a single Shapely object is constructed: bbox tests the record's stored bounding box, mask tests against a real geometry, and where is an OGR SQL predicate over the .dbf columns. On a Shapefile with a .qix index the bounding-box test becomes an index seek rather than a scan, which is the difference between reading four hundred features and reading four million.

import geopandas as gpd
from shapely.geometry import box

# Study area expressed in the LAYER's CRS — pyogrio does not reproject a raw tuple
study_bbox = (390000, 5810000, 405000, 5825000)   # EPSG:25832 metres

downtown_parcels = gpd.read_file(
    "national_parcels.shp",
    engine="pyogrio",
    bbox=study_bbox,                          # driver-level extent filter
    where="area_m2 > 500 AND land_use <> 'ROAD'",   # OGR SQL over the .dbf
    columns=["parcel_id", "land_use", "area_m2"],
)

# A GeoSeries mask is reprojected for you; a bare tuple is NOT.
floodplain_boundary = gpd.read_file("floodplain.geojson").to_crs("EPSG:25832")
at_risk = gpd.read_file(
    "national_parcels.shp",
    engine="pyogrio",
    mask=floodplain_boundary.geometry,        # exact geometry, not just its envelope
)

print(f"{len(downtown_parcels):,} of the layer's features were materialised")

The distinction between bbox and mask costs people real time. bbox is a rectangle test against stored envelopes and is nearly free; mask intersects against a geometry and therefore reads and tests candidate records, so it is exact but slower. Chaining them — a bbox to narrow, then a Shapely predicate in memory to refine — is usually faster than a mask alone. The trap in both is the coordinate system: a raw tuple is interpreted in the layer's CRS with no conversion, so passing degrees to a UTM layer returns an empty frame rather than an error. Passing a GeoSeries instead lets GeoPandas reproject the filter for you, which is why the second call above is the safer habit.

For sampling a large delivery — the first thing anyone does with an unfamiliar file — skip_features and max_features page through it without touching the rest, and pyogrio.read_info answers the structural questions (feature count, field names and types, CRS, geometry type) by reading headers only.

import pyogrio

info = pyogrio.read_info("national_parcels.shp")
print(info["features"], info["geometry_type"], info["crs"])
# 4127653 Polygon EPSG:25832

sample = gpd.read_file("national_parcels.shp", engine="pyogrio",
                       skip_features=0, max_features=1000)

Geometry & Data Processing Details

Raw vector inputs routinely carry self-intersections, duplicate vertices, ring-order errors, and null geometries that break spatial predicates the moment they are used. Repair them at ingestion with shapely.make_valid, which uses GEOS to return a valid geometry of the appropriate type rather than throwing. This is the same repair logic covered in depth under Topology Validation & Repair; doing a first pass here means malformed records never reach the join stage.

from shapely import make_valid
import geopandas as gpd

# 1. Repair invalid geometries with GEOS-backed validation (vectorized in Shapely 2.0)
mask_invalid = ~parcels.geometry.is_valid
parcels.loc[mask_invalid, "geometry"] = parcels.loc[mask_invalid, "geometry"].apply(make_valid)

# 2. Drop null or empty geometries that survive repair
parcels = parcels[parcels.geometry.notna() & ~parcels.geometry.is_empty].copy()

# 3. Build the spatial index once so downstream bbox queries and joins are fast
_ = parcels.sindex        # triggers R-tree construction

assert parcels.geometry.is_valid.all(), "invalid geometry survived repair"

Two format-specific processing details matter here. First, the Shapefile .dbf table cannot store datetime, bool, or full float64 precision — dates come back as strings or get truncated, so cast temporal fields to ISO-8601 explicitly and round floats before you rely on them. Second, Shapefiles freely mix single- and multi-part geometries in one layer; explode them to a uniform part type when a downstream step expects singleparts.

# Coerce mixed single/multi geometries to a uniform part type
if parcels.geom_type.nunique() > 1:
    parcels = parcels.explode(index_parts=False).reset_index(drop=True)

# Shapefile field names are capped at 10 characters — normalise before re-export
parcels.columns = [c[:10] if c != "geometry" else c for c in parcels.columns]

For polygon repair specifically, prefer make_valid over the old buffer(0) trick: a zero-width buffer can silently alter area and drop interior rings, whereas make_valid preserves the original geometry's structure. Use explain_validity to log why a record failed before deciding how to handle it.

A third detail is dimensionality, and it is the one that survives longest undetected. A PolygonZ layer read into GeoPandas keeps its Z ordinate, so every coordinate is a triple; .area and .length ignore the third value, to_file writes it back, and a GeoJSON export emits three-element positions that a web renderer either ignores or misreads as a styling hint. Unless the elevation is part of the analysis, flatten on read with pyogrio's force_2d=True rather than discovering the extra ordinate three steps downstream. Measure (M) values need no such flag because GEOS has no measure dimension at all — they are dropped on parse, permanently, which means a Shapefile is the only place that data still exists.

import geopandas as gpd
import shapely

parcels = gpd.read_file("parcels_3d.shp", engine="pyogrio", force_2d=True)
assert not shapely.has_z(parcels.geometry.values).any(), "Z survived the read"

# Ring orientation: RFC 7946 wants exterior rings counter-clockwise; Shapefile
# writes them clockwise. Normalise before any GeoJSON export.
# shapely.ops.orient handles Polygon and MultiPolygon; sign=1.0 means CCW exterior
from shapely.ops import orient

parcels["geometry"] = parcels.geometry.apply(lambda g: orient(g, sign=1.0))

Attribute typing has one more Shapefile-specific trap that no repair pass catches. The .dbf stores every value as fixed-width ASCII with a declared width and precision, and it has no concept of NULL — a missing number is written as blanks, which readers interpret as zero, empty string, or NaN depending on the field type and the driver. A population column where "unknown" and "zero" are the same byte pattern cannot be repaired after the fact; the only fix is to carry an explicit sentinel or a companion "is_known" flag into the file, or to stop using the format for anything with genuine nulls. Width matters too: a double field declared with width 10 and precision 2 truncates rather than rounds on write, so an area of 12345678.91 comes back as a different number than the one you computed.

CRS Alignment & Projection Pipeline

The two formats disagree about coordinates by design, and that disagreement is the single most common source of silent errors at this stage. GeoJSON is fixed by RFC 7946 to WGS84 (EPSG:4326) in longitude, latitude order; a Shapefile can be in any CRS its .prj declares — or none at all. Reconcile both to one authoritative projection before measuring anything, using the transformation rules detailed under Coordinate Reference System Transformations and, at the library level, Coordinate Systems with PyProj.

import geopandas as gpd

# GeoJSON is ALWAYS lon/lat WGS84 by spec — read it, then reproject for metric work
zoning = gpd.read_file("zoning.geojson", engine="pyogrio")   # EPSG:4326
assert zoning.crs.to_epsg() == 4326

# Reproject to a LOCAL metric CRS before any area/distance/overlay.
# Never do metric analysis in EPSG:4326 (degrees) or EPSG:3857 (scale distorts with latitude).
zoning_m = zoning.to_crs("EPSG:25832")     # ETRS89 / UTM 32N — units are metres
assert zoning_m.crs.axis_info[0].unit_name == "metre"

# Align the Shapefile-derived layer to the SAME target so a join can match
parcels_m = parcels.to_crs("EPSG:25832")

The .prj itself deserves scrutiny before it is trusted. ESRI writes WKT1 with its own spellings — GCS_North_American_1983, Lambert_Conformal_Conic, a Transverse_Mercator with no authority code anywhere in the string — and pyproj parses it into a perfectly usable CRS object that cannot name itself. crs.to_epsg() then returns None, which breaks every downstream assertion written as crs.to_epsg() == 25832 and every database write that needs an SRID. The fix is to lower the matching threshold rather than to hard-code a guess: pyproj scores candidate registry entries and to_epsg(min_confidence=25) accepts a good-but-not-exact match, which is almost always the right one for a national grid.

import geopandas as gpd
from pyproj import CRS

parcels = gpd.read_file("vendor_delivery.shp", engine="pyogrio")

code = parcels.crs.to_epsg()                     # None on an ESRI WKT1 .prj
if code is None:
    code = parcels.crs.to_epsg(min_confidence=25)   # fuzzy registry match

print(code, parcels.crs.name, parcels.crs.is_projected)
# 25832 ETRS89 / UTM zone 32N True

# Re-tag with the authority code so the SRID is explicit downstream
parcels = parcels.set_crs(CRS.from_epsg(code), allow_override=True)

Confirm the match before accepting it — compare parcels.crs.name and the datum against the provider's metadata, because a confident-looking fuzzy match between two datums that differ by a metre-scale shift is exactly the error a bounds check cannot see. When the two candidates differ only by datum, the grid-shift mechanics that decide the answer are covered in Coordinate Systems with PyProj.

Two rules prevent the classic failures. Use set_crs only to label a genuinely unlabelled layer (the .prj-less Shapefile) and to_crs to actually move coordinates — reaching for set_crs to fix visible misalignment just relabels wrong data as right. And when you build a raw pyproj.Transformer for point work, pass always_xy=True, because EPSG:4326 is defined latitude-first by its authority while Python, GeoJSON, and web maps all assume longitude-first; omitting it swaps every coordinate and lands your data in the ocean off West Africa.

Production Export & Integration

Export format should match the consumer. For analytical and cloud pipelines, serialise to GeoParquet — it is columnar, typed, compressed, and preserves the CRS, unlike the Shapefile's lossy .dbf. The comparison and trade-offs live in GeoParquet vs Shapefile for Storage under Cloud-Native Geospatial Formats. For web mapping, emit lean RFC 7946 GeoJSON in EPSG:4326 with trimmed coordinate precision and stripped null properties.

import json
import geopandas as gpd
from shapely import set_precision

# --- Analytical target: GeoParquet (keeps CRS + dtypes, ~10x smaller than Shapefile) ---
parcels_m.to_parquet("parcels.parquet")

# --- Web target: optimised RFC 7946 GeoJSON ---
web = parcels_m.to_crs("EPSG:4326").copy()

# Snap coordinates to ~1e-6 deg (~11 cm) to shrink payload without visible loss
web["geometry"] = set_precision(web.geometry.values, grid_size=1e-6)

feature_collection = json.loads(web.to_json(drop_id=True))
for feature in feature_collection["features"]:
    feature["properties"] = {k: v for k, v in feature["properties"].items() if v is not None}

with open("parcels_web.geojson", "w", encoding="utf-8") as fh:
    json.dump(feature_collection, fh, separators=(",", ":"))   # no whitespace

Hand-building the GeoJSON as above gives you exact control, but GDAL will do the same work at the driver level if you let it, and its options are worth knowing because they are enforced rather than advisory. RFC7946=YES makes the writer reproject to WGS84, emit right-hand-rule ring orientation, split geometries at the antimeridian, and round coordinates to seven decimal places — the whole conformance checklist in one flag. COORDINATE_PRECISION overrides that rounding, and WRITE_BBOX adds the collection-level bounding box that some tile builders expect.

import geopandas as gpd

web = parcels_m.to_crs("EPSG:4326")

# Driver-enforced conformance instead of hand-rolled post-processing
web.to_file(
    "parcels_web.geojson",
    driver="GeoJSON",
    engine="pyogrio",
    layer_options={"RFC7946": "YES", "COORDINATE_PRECISION": "6", "WRITE_BBOX": "YES"},
)

# Newline-delimited: one feature per line, streamable and appendable
web.to_file("parcels_stream.geojsonl", driver="GeoJSONSeq", engine="pyogrio")

The newline-delimited variant is the underrated one. Because each line is an independent feature, a consumer can process the file with a for line in fh loop at constant memory, a producer can append to it without rewriting the document, and a failure part-way through leaves a truncated but still parseable file rather than an unclosed bracket. It is not RFC 7946 — there is no FeatureCollection wrapper — so it belongs in pipelines, not in browser downloads. Where the consumer needs both streamability and a spatial index, FlatGeobuf is the better single-file target: one binary file, a packed Hilbert R-tree in the header, and partial reads over HTTP range requests, with none of the Shapefile's sidecar coupling or field-name ceiling.

Expect roughly an order of magnitude between these targets on the same data. A million-parcel layer that occupies about 250 MB across a Shapefile family lands near 400 MB as pretty-printed GeoJSON, around 180 MB as minified GeoJSON at six-decimal precision, and well under 40 MB as compressed GeoParquet — and coordinate precision, not the container, is the single biggest lever on the text formats. Seven decimal places of longitude resolves about a centimetre; parcel boundaries surveyed to the nearest decimetre carry four digits of pure noise per ordinate, and trimming them is lossless in every sense that matters.

A short pre-export checklist keeps outputs trustworthy: assert the CRS is what you claim (gdf.crs.to_epsg()), confirm gdf.geometry.is_valid.all(), and for anything over ~50 MB destined for a browser, reach for tiled PMTiles or a PostGIS/GeoParquet backend instead of a monolithic GeoJSON download. When loading into PostGIS, write with gdf.to_postgis("parcels", engine, if_exists="replace") after confirming the geometry column's SRID matches the target table.

Windows & Platform Edge Cases & Debugging

Most parsing failures on this stage are environmental or encoding-related, not logical, and fall into a few recurring causes:

• Non-ASCII attributes come back as mojibake (ü, é). The .dbf was written in a legacy codepage with no .cpg sidecar. Read with the correct encoding explicitly — gpd.read_file("legacy.shp", encoding="cp1252") (or latin1) — rather than trusting the default utf-8.
• .prj-less Shapefile reads back as crs=None. GeoPandas cannot infer the CRS; set_crs the known authority code as the very first step, before any reprojection or join.
• missing .shx file / unable to open .shx. The offset index was lost in a partial copy. Set SHAPE_RESTORE_SHX=YES in the environment (os.environ["SHAPE_RESTORE_SHX"] = "YES") so GDAL rebuilds it from the .shp.
• CRSError or missing PROJ grids on Windows. The PROJ_DATA (older builds: PROJ_LIB) variable points at a stale directory when GDAL, Fiona, and pyproj come from mixed channels. Confirm with pyproj.datadir.get_data_dir() and reinstall the stack from one conda-forge environment.
• GeoJSON coordinates look transposed / points fall in the wrong hemisphere. The document was written [lat, lon] in violation of RFC 7946, or a Transformer was built without always_xy=True. Validate structure first — the recipe is Best Practices for GeoJSON Validation.
• Field names silently truncated or duplicated on Shapefile export. Two columns collide once clipped to 10 characters. Rename before writing, or export to GeoParquet/GeoJSON where the 10-char cap does not apply. The full automation is in Automating Shapefile Cleanup with Python.
• Reading a zipped delivery fails or finds no layer. Do not unzip it — read through GDAL's virtual filesystem: gpd.read_file("/vsizip/deliveries/parcels.zip/parcels.shp", engine="pyogrio"). The inner path is required when the archive holds more than one Shapefile family, and the fiona-style zip:// prefix is not interchangeable with /vsizip/.
• .cpg says one thing, the bytes say another. Vendor .cpg files are frequently wrong. The GDAL config option SHAPE_ENCODING overrides it for the whole process (os.environ["SHAPE_ENCODING"] = "cp1252"), and setting it to the empty string disables recoding entirely so you get raw bytes to inspect. Decide from a sample of known accented values, not from the sidecar.
• PermissionError writing a Shapefile on Windows. Another process holds the .shp or .dbf open — a desktop GIS session, a file-preview handler, or an antivirus scanner mid-scan. Windows locks are mandatory, unlike the advisory locks on Linux, so the same script that works in CI fails on a workstation. Write to a temporary name and rename afterwards.
• Paths break past 260 characters or on non-ASCII directory names. Windows' legacy MAX_PATH limit still bites deeply nested output trees, and GDAL expects UTF-8 filenames regardless of the console codepage. Keep ingestion roots short and prefer ASCII directory names for anything a driver will open by path — the failure surfaces as a file-not-found on a file that plainly exists.
• A GeoJSON with 3D coordinates renders flat but joins wrong. The Z ordinate survived the read and is being carried through predicates that ignore it; re-read with force_2d=True so the geometry matches what the analysis assumes.

Frequently Asked Questions

Should I use pyogrio or fiona to read Shapefiles? Use pyogrio for anything but the smallest files — its multi-threaded, Arrow-backed reads are 5–10× faster and it is the GeoPandas 1.0+ default. Keep fiona only if you need record-by-record streaming with per-feature Python callbacks that pyogrio's batch model does not expose.

Why does my Shapefile lose its coordinate system after reading? The dataset shipped without a .prj file, so there is no CRS to read. GeoPandas returns crs=None; you must set_crs the known authority code (from the data provider's metadata) before reprojecting. Guessing EPSG:4326 when the coordinates are actually a projected grid will silently corrupt every later measurement.

How do I fix latin-1/cp1252 garbled text in Shapefile attributes? Pass the source encoding explicitly: gpd.read_file("data.shp", encoding="cp1252"). The garbling happens because the .dbf lacks a .cpg sidecar declaring its codepage, so the reader defaults to UTF-8 and mis-decodes accented characters. Re-export with a .cpg or move to GeoJSON/GeoParquet, which are UTF-8 by definition.

Is GeoJSON always in WGS84? Per RFC 7946, yes — coordinates must be WGS84 (EPSG:4326) in longitude, latitude order. Older pre-2016 GeoJSON with a crs member is non-conformant; treat its declared CRS as advisory, verify against known control points, and reproject to a local metric CRS before any area or distance work.

When should I convert a Shapefile to GeoParquet instead of keeping it? Whenever the data is analytical, larger than a few hundred MB, or needs to preserve real dtypes and precision. GeoParquet is columnar, compressed, CRS-aware, and free of the Shapefile's 10-character field names and 2 GB size ceiling. Keep Shapefiles only for interchange with tools that cannot read anything else.

How do I keep memory flat when parsing a national-scale Shapefile? Read in batches with pyogrio.read_dataframe(..., batch_size=50_000, as_iterator=True), prune columns with the columns argument so unused attributes never load, and filter or repair each batch before concatenating. Peak memory then tracks one batch, not the whole file.