geopandas-consistency

GeoPandas — consistent, modern idioms

GeoPandas is stable, but generated geospatial code has a signature failure: doing math in degrees. Models also confuse declaring a CRS with reprojecting, loop over geometries that Shapely 2.x handles vectorized, and use pre-1.0 argument names. This skill pins the canonical idiom set — GeoPandas 1.x with Shapely 2.x — built around CRS discipline.

Canonical idioms — always X, never Y

Always	Never	Why
gdf = gdf.set_crs(4326) to declare, gdf = gdf.to_crs(3857) to reproject	using one for the other	set_crs relabels metadata (coords unchanged); to_crs transforms coordinates. Confusing them silently shifts the world.
reproject to a projected CRS (gdf.to_crs(gdf.estimate_utm_crs())) before .area, .length, .buffer, .distance	metric ops in EPSG:4326	Degrees aren't meters: areas/buffers in geographic CRS are dimensionally wrong (GeoPandas warns; the numbers are garbage).
gpd.sjoin(a, b, predicate="within")	op="within"	op= was removed; predicate= is the modern name.
check a.crs == b.crs (or to_crs(b.crs)) before joins/overlays	joining frames with mismatched CRS	Mismatch either raises or — with one side unset — matches nothing/nonsense.
vectorized GeoSeries ops (gdf.geometry.buffer(100), .centroid, .union_all())	gdf.geometry.apply(lambda g: g.buffer(100)) loops	Shapely 2.x is vectorized C; apply-loops are 10–100× slower for identical results.
union_all()	cascaded_union, deprecated unary_union habits	cascaded_union is removed (Shapely 2); union_all() is the modern spelling.
gpd.points_from_xy(df.lon, df.lat)	Point(lat, lon) ordering	Coordinates are (x=lon, y=lat); swapping puts Paris in the Indian Ocean — and nothing errors.
gdf.to_file("out.gpkg") (GeoPackage) as working format	shapefile round-trips by default	Shapefiles truncate column names to 10 chars, mangle dtypes/encodings, and are multi-file.
dissolve(by="region", aggfunc=...)	groupby + manual geometry unioning	dissolve is the canonical group-and-union with attribute aggregation.
keep one active geometry column; set_geometry when switching	overwriting gdf["geometry"] with mixed leftovers	Multiple geometry columns are allowed; ops use only the active one.

House style — "stores within 5 km of a depot":

import geopandas as gpd

stores = gpd.read_file("stores.gpkg")
depots = gpd.read_file("depots.gpkg").to_crs(stores.crs)

metric = stores.estimate_utm_crs()
stores_m = stores.to_crs(metric)
depots_m = depots.to_crs(metric)

near = gpd.sjoin(
    stores_m,
    depots_m.assign(geometry=depots_m.buffer(5_000)),  # meters, because projected
    predicate="within",
    how="inner",
)

Pitfalls that produce silently wrong results

• Area/length in EPSG:4326 returns square-degrees — values look plausible scaled, and downstream math is silently wrong. For global data, an equal-area CRS (e.g. EPSG:6933) for areas; UTM for local work; geographic_crs warning is a hard stop.
• sjoin duplicates rows when geometries match multiple partners (and appends index_right); de-dup or aggregate deliberately. Run a second sjoin only after dropping the leftover index_right column (name collision raises).
• sjoin_nearest max_distance is in CRS units — degrees if you forgot to project, which silently changes the search radius by ~111 km per unit.
• overlay vs sjoin: sjoin attaches attributes (geometries unchanged); overlay cuts geometries (intersection/union/difference). Using sjoin where clipping was intended keeps full polygons and inflates areas.
• Invalid geometries (self-intersections) make predicates lie or raise: gdf.geometry = gdf.geometry.make_valid() before overlay-type work.
• .centroid of multipolygons can fall outside the shape — use representative_point() for label/contains purposes.
• Reading CSVs of lat/lon: build with points_from_xy then set_crs(4326) — data read as plain pandas has no CRS at all, and forgetting set_crs poisons every later to_crs.
• Mixed geometry types in one column (points + polygons) break many ops quietly — gdf.geom_type.value_counts() first.

Version notes

Target GeoPandas 1.x (idioms identical in late 0.x) on Shapely 2.x and pyogrio I/O (default engine in 1.0, much faster than fiona). Breaking lines: op=→predicate=, Shapely 2 removing cascaded_union and making scalar geometry ops vectorized, unary_union property → union_all() method, append-era pandas idioms gone (it's a pandas subclass — pandas-consistency rules apply too).

Workflow

1. On load, immediately check gdf.crs, gdf.geom_type.value_counts(), and validity; declare missing CRS with set_crs only when you know what the coords are.
2. Plan the CRS per operation: geographic for storage/interchange (4326), projected (UTM/equal-area) for measurement; reproject explicitly at the boundary.
3. Use vectorized GeoSeries methods; sjoin/sjoin_nearest/overlay/dissolve/clip for set operations — chosen by whether geometries must be cut.
4. Write GeoPackage for work products; GeoJSON (always 4326 per spec) for web interchange.
5. When reviewing, hunt: set_crs/to_crs confusion, math-in-degrees, op=, apply-loops, Point(lat, lon), shapefile round-trips, unchecked sjoin row explosion.

For the fuller CRS guide, join/overlay matrix, validity handling, and I/O engine notes, read references/geopandas-patterns.md.