"""build_geometry.py  (AUTHOR SIDE ONLY, students never run this)

Converts the Census cartographic boundary shapefile into a small, simplified
GeoJSON that the student-side mapping code can draw with matplotlib alone.

This is the split that keeps the promise made in the front matter. Reading a
shapefile needs a GIS stack. Drawing a polygon does not. So the authors pay
the GIS dependency once, here, and ship the result.

    python build_geometry.py

Writes ../snapshots/ms_tracts.geojson and ms_counties.geojson.
"""
from __future__ import annotations

import hashlib
import io
import json
import zipfile
from pathlib import Path

import geopandas as gpd
import requests

SNAP = Path(__file__).resolve().parent.parent / "snapshots"
CB = "https://www2.census.gov/geo/tiger/GENZ2023/shp/cb_2023_{fips}_{lvl}_500k.zip"
STATE = "28"

# Albers Equal Area Conic, USGS. Any CONUS-scale map is drawn in this and
# never in raw longitude and latitude, which distorts area at this extent.
EPSG_ALBERS = 5070
TOLERANCE_M = 250   # simplification tolerance, meters


def fetch(level: str) -> gpd.GeoDataFrame:
    # Tracts publish one file per state; counties publish one national file.
    url = CB.format(fips="us" if level == "county" else STATE, lvl=level)
    print(f"  downloading {url.rsplit('/', 1)[-1]} ...", flush=True)
    r = requests.get(url, timeout=600)
    r.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(r.content)) as z:
        name = next(n for n in z.namelist() if n.endswith(".shp"))
        z.extractall(SNAP / "_shp")
    return gpd.read_file(SNAP / "_shp" / name)


def export(level: str, out_name: str, id_field: str) -> None:
    gdf = fetch(level)
    if level == "county":
        gdf = gdf[gdf["STATEFP"] == STATE]
    gdf = gdf.to_crs(EPSG_ALBERS)
    gdf["geometry"] = gdf.geometry.simplify(TOLERANCE_M, preserve_topology=True)
    gdf = gdf[[id_field, "geometry"]].rename(columns={id_field: "id"})

    # Hand-roll the GeoJSON so coordinates round to whole meters. Albers is
    # metric, so rounding costs 1 meter of precision and roughly halves the
    # file, which matters for something shipped in a repository.
    feats = []
    for _, row in gdf.iterrows():
        geom = row.geometry
        polys = ([geom] if geom.geom_type == "Polygon"
                 else list(geom.geoms))
        rings = [[[[round(x), round(y)] for x, y in p.exterior.coords]]
                 for p in polys]
        feats.append({"type": "Feature",
                      "properties": {"id": row["id"]},
                      "geometry": {"type": "MultiPolygon",
                                   "coordinates": rings}})
    obj = {"type": "FeatureCollection", "crs": f"EPSG:{EPSG_ALBERS}",
           "features": feats}
    path = SNAP / out_name
    path.write_text(json.dumps(obj, separators=(",", ":")))
    digest = hashlib.sha256(path.read_bytes()).hexdigest()
    print(f"{out_name:24s} features={len(feats):5d} "
          f"size={path.stat().st_size/1e6:5.2f} MB  sha256={digest[:16]}")


if __name__ == "__main__":
    SNAP.mkdir(parents=True, exist_ok=True)
    export("tract", "ms_tracts.geojson", "GEOID")
    export("county", "ms_counties.geojson", "GEOID")
    print("\nStudents read these with json and draw them with matplotlib.")
