Turn Geographic Data into Maps and Statistics SKILL.md

Slice and map geospatial data without a GIS license. Use for shapefiles, GeoPackage, GeoJSON, or GeoTIFF rasters; zonal statistics like mean elevation per county; buffers and spatial joins; or publication maps with choropleths, hillshade, insets, and zoom breakouts. Downloads come from vetted official sources such as USGS, NOAA, and the US Census Bureau.
Back to pack
---
name: geospatial-analyst
description: Slice and map geospatial data without a GIS license. Use for shapefiles, GeoPackage, GeoJSON, or GeoTIFF rasters; zonal statistics like mean elevation per county; buffers and spatial joins; or publication maps with choropleths, hillshade, insets, and zoom breakouts. Downloads come from vetted official sources such as USGS, NOAA, and the US Census Bureau.
category: Make sense of your data
type: bundled
requires:
  packages:
    - geopandas>=1.0
    - rasterio>=1.4
    - pyogrio>=0.10
    - shapely>=2.0
    - pyproj>=3.7
    - mapclassify>=2.6
setup:
  - "No connector is required; the skill runs on the local Python geo stack."
  - "First use prompts to install geopandas, rasterio, and friends into the run_code environment; approve once."
  - "Datasets download from vetted official sources (USGS, NOAA, Census, Natural Earth) with a recorded provenance file."
---


# Geospatial Analyst

An ArcGIS-class analysis and cartography workflow built on the open-source
Python geo stack (GeoPandas/Rasterio/Shapely), using redeep's normal tools:
`run_code` for exploration, `bash` for the bundled scripts, `analyze_image` to
inspect rendered maps. Everything is reversible: data, layers, and figures live
in a **map project** directory driven by a `layers.json` manifest.

## The map project (undo/redo backbone)

```
map_project/
  layers.json      <- the manifest: canvas + ordered layers + history
  data/            <- downloaded datasets, each with _provenance.json
  renders/         <- outputs (map.png, stats JSON/CSV)
```

Start every analysis/map by creating one:

```bash
python3 "<skill>/scripts/geokit.py" init map_project --crs EPSG:4326 --title "..."
```

Layers are the unit of undo. Each add/restyle/show/hide/remove is recorded in
`layers.json` history; reverting a change = editing the manifest, re-rendering
= `render_map.py`. **Never** composite everything into one script — add one
layer per geokit command so any element (basemap, choropleth, markers, inset)
can be independently styled, hidden, or deleted without touching others.

Typical z order: hillshade (1) < raster (5) < admin polygons (10) < choropleth
(20) < lines (30) < points (40) < labels (50) < graticule (60) <
scalebar/north_arrow (70) < inset/breakout (80).

## Workflow

1. **Acquire data — vetted sources only.** `fetch_data.py list` shows the
   registry; `fetch_data.py get DATASET --out map_project/data` downloads and
   records provenance (URL, SHA-256, date, license) in `_provenance.json`.
   For elevation use `get usgs_3dep --bbox xmin,ymin,xmax,ymax` (TNM API; pick
   a small bbox — tiles are ~1x1 degree). If the user's needed dataset is not
   in the registry, prefer .gov/.edu/official-CDN hosts and add it to
   `VETTED_SOURCES` in geokit.py after verifying the host, citing the agency.
2. **Explore & slice** in `run_code`: `import geopandas, rasterio`; read with
   `gpd.read_file` / `rasterio.open`; profile columns, CRS, bounds, dtypes.
   Reproject with `to_crs`, clip with `gdf.clip`, buffer with `buffer` — but only
   after projecting to a metric CRS like EPSG:3857 or local UTM (**never buffer
   in degrees**). Save intermediates under `data/` with pyogrio
   (`gdf.to_file(..., engine="pyogrio")`).
3. **Zonal statistics** (raster stats per polygon):
   `zonal_stats.py --raster dem.tif --zones counties.gpkg --stats mean min max
   std count --id-column GEOID --out-json stats.json --out-csv stats.csv`.
   Zones are reprojected to the raster CRS (the raster grid is never
   resampled); zones with no valid pixels report null, not zero. Join stats
   back to zones in `run_code` for choropleths (`zones.merge(stats, on=...)`).
4. **Build the map layer-by-layer** with `geokit.py add`. Style keys: `column`
   (attribute to map), `cmap`, `color`, `alpha`, `scheme`
   (quantile/equal_interval/fisherjenks) + `k`, `categorical: true`,
   `markersize`, `edgecolor`, `linewidth`, `label_column` (annotate features).
   Then `render_map.py map_project` -> `renders/map.png`. Inspect with
   `analyze_image` and iterate: `geokit.py restyle --id choropleth --patch
   '{"style":{"cmap":"magma"}}'` then re-render.
5. **Insets & zoom breakouts**: `geokit.py inset --kind breakout --rect
   '[0.65,0.62,0.32,0.35]' --bbox '[xmin,miny,maxx,maxy]' --layers '[...]'` —
   the renderer highlights the bbox on the main map and draws the detail
   panel; `--kind inset` is a plain secondary map (e.g. a locator country).
6. **Publish-safe outputs**: renders are 300 dpi at journal column widths
   (7.48 in double-column default). Verify legibility with `analyze_image` at
   scale: titles >= 10pt, annotations >= 6pt, distinguishable colors, readable
   legend, scalebar + north arrow for topographic-style figures. The figure
   footer auto-cites each layer's data sources from provenance records.
7. **R option**: for raster-first workflows, R's `sf`/`terra` offer the same
   capabilities (see references/r-spatial.md). Prefer Python unless the user
   asks for R or is already in R.

## Script reference

| Script | Role |
|---|---|
| `geokit.py` | Project manifest CLI: init/add/inset/remove/show/hide/restyle/canvas/list/validate; vetted-source registry; importable from `run_code` (`import geokit`) |
| `fetch_data.py` | Download from vetted official sources with provenance; `list` / `get` |
| `zonal_stats.py` | Per-polygon raster statistics -> JSON + CSV |
| `render_map.py` | Manifest -> publication PNG; layers in z order, legends, insets/breakouts, provenance footer |

All scripts print JSON; consume their output programmatically. Scripts are
black boxes: run them, don't read them. If one fails, its JSON error names the
fix (e.g. missing geo stack -> install via the skill's `requires.packages`).

## Conventions & gotchas

- **CRS discipline**: know every layer's CRS (`gdf.crs`, `src.crs`); reproject
  vectors to the map CRS before styling/plotting; project to metric CRS for
  area/distance work. Mixing CRSes silently misplaces data by kilometers.
- **Degrees are not meters**: buffering/density in EPSG:4326 produces nonsense;
  use EPSG:3857 or local UTM (EPSG:326xx/327xx) for geometry ops.
- **Nodata**: always check `src.nodata`; mask before stats; categorical-code
  rasters need `categorical: true` styling, not continuous colormaps.
- **Attribute joins for choropleths**: match on exact keys (FIPS/GEOID), watch
  dtype mismatches (str vs int), and sanity-check row counts before mapping —
  a silent inner-join loss shows up as holes in the map.
- **Publication style**: subtle graticules, desaturated admin borders,
  perceptually-uniform colormaps (viridis/magma/cividis) for continuous data,
  colorblind-safe categorical palettes, compact legends. When in doubt, render
  then `analyze_image` the PNG.
- **Provenance**: every figure is only as credible as its data lineage — cite
  agency + product + date in the footer (auto via `_provenance.json`); record
  rolling feeds (e.g. USGS earthquake GeoJSON) with their retrieval date.
- **File formats**: prefer GeoPackage (`.gpkg`) for project intermediates
  (single file, no sidecar issues, multiple layers); GeoJSON for interchange;
  keep original shapefile component sets together under `data/<dataset>/`.

## References (read on demand)

- `references/cheatsheet.md` — dense GeoPandas/Rasterio/pyproj recipe card
- `references/data-sources.md` — the vetted-source registry in depth: catalog
  URLs, API patterns, how to vet a new source
- `references/r-spatial.md` — the same workflow in R (sf/terra/ggplot2)
- `references/cartography.md` — publication figure checklist, colormap/scheme
  choice, inset/breakout design patterns, ArcGIS-equivalent recipes