Europa/Nuvolari/tool/generate_places.py

"""Generate the searchable list of Piedmont municipalities.

The app is region-scoped, so a bundled list beats a geocoding service on every
axis that matters here: the search is instant, works with no network, needs no
API key, has no rate limit to respect, and cannot return a result in a place the
app has no radar for.

Source: Istat "Confini delle unità amministrative a fini statistici",
generalised boundaries, published under CC BY 4.0. The shapefile is in WGS84 /
UTM zone 32N, so each municipality's representative point is converted back to
latitude and longitude here rather than at run time.

Dependencies: `pyshp` only. The inverse transverse Mercator is implemented below
instead of pulling in pyproj, because that is the single piece of geodesy this
script needs and a whole projection library is a heavy thing to make everyone
install for it.

Usage:
    python -m pip install pyshp
    python tool/generate_places.py
"""

from __future__ import annotations

import io
import json
import math
import pathlib
import sys
import urllib.request
import zipfile

import shapefile  # pyshp

SOURCE_URL = (
    "https://www.istat.it/storage/cartografia/confini_amministrativi/"
    "generalizzati/2025/Limiti01012025_g.zip"
)
SHAPE_BASE = "Com01012025_g/Com01012025_g_WGS84"

# Istat region code for Piemonte.
REGION_CODE = 1

# Istat province codes to the two-letter abbreviations people actually recognise.
PROVINCES = {
    1: "TO",
    2: "VC",
    3: "NO",
    4: "CN",
    5: "AT",
    6: "AL",
    96: "BI",
    103: "VB",
}

OUT_PATH = (
    pathlib.Path(__file__).resolve().parent.parent
    / "app"
    / "assets"
    / "regions"
    / "piemonte_comuni.json"
)

CACHE_PATH = pathlib.Path(__file__).resolve().parent / ".cache" / "istat_limiti.zip"

# --- WGS84 / UTM 32N inverse -------------------------------------------------

_A = 6378137.0  # WGS84 semi-major axis
_F = 1 / 298.257223563  # flattening
_K0 = 0.9996  # UTM scale factor
_FALSE_EASTING = 500000.0
_CENTRAL_MERIDIAN = math.radians(9.0)  # zone 32


def utm32n_to_wgs84(easting: float, northing: float) -> tuple[float, float]:
    """Converts UTM zone 32N metres to (longitude, latitude) in degrees.

    Snyder's inverse transverse Mercator series, accurate to a few centimetres
    inside the zone — orders of magnitude finer than a municipality centroid
    needs to be.
    """
    e2 = _F * (2 - _F)
    e_prime2 = e2 / (1 - e2)

    x = easting - _FALSE_EASTING
    m = northing / _K0

    mu = m / (_A * (1 - e2 / 4 - 3 * e2**2 / 64 - 5 * e2**3 / 256))
    e1 = (1 - math.sqrt(1 - e2)) / (1 + math.sqrt(1 - e2))

    phi1 = (
        mu
        + (3 * e1 / 2 - 27 * e1**3 / 32) * math.sin(2 * mu)
        + (21 * e1**2 / 16 - 55 * e1**4 / 32) * math.sin(4 * mu)
        + (151 * e1**3 / 96) * math.sin(6 * mu)
        + (1097 * e1**4 / 512) * math.sin(8 * mu)
    )

    sin_phi1 = math.sin(phi1)
    cos_phi1 = math.cos(phi1)
    tan_phi1 = math.tan(phi1)

    c1 = e_prime2 * cos_phi1**2
    t1 = tan_phi1**2
    n1 = _A / math.sqrt(1 - e2 * sin_phi1**2)
    r1 = _A * (1 - e2) / (1 - e2 * sin_phi1**2) ** 1.5
    d = x / (n1 * _K0)

    latitude = phi1 - (n1 * tan_phi1 / r1) * (
        d**2 / 2
        - (5 + 3 * t1 + 10 * c1 - 4 * c1**2 - 9 * e_prime2) * d**4 / 24
        + (61 + 90 * t1 + 298 * c1 + 45 * t1**2 - 252 * e_prime2 - 3 * c1**2)
        * d**6
        / 720
    )

    longitude = _CENTRAL_MERIDIAN + (
        d
        - (1 + 2 * t1 + c1) * d**3 / 6
        + (5 - 2 * c1 + 28 * t1 - 3 * c1**2 + 8 * e_prime2 + 24 * t1**2)
        * d**5
        / 120
    ) / cos_phi1

    return math.degrees(longitude), math.degrees(latitude)


# --- geometry ----------------------------------------------------------------


def ring_area_and_centroid(ring: list[tuple[float, float]]) -> tuple[float, float, float]:
    """Signed area and centroid of a closed ring, by the shoelace formula."""
    area2 = 0.0
    cx = 0.0
    cy = 0.0
    for i in range(len(ring) - 1):
        x0, y0 = ring[i]
        x1, y1 = ring[i + 1]
        cross = x0 * y1 - x1 * y0
        area2 += cross
        cx += (x0 + x1) * cross
        cy += (y0 + y1) * cross
    if area2 == 0:
        # Degenerate ring: fall back to the mean of its vertices.
        xs = [p[0] for p in ring]
        ys = [p[1] for p in ring]
        return 0.0, sum(xs) / len(xs), sum(ys) / len(ys)
    return area2 / 2, cx / (3 * area2), cy / (3 * area2)


def representative_point(shape: shapefile.Shape) -> tuple[float, float]:
    """A point inside the municipality, in the shapefile's own coordinates.

    Uses the centroid of the **largest** ring rather than of the whole shape.
    Many Piedmont municipalities are multipart — a main body plus exclaves — and
    averaging those together can land the point on a neighbour.
    """
    parts = list(shape.parts) + [len(shape.points)]
    best_area = -1.0
    best = (shape.points[0][0], shape.points[0][1])

    for start, end in zip(parts, parts[1:]):
        ring = [(p[0], p[1]) for p in shape.points[start:end]]
        if len(ring) < 4:
            continue
        if ring[0] != ring[-1]:
            ring.append(ring[0])
        area, cx, cy = ring_area_and_centroid(ring)
        if abs(area) > best_area:
            best_area = abs(area)
            best = (cx, cy)

    return best


# --- main --------------------------------------------------------------------


def fetch_archive() -> bytes:
    if CACHE_PATH.exists():
        print(f"using cached {CACHE_PATH}")
        return CACHE_PATH.read_bytes()

    print(f"downloading {SOURCE_URL}")
    with urllib.request.urlopen(SOURCE_URL, timeout=300) as response:
        data = response.read()
    CACHE_PATH.parent.mkdir(parents=True, exist_ok=True)
    CACHE_PATH.write_bytes(data)
    print(f"cached {len(data) / 1_000_000:.1f} MB at {CACHE_PATH}")
    return data


def main() -> int:
    archive = zipfile.ZipFile(io.BytesIO(fetch_archive()))
    reader = shapefile.Reader(
        shp=io.BytesIO(archive.read(SHAPE_BASE + ".shp")),
        dbf=io.BytesIO(archive.read(SHAPE_BASE + ".dbf")),
        shx=io.BytesIO(archive.read(SHAPE_BASE + ".shx")),
        # UTF-8, confirmed by decoding a known name down to its codepoints:
        # "Aglie" ends in a single U+00E8 under utf-8, and in U+00C3 U+00A8
        # under latin-1. A terminal renders both identically, so do not
        # settle this by looking at printed output -- the guard below is
        # what actually holds the line.
        encoding="utf-8",
    )

    places = []
    unknown_provinces = set()

    for record, shape in zip(reader.iterRecords(), reader.iterShapes()):
        if record["COD_REG"] != REGION_CODE:
            continue

        province = PROVINCES.get(record["COD_PROV"])
        if province is None:
            unknown_provinces.add(record["COD_PROV"])
            continue

        easting, northing = representative_point(shape)
        longitude, latitude = utm32n_to_wgs84(easting, northing)

        places.append(
            {
                "name": record["COMUNE"],
                "province": province,
                "istat": record["PRO_COM_T"],
                "lat": round(latitude, 5),
                "lng": round(longitude, 5),
            }
        )

    # Latin-1 is an assumption about someone else's file. If Istat ever
    # re-encodes it, the accents turn into these tell-tale sequences rather
    # than into an error, so check for them explicitly. Written as escapes
    # so this guard cannot itself be broken by an encoding accident.
    MOJIBAKE_MARKERS = ("Ã", "Â", "�")
    mis_decoded = [
        place["name"]
        for place in places
        if any(marker in place["name"] for marker in MOJIBAKE_MARKERS)
    ]
    if mis_decoded:
        print(f"ERROR: names look mis-decoded: {mis_decoded[:5]}", file=sys.stderr)
        return 1

    if unknown_provinces:
        print(f"ERROR: unmapped province codes {sorted(unknown_provinces)}", file=sys.stderr)
        return 1

    places.sort(key=lambda place: place["name"])

    OUT_PATH.parent.mkdir(parents=True, exist_ok=True)
    OUT_PATH.write_text(
        json.dumps(
            {
                "region": "piemonte",
                "source": "Istat — Confini delle unità amministrative a fini statistici, 1 gennaio 2025",
                "license": "CC BY 4.0",
                "sourceUrl": SOURCE_URL,
                "places": places,
            },
            ensure_ascii=False,
            separators=(",", ":"),
        )
        + "\n",
        encoding="utf-8",
    )

    size_kb = OUT_PATH.stat().st_size / 1024
    print(f"{len(places)} comuni written to {OUT_PATH} ({size_kb:.0f} KB)")

    by_province: dict[str, int] = {}
    for place in places:
        by_province[place["province"]] = by_province.get(place["province"], 0) + 1
    print("by province:", dict(sorted(by_province.items())))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())