Geospatial Analysis skill for Claude Code

Spatial work is correct when every layer has a declared CRS, measurements (distance, area, buffer) happen in a projected CRS appropriate to the region, joins use a spatial index, invalid geometries are repaired or reported rather than dropped silently, and the output map was looked at.

CRS decisions

Task	CRS
Storage, interchange, web maps, GPS input	EPSG:4326 (lon/lat degrees)
Distances, buffers within a city/region	local UTM zone (gdf.estimate_utm_crs()) or national grid (EPSG:27700 UK, EPSG:2154 France)
Areas and densities across large regions	equal-area (EPSG:6933 global, EPSG:5070 CONUS, EPSG:3035 Europe)
Web tiles display	EPSG:3857 - never measure in it (distorts distance and area heavily away from the equator)
Global great-circle distances	geodesic (PostGIS geography, pyproj.Geod)

Axis order trap: EPSG:4326 is formally lat/lon, but GeoJSON, Shapely, GeoPandas, and PostGIS use x=lon, y=lat. When points land in the ocean off Africa or in Antarctica, suspect swapped coordinates or (0, 0) placeholders.

Workflow

- [ ] 1. Load and inspect each layer: CRS, geometry types, bounds, null/empty/invalid counts
- [ ] 2. Set missing CRS (only if known) and reproject to a common CRS
- [ ] 3. Repair invalid geometries; record how many
- [ ] 4. Do the operation (join, buffer, overlay, aggregate) in the right CRS
- [ ] 5. Sanity-check results against known facts (counts, totals, a spot check on a map)
- [ ] 6. Write outputs in an appropriate format and CRS; render the map

1-3. Load, inspect, repair

import geopandas as gpd
import pandas as pd

zones = gpd.read_file("zones.gpkg")                    # pyogrio engine by default in GeoPandas 1.x
pts = pd.read_csv("stores.csv")
stores = gpd.GeoDataFrame(pts, geometry=gpd.points_from_xy(pts.lon, pts.lat), crs="EPSG:4326")

for name, g in {"zones": zones, "stores": stores}.items():
    print(name, g.crs, g.geom_type.value_counts().to_dict(), g.total_bounds,
          "null", g.geometry.isna().sum(), "empty", g.is_empty.sum(), "invalid", (~g.is_valid).sum())

zones = zones.to_crs(stores.crs) if zones.crs != stores.crs else zones
bad = ~zones.is_valid
zones.loc[bad, "geometry"] = zones.loc[bad, "geometry"].make_valid()

• set_crs declares what the coordinates already are; to_crs transforms them. Calling set_crs to “convert” corrupts the data.
• Never guess a missing CRS. Check the source metadata (.prj, documentation); bounds in the hundreds of thousands suggest a projected CRS, -180..180 / -90..90 suggests degrees.
• make_valid can change geometry types (polygon to GeometryCollection); filter or explode afterwards if downstream expects polygons. buffer(0) is the older trick and can drop parts.
• Check bounds against the expected region; drop or flag points at (0, 0).

4. Operations

Spatial join (points to polygons):

joined = gpd.sjoin(stores, zones[["zone_id", "geometry"]], how="left", predicate="within")
unmatched = joined["zone_id"].isna().sum()
dupes = joined.index.duplicated().sum()   # points on shared borders match twice with "intersects"

Predicate	Use
within / contains	point-in-polygon assignment
intersects	any overlap; boundary points match multiple polygons
dwithin (GeoPandas 1.x with distance=)	within a distance, in the CRS units

Nearest feature with a cap:

utm = stores.estimate_utm_crs()
nearest = gpd.sjoin_nearest(stores.to_crs(utm), stations.to_crs(utm),
                            how="left", max_distance=2000, distance_col="dist_m")

Buffers, distances, areas - always in meters in a projected CRS:

stores_m = stores.to_crs(utm)
stores_m["catchment"] = stores_m.buffer(1000)                    # 1 km
zones_eq = zones.to_crs("EPSG:6933")
zones["area_km2"] = zones_eq.area / 1e6

Overlay and dissolve:

parts = gpd.overlay(parcels, flood_zones, how="intersection", keep_geom_type=True)
by_region = zones.dissolve(by="region", aggfunc={"population": "sum"})

Area-weighted interpolation (population from zones onto a grid) needs equal-area CRS and overlay with the share of each zone’s area; do not assign by centroid when polygons straddle.

Hex aggregation for maps and density: H3 (h3 package, latlng_to_cell(lat, lng, res)) gives equal-ish cells; resolution 7 is ~5 km2, 8 is ~0.7 km2, 9 is ~0.1 km2.

Performance

• GeoPandas joins use the sindex (STRtree) automatically; for custom loops, query gdf.sindex.query(geoms, predicate="intersects") in bulk instead of iterating rows.
• Shapely 2 functions are vectorized over arrays; never apply a Python lambda per row for distance, contains, or buffer.
• Simplify complex polygons before joins if precision allows (simplify(tolerance_m, preserve_topology=True) in a projected CRS).
• Above a few million features, push the work into PostGIS or DuckDB’s spatial extension.
• Prefer GeoParquet (to_parquet) for intermediate data: fast, typed, keeps CRS. Shapefile truncates column names to 10 characters and has a 2 GB limit; use GeoPackage for interchange with desktop GIS.

PostGIS

CREATE INDEX CONCURRENTLY stores_geom_gist ON stores USING gist (geom);

-- Within 1 km, using geography for meters on lon/lat data
SELECT s.id FROM stores s
WHERE ST_DWithin(s.geom::geography, ST_SetSRID(ST_MakePoint($1, $2), 4326)::geography, 1000);

-- Point in polygon join
SELECT z.zone_id, count(*) FROM stores s JOIN zones z ON ST_Within(s.geom, z.geom) GROUP BY 1;

-- K nearest using the index operator
SELECT id, geom <-> ST_SetSRID(ST_MakePoint($1, $2), 4326) AS d
FROM stations ORDER BY d LIMIT 5;

• ST_DWithin uses the index; ST_Distance(...) < 1000 in WHERE does not.
• Casting the column to geography in the query can bypass a geometry index; index the expression (USING gist ((geom::geography))) or store a geography column if most queries are metric.
• Huge polygons (countries, coastlines) make joins slow; ST_Subdivide them into a helper table.
• Repair with ST_MakeValid; check with ST_IsValidReason.

6. Maps

• Interactive check: gdf.explore(column="value", tooltip=["name", "value"]) (folium) - reproject to EPSG:4326 first if needed.
• Static: gdf.plot(column="value", legend=True, scheme="quantiles") (needs mapclassify), add basemap with contextily in EPSG:3857.
• Choropleths map rates or densities (per km2, per capita), not raw counts, or large areas dominate visually.
• Label the CRS, data source, and date on published maps.

Gotchas

• Mixing CRSs in one operation raises a warning in GeoPandas but proceeds in many tools; treat the warning as an error.
• distance between lon/lat geometries returns degrees; one degree of longitude varies from 111 km at the equator to 0 at the poles.
• Antimeridian (180 degrees) crossing polygons break in planar operations; split them.
• Geocoding quality: keep the geocoder’s confidence/match type and exclude low-precision matches (country or city centroid) from fine-grained analysis.
• Location data is often personal data; aggregate before sharing (see data-anonymization).

Verification

• Every layer prints a CRS; operations ran in the planned CRS (assert g.crs == expected).
• Join diagnostics: unmatched count and duplicate matches reported and explained.
• Totals preserved where they should be (population before and after interpolation within 1%; area of dissolved layer equals the sum of parts).
• Spot-check five features on an interactive map against a basemap.