#!/usr/bin/env python3
"""Builds the app's map of Japan's prefectures from Natural Earth.

Natural Earth's vector data is in the public domain. This script downloads
its Admin 1 (states, provinces) layer at 1:10m, release v5.1.2, checks the
download against the hash recorded below, keeps Japan's 47 prefectures, and
reduces them to a compact outline for display: Douglas-Peucker simplification
(0.006 degrees), islands too small to see dropped (each prefecture keeps its
largest), coordinates rounded to 4 decimals. The result is an approximation,
not survey-accurate data.

Each region is named by its ISO 3166-2 code (JP-01 ... JP-47), with Natural
Earth's English and Japanese names.

Run from the repository root with Python 3 and no other packages:

    python3 tool/maps/japan_prefectures.py

It rewrites halkyon_learn_flutter/assets/maps/japan_prefectures.geojson.
"""

import hashlib
import json
import math
import sys
import urllib.request
from pathlib import Path

SOURCE = (
    "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/"
    "v5.1.2/geojson/ne_10m_admin_1_states_provinces.geojson"
)
SOURCE_SHA256 = "22d0e3ad85eb3e27f17cabf8ba2d50e554fbc27a87796ff891d958185da62fb5"
TARGET = Path("halkyon_learn_flutter/assets/maps/japan_prefectures.geojson")

TOLERANCE = 0.006  # degrees
MIN_AREA = 0.0005  # square degrees; smaller islands are dropped
DECIMALS = 4


def simplify(points, tolerance):
    """Douglas-Peucker on an open polyline of (x, y) pairs."""
    if len(points) < 3:
        return points
    first, last = points[0], points[-1]
    dx, dy = last[0] - first[0], last[1] - first[1]
    length = math.hypot(dx, dy)
    farthest, index = -1.0, 0
    for i in range(1, len(points) - 1):
        px, py = points[i]
        if length == 0:
            distance = math.hypot(px - first[0], py - first[1])
        else:
            distance = abs(dy * px - dx * py + last[0] * first[1]
                           - last[1] * first[0]) / length
        if distance > farthest:
            farthest, index = distance, i
    if farthest <= tolerance:
        return [first, last]
    left = simplify(points[: index + 1], tolerance)
    right = simplify(points[index:], tolerance)
    return left[:-1] + right


def simplify_ring(ring):
    """A closed ring simplified, rounded, without repeated points; None when
    it collapses below a triangle."""
    # Split at the point farthest from the start, so a closed ring (first
    # point = last point) simplifies as two open halves.
    start = ring[0]
    far = max(range(len(ring)), key=lambda i: math.hypot(
        ring[i][0] - start[0], ring[i][1] - start[1]))
    halves = simplify(ring[: far + 1], TOLERANCE)[:-1] + simplify(
        ring[far:], TOLERANCE)
    rounded = []
    for x, y in halves:
        point = [round(x, DECIMALS), round(y, DECIMALS)]
        if not rounded or rounded[-1] != point:
            rounded.append(point)
    if rounded[0] != rounded[-1]:
        rounded.append(rounded[0])
    return rounded if len(rounded) >= 4 else None


def area(ring):
    """The absolute area of a ring, in square degrees (shoelace)."""
    return abs(sum(x1 * y2 - x2 * y1
                   for (x1, y1), (x2, y2) in zip(ring, ring[1:]))) / 2


def polygons_of(geometry):
    if geometry["type"] == "Polygon":
        return [geometry["coordinates"]]
    if geometry["type"] == "MultiPolygon":
        return geometry["coordinates"]
    raise ValueError(f"Unexpected geometry {geometry['type']}")


def reduce(geometry):
    polygons = polygons_of(geometry)
    largest = max(range(len(polygons)), key=lambda i: area(polygons[i][0]))
    kept = []
    for i, polygon in enumerate(polygons):
        if i != largest and area(polygon[0]) < MIN_AREA:
            continue
        exterior = simplify_ring(polygon[0])
        if exterior is None:
            continue
        holes = [hole for hole in map(simplify_ring, polygon[1:])
                 if hole is not None]
        kept.append([exterior, *holes])
    return {"type": "MultiPolygon", "coordinates": kept}


def main():
    print(f"Downloading {SOURCE}", file=sys.stderr)
    data = urllib.request.urlopen(SOURCE).read()
    digest = hashlib.sha256(data).hexdigest()
    if SOURCE_SHA256 is not None and digest != SOURCE_SHA256:
        sys.exit(f"The download's SHA-256 is {digest}, not {SOURCE_SHA256}.")
    print(f"SHA-256 {digest}", file=sys.stderr)
    source = json.loads(data)
    features = []
    for feature in source["features"]:
        properties = feature["properties"]
        if properties.get("adm0_a3") != "JPN":
            continue
        features.append({
            "type": "Feature",
            "properties": {
                "code": properties["iso_3166_2"],
                "name": properties["name"],
                "name_ja": properties["name_ja"],
            },
            "geometry": reduce(feature["geometry"]),
        })
    features.sort(key=lambda feature: feature["properties"]["code"])
    codes = [feature["properties"]["code"] for feature in features]
    if codes != [f"JP-{n:02d}" for n in range(1, 48)]:
        sys.exit(f"Expected the 47 prefectures JP-01 ... JP-47, got {codes}")
    TARGET.write_text(json.dumps(
        {"type": "FeatureCollection", "features": features},
        ensure_ascii=False, separators=(",", ":")) + "\n", encoding="utf-8")
    print(f"Wrote {TARGET} ({TARGET.stat().st_size} bytes)", file=sys.stderr)


if __name__ == "__main__":
    main()