"""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())