hexcount

Aggregate geographic points into a hexagonal-grid choropleth. Feed it a GeoJSON of points; it bins them into a pointy-top hex grid, counts how many land in each cell, and writes a GeoJSON choropleth plus an optional preview PNG. The GeoJSON output is then styled and exported as a finished map in QGIS.

San Francisco street-tree density: ~195,000 trees from DataSF binned into a 250 m hex grid, styled in QGIS over an OpenStreetMap basemap.

The core (hexcount.core, hexcount.geojson) is pure standard library — no GeoPandas, no Shapely, no NumPy. Plotting and the CLI are opt-in extras, so the engine that does the geometry has zero runtime dependencies and a fully hermetic test suite.

Install

python -m venv .venv && source .venv/bin/activate
pip install -e ".[dev]"      # everything, including tooling
# or, leaner installs:
pip install -e ".[cli]"      # CLI only
pip install -e ".[viz]"      # plotting only

Use

# Bin points into a hex grid (size is center-to-corner, in input units).
hexcount aggregate examples/points.geojson --size 2 --out grid.geojson

# Optional quick preview PNG (needs the [viz] extra).
hexcount plot grid.geojson --out grid.png --cmap viridis

grid.geojson opens directly in QGIS, where the count attribute drives a graduated (choropleth) style.

Reproducing the San Francisco map

The hero image above is built from the live DataSF "Street Tree List":

# 1. Download the dataset as GeoJSON (~190k trees).
curl -L -o sf_trees.geojson \
  'https://data.sfgov.org/resource/tkzw-k3nq.geojson?$limit=300000'

# 2. Drop bad geocodes and project WGS84 -> Web Mercator (EPSG:3857).
python scripts/prep_trees.py sf_trees.geojson sf_trees_3857.geojson

# 3. Bin into 250 m hexes.
hexcount aggregate sf_trees_3857.geojson --size 250 --out sf_grid.geojson

Then load sf_grid.geojson in QGIS, set its CRS to EPSG:3857, and apply a graduated style on count. The points are reprojected to Web Mercator before binning so the hexagons are equal-aspect and align with web basemap tiles — binning raw lon/lat degrees would stretch the cells, since a degree of longitude is shorter than a degree of latitude at San Francisco's latitude.

How it works

A point is mapped to a hex by converting its planar (x, y) to fractional axial coordinates (q, r), then cube-rounding to the nearest valid cell. Cube rounding rounds each of the three cube coordinates independently, then repairs the one with the largest rounding error so the x + y + z == 0 invariant holds — this is what guarantees every point lands in exactly one hex with no gaps or overlaps. Cells are emitted sorted by axial address so output files and tests are byte-stable across runs.

Design notes

• Pure-stdlib core. All geometry and aggregation live behind a dependency-free API, so the test suite is hermetic and the package installs anywhere. Heavy libraries (matplotlib) sit behind the [viz] extra and import at call time, never at module load.
• Ground-truth tests. Fixtures are small enough to hand-count, so the asserted hex counts are real expected values rather than recorded snapshots.
• GIS-standard output. Output is ordinary GeoJSON polygons with a numeric count property, so the result is consumable by any GIS tool, not just this one — the QGIS step is a clean handoff, not a lock-in.

Project layout

src/hexcount/
  core.py      # hex geometry + aggregation (stdlib only)
  geojson.py   # GeoJSON read/write (stdlib only)
  viz.py       # matplotlib preview (viz extra)
  cli.py       # Typer CLI (cli extra)
scripts/
  prep_trees.py  # DataSF tree GeoJSON -> Web Mercator points
tests/           # hermetic unit tests
examples/        # synthetic points dataset

License

MIT