Geospatial Analysis skill for Codex

Correct spatial work knows the CRS of every layer, measures in a projected or geodesic system rather than degrees, validates geometries before operating on them, checks that spatial joins produced the expected number of matches, and uses a spatial index for anything beyond a few thousand features.

Inspect first

1. Read AGENTS.md for the stack (GeoPandas, PostGIS, DuckDB spatial) and data location.
2. For every layer, print CRS, geometry types, bounds, count and validity:

import geopandas as gpd

g = gpd.read_file("data/districts.gpkg")          # also .shp, .geojson; GeoParquet: gpd.read_parquet
print(g.crs, g.geom_type.value_counts().to_dict(), g.total_bounds, len(g))
print("invalid:", (~g.is_valid).sum(), "empty:", g.is_empty.sum(), "null:", g.geometry.isna().sum())

ogrinfo -so -al data/districts.gpkg     # GDAL: layer CRS, extent, field types, without Python

3. Bounds tell you the units: values within -180..180 and -90..90 are degrees (likely EPSG:4326); values in the hundreds of thousands or millions are metres in a projected CRS.
4. Ask for or find the source CRS when a file has none (g.crs is None). Never guess silently; a Shapefile without .prj needs the provider to say what it is.

CRS rules

Task	CRS
Storage and exchange (GeoJSON, APIs, web maps)	EPSG:4326 (lon, lat order in GeoJSON)
Web map tiles display	EPSG:3857 (display only, never for measurement)
Areas, buffers, distances within a city, region or country	Local projected CRS: UTM zone (g.estimate_utm_crs()) or the national grid (EPSG:27700 UK, EPSG:25832 Germany, EPSG:2154 France)
Equal-area statistics over continents	An equal-area CRS such as EPSG:6933 or a Lambert azimuthal equal-area
Distances across continents	Geodesic: PostGIS geography, or pyproj.Geod

g = g.set_crs("EPSG:4326")                 # declare a missing CRS (no coordinate change)
g_m = g.to_crs(g.estimate_utm_crs())       # transform coordinates for measurement
g_m["area_km2"] = g_m.area / 1e6

set_crs labels; to_crs transforms. Using set_crs to “convert” is a classic bug that moves data to the wrong continent. Always bring both layers to the same CRS before any overlay or join; GeoPandas raises or warns on mismatch, do not suppress it.

Fix geometries first

import shapely
from shapely import make_valid

bad = ~g.is_valid
print(pd.Series(shapely.is_valid_reason(g.loc[bad].geometry.values)).value_counts().head())
g.loc[bad, "geometry"] = make_valid(g.loc[bad, "geometry"].values)
g = g[~g.is_empty & g.geometry.notna()]

The reasons (self-intersection, ring orientation) tell you whether the source is broken or merely messy. make_valid can return GeometryCollections; filter to polygons with g.explode() plus a type filter when an operation requires polygons. Avoid the old buffer(0) trick on large datasets: it can drop parts of bow-tie polygons.

Spatial joins

points = gpd.GeoDataFrame(
    stores, geometry=gpd.points_from_xy(stores["lon"], stores["lat"]), crs="EPSG:4326"
).to_crs(districts.crs)

joined = gpd.sjoin(points, districts[["district_id", "geometry"]], how="left", predicate="within")
print("unmatched:", joined["district_id"].isna().sum())
print("multi-matched:", joined.index.duplicated().sum())

Need	Call
Point in polygon	sjoin(predicate="within") (points left)
Polygons touching or overlapping	sjoin(predicate="intersects")
Nearest feature, optionally within a radius	sjoin_nearest(right, how="left", max_distance=5000, distance_col="dist_m") in a projected CRS
Area-weighted share of overlap	gpd.overlay(a, b, how="intersection") then area / original area
Merge polygons by attribute	g.dissolve(by="region", aggfunc={"population": "sum"})
Catchment within N metres	g_m.buffer(1000) in metres CRS, then join

• Points on shared boundaries match two polygons with intersects; within excludes boundary points. Decide the tie rule and de-duplicate explicitly.
• sjoin_nearest in EPSG:4326 measures in degrees and warns; project first.
• Swapped lat/lon puts points in the ocean or Antarctica; plot or check bounds.

PostGIS

create index if not exists idx_stores_geom on stores using gist (geom);
analyze stores;

-- stores within 2 km of each customer, metres via geography
select c.customer_id, s.store_id,
       st_distance(c.geom::geography, s.geom::geography) as dist_m
from customers c
join stores s
  on st_dwithin(c.geom::geography, s.geom::geography, 2000);

-- point in polygon with index use
select d.district_id, count(*)
from events e join districts d on st_intersects(d.geom, e.geom)
group by 1;

• Use ST_DWithin, not ST_Distance(...) < x, so the GiST index is used.
• Mixed SRIDs fail; transform with ST_Transform(geom, 25832) or store one SRID per column (geometry(Point, 4326)).
• Very large polygons (countries, coastlines) make joins slow: split with ST_Subdivide(geom, 256) into a helper table and index it.
• ST_IsValid / ST_MakeValid before overlays; check ST_IsValidReason for the cause.
• Creating indexes or helper tables on a shared production database is a schema change; confirm with the user and prefer create index concurrently on busy tables.

Formats

Format	Use
GeoParquet	Default for analysis and large data: fast, typed, keeps CRS (to_parquet)
GeoPackage	Single-file exchange with GIS desktop users (QGIS), multiple layers
GeoJSON	Web and APIs; EPSG:4326 only per RFC 7946; large files are slow
Shapefile	Legacy input only: 10-char field names, 2 GB limit, sidecar files

Reading remote files or tiles needs network; in a sandbox with network off, ask for local copies.

Gotchas

• Area or length in EPSG:4326 returns square degrees; GeoPandas warns. Project first.
• Web Mercator inflates areas toward the poles (Scandinavia, Canada) by large factors.
• Spatial joins can multiply rows; aggregate after joining or de-duplicate by the left index.
• dissolve on invalid geometries fails or produces slivers; fix validity first.
• Precision: coordinates with 15 decimals are false precision; 6 decimals in degrees is about 0.1 m.
• Census and admin boundaries change between vintages; record the boundary year.

Verification

• Plot each layer over a basemap or each other (ax = a.plot(); b.boundary.plot(ax=ax)) and save the PNG; misplaced layers are obvious.
• Unmatched and multi-matched counts after every join are reported and explained.
• One area or distance checked against an external reference (official district area, a measured route) within a few percent.
• Total of a dissolved attribute equals the total before dissolving.

Report template

**Layers:** stores (1,284 points, EPSG:4326), districts (412 polygons, EPSG:25832, 2024 boundaries)
**Fixes:** 7 invalid district polygons repaired with make_valid; 3 stores with swapped lat/lon corrected
**Join:** stores within districts: 1,271 matched, 13 outside all districts (listed in outputs/unmatched.csv), 0 on boundaries
**Measures:** district areas in EPSG:25832; Mitte = 39.5 km2 vs official 39.47 km2
**Outputs:** outputs/stores_by_district.parquet (GeoParquet), figures/coverage.png