path: root/scripts/update-osm-parks.py

#!/usr/bin/env python3
"""小平市の公園ポリゴンをOpenStreetMapから取得して更新する。

使い方:
    py -X utf8 scripts/update-osm-parks.py

出力:
    src/data/kodaira-parks-osm.json
    src/data/kodaira-boundary.json（小平市の境界線）

OpenStreetMapのデータは ODbL 1.0。地図表示時は
「© OpenStreetMap contributors」の帰属表示が必要。
"""
import difflib
import json
import math
import os
import re
import sys
import time
import unicodedata
import urllib.request

CITY = "小平市"
OURS = "src/data/kodaira-parks.json"
OUT = "src/data/kodaira-parks-osm.json"
BOUNDARY_OUT = "src/data/kodaira-boundary.json"
OVERRIDES = "src/data/park-name-overrides.json"

# OSMの位置補正（メートル）。地図側のLeafletペインに付いていた
# 16pxの余白をCSSで解消したため、現在は0（補正なし）。
# 実位置との差が残る場合だけ正の値を入れる。
NORTH_SHIFT_M = 0
LAT_PER_M = 1 / 111320

# 名前が一致しないポリゴンの自動照合。誤リンクを避けるため保守的に絞る。
NEAR_DIST_M = 20  # 公園の点がポリゴンからこの距離以内なら同一とみなす
NEAR_NAME_MIN = 0.5  # 近接照合で、名前付きポリゴンに求める名前類似度
FUZZY_NAME_MIN = 0.8  # 正規化後の名前類似度がこれを超えたら候補（厳密に大）
FUZZY_MAX_DIST_M = 500  # 名前類似だけで照合する場合の距離上限
# ポリゴン面積と台帳面積の比。極端に違うものは別施設とみなす
AREA_RATIO_MIN = 0.2
AREA_RATIO_MAX = 5.0

# 「緑地（OSM 1234567）」は OSM の ID を名前に含んでいる
OSM_ID_NAME_RE = re.compile(r"^緑地（OSM (\d+)）$")

# 数字・方位・新旧など「位置の違い」を表す文字。これだけが違う名前は
# 別公園の可能性が高い（例: 七小西第1公園 / 七小西第2公園）。
POSITIONAL_CHARS = set("0123456789０１２３４５６７８９一二三四五六七八九十第東西南北上中下前後新")

ENDPOINTS = [
    "https://overpass-api.de/api/interpreter",
    "https://overpass.kumi.systems/api/interpreter",
    "https://overpass.private.coffee/api/interpreter",
]

# leisure=park だけでは取りこぼす児童遊園・広場・緑地も取得する。
# 照合は下の厳しめの条件（完全一致→近接→名前類似）で行う。
QUERY = f"""[out:json][timeout:90];
area["name"="{CITY}"]["boundary"="administrative"]->.a;
(
  way["leisure"~"^(park|playground|garden|recreation_ground|village_green)$"](area.a);
  relation["leisure"~"^(park|playground|garden|recreation_ground|village_green)$"](area.a);
  way["landuse"~"^(recreation_ground|village_green)$"](area.a);
  relation["landuse"~"^(recreation_ground|village_green)$"](area.a);
  way["landuse"="forest"]["name"~"コゲラの森|こげらの森"](area.a);
  relation["landuse"="forest"]["name"~"コゲラの森|こげらの森"](area.a);
);
out geom;
"""

BOUNDARY_QUERY = f"""[out:json][timeout:60];
rel["name"="{CITY}"]["boundary"="administrative"];
out geom;
"""

# 名称が一致しないが、住所・座標から同一と確認できた公園。
ALIAS = {
    "小平中央公園": "中央公園",
    "じゃぶじゃぶ公園": "天神じゃぶじゃぶ公園",
    "リサイクルセンター東公園": "リサイクルセンター南公園",
    "小平市立薬用植物園前公園": "薬用植物園前公園",
    "花南ポケットパーク": "花南ポケットパーク公園",
    "花５四季公園": "花五四季公園",
    "天神町のびのび広場": "天神町のびのび広場公園",
}


def norm(s):
    s = "".join(unicodedata.normalize("NFKC", str(s or "")).split())
    return s.replace("ヶ", "ケ").replace("ツ", "ッ")


def similarity(a, b):
    return difflib.SequenceMatcher(None, a, b).ratio()


def strip_positional(s):
    return "".join(c for c in s if c not in POSITIONAL_CHARS)


def rings(geom):
    if geom["type"] == "Polygon":
        return [geom["coordinates"][0]]
    return [poly[0] for poly in geom["coordinates"]]


def point_in_ring(lon, lat, ring):
    inside = False
    for (x1, y1), (x2, y2) in zip(ring, ring[1:]):
        if (y1 > lat) != (y2 > lat):
            if lon < (x2 - x1) * (lat - y1) / (y2 - y1) + x1:
                inside = not inside
    return inside


def segment_dist_m(lon, lat, a, b):
    """点から線分までの距離（メートル）。"""
    k = math.cos(math.radians(lat))
    x1, y1 = a[0] * k, a[1]
    x2, y2 = b[0] * k, b[1]
    px, py = lon * k, lat
    dx, dy = x2 - x1, y2 - y1
    if dx == 0 and dy == 0:
        t = 0.0
    else:
        t = max(0.0, min(1.0, ((px - x1) * dx + (py - y1) * dy) / (dx * dx + dy * dy)))
    return math.hypot(px - (x1 + t * dx), py - (y1 + t * dy)) * 111320


def polygon_dist_m(geom, lat, lon):
    """公園の代表点からポリゴンまでの距離。内側なら0。"""
    rs = rings(geom)
    for ring in rs:
        if point_in_ring(lon, lat, ring):
            return 0.0
    return min(
        (segment_dist_m(lon, lat, a, b) for ring in rs for a, b in zip(ring, ring[1:])),
        default=float("inf"),
    )


def polygon_area_m2(geom):
    total = 0.0
    for ring in rings(geom):
        if len(ring) < 4:
            continue
        lat0 = sum(p[1] for p in ring) / len(ring)
        kx = 111320 * math.cos(math.radians(lat0))
        ky = 111320
        s = 0.0
        for (x1, y1), (x2, y2) in zip(ring, ring[1:] + ring[:1]):
            s += (x1 * kx) * (y2 * ky) - (x2 * kx) * (y1 * ky)
        total += abs(s) / 2
    return total


def fetch(query):
    last = None
    for ep in ENDPOINTS:
        try:
            req = urllib.request.Request(
                ep,
                data=query.encode("utf-8"),
                headers={"User-Agent": "yasutakeyohei.com park map updater"},
            )
            with urllib.request.urlopen(req, timeout=120) as r:
                data = json.loads(r.read().decode("utf-8"))
            if data.get("elements"):
                return data
            last = "0 elements"
        except Exception as e:  # noqa: BLE001
            last = e
        print(f"  retry: {ep}: {last}", file=sys.stderr)
        time.sleep(2)
    raise SystemExit(f"Overpass API failed: {last}")


def pt(g):
    return (round(g["lon"], 6), round(g["lat"] + NORTH_SHIFT_M * LAT_PER_M, 6))


def ring_from_ways(ways):
    rings = []
    remaining = [list(w) for w in ways if len(w) >= 2]
    while remaining:
        ring = remaining.pop(0)
        while ring[0] != ring[-1] and remaining:
            found = False
            for i, w in enumerate(remaining):
                if w[0] == ring[-1]:
                    ring.extend(w[1:]); remaining.pop(i); found = True; break
                if w[-1] == ring[-1]:
                    ring.extend(list(reversed(w))[1:]); remaining.pop(i); found = True; break
                if w[0] == ring[0]:
                    ring = list(reversed(w))[:-1] + ring; remaining.pop(i); found = True; break
                if w[-1] == ring[0]:
                    ring = w[:-1] + ring; remaining.pop(i); found = True; break
            if not found:
                break
        if ring[0] != ring[-1]:
            ring.append(ring[0])
        if len(ring) >= 4:
            rings.append(ring)
    return rings


def build_boundary():
    """小平市の行政境界をGeoJSONにする。"""
    data = fetch(BOUNDARY_QUERY)
    if not data.get("elements"):
        print("  boundary relation not found", file=sys.stderr)
        return None
    rel = data["elements"][0]
    ways = []
    for m in rel.get("members", []):
        if m.get("role") == "outer" and m.get("geometry"):
            ways.append([pt(p) for p in m["geometry"]])
    rings = ring_from_ways(ways)
    if not rings:
        print("  boundary ring not assembled", file=sys.stderr)
        return None
    if len(rings) == 1:
        geom = {"type": "Polygon", "coordinates": [rings[0]]}
    else:
        geom = {"type": "MultiPolygon", "coordinates": [[r] for r in rings]}
    props = {
        "name": rel.get("tags", {}).get("name", CITY),
        "boundary": "administrative",
        "admin_level": rel.get("tags", {}).get("admin_level", ""),
    }
    fc = {"type": "FeatureCollection", "features": [{"type": "Feature", "properties": props, "geometry": geom}]}
    return fc, len(rings), sum(len(r) for r in rings)


def main():
    ours = json.load(open(OURS, encoding="utf-8"))
    our_names = {norm(r["公園名"]): r["公園名"] for r in ours}
    alias_norm = {norm(k): v for k, v in ALIAS.items()}

    # 市の公園の代表点（緯度・経度）と台帳面積
    our_points = []
    our_areas = {}
    our_id_names = {}
    for r in ours:
        name = r["公園名"]
        m = OSM_ID_NAME_RE.match(name)
        if m:
            our_id_names[int(m.group(1))] = name
        try:
            our_areas[name] = float(r.get("面積㎡") or "")
        except ValueError:
            pass
        try:
            lat = float(r.get("緯度") or "")
            lon = float(r.get("経度") or "")
        except ValueError:
            continue
        our_points.append((name, lat, lon))

    # 前回の出力（新規に照合できたもののレポートと、部分取得時の補完に使う）
    previous_matched = set()
    previous_features = {}
    if os.path.exists(OUT):
        try:
            prev = json.load(open(OUT, encoding="utf-8"))
            for f in prev.get("features", []):
                p = f.get("properties", {})
                key = f"{p.get('osm_type')}/{p.get('osm_id')}"
                previous_features[key] = f
                if p.get("matched"):
                    previous_matched.add(key)
        except (ValueError, KeyError):
            pass

    # OSM ID単位の手動上書き（osm_type/osm_id -> 市の公園名）
    overrides = {}
    if os.path.exists(OVERRIDES):
        raw = json.load(open(OVERRIDES, encoding="utf-8"))
        for key, entry in raw.items():
            city = entry.get("city", "") if isinstance(entry, dict) else str(entry)
            city = str(city).strip()
            if not city:
                continue
            canonical = our_names.get(norm(city))
            if canonical:
                overrides[key] = canonical
            else:
                print(f"  warning: override city not found: {key} -> {city}", file=sys.stderr)

    def exact_match(name):
        k = norm(name)
        if k in our_names:
            return our_names[k], "exact"
        if k in alias_norm:
            return alias_norm[k], "alias"
        # 「〇〇児童公園」→「〇〇公園」のような表記差を吸収
        k2 = k.replace("児童", "")
        if k2 in our_names:
            return our_names[k2], "exact"
        return "", ""

    def area_ok(park, osm_area):
        official = our_areas.get(park)
        if not official:
            return True
        return AREA_RATIO_MIN <= osm_area / official <= AREA_RATIO_MAX

    def fuzzy_match(name, geom):
        """ポリゴンが重なる、または名前が十分似ている市の公園を返す。"""
        dists = [(polygon_dist_m(geom, q[1], q[2]), q[0]) for q in our_points]
        near = sorted(dists)
        near_dist, near_park = near[0]
        osm_area = polygon_area_m2(geom)
        if near_dist <= NEAR_DIST_M:
            # 複数の公園が重なる場合はどちらか判定できないので見送る
            if len(near) > 1 and near[1][0] <= NEAR_DIST_M:
                return "", "", near_dist, 0.0
            # 別名の公園は、ポリゴンが同じ場所にあっても結びつけない
            if name and similarity(norm(name), norm(near_park)) < NEAR_NAME_MIN:
                return "", "", near_dist, 0.0
            # 極端に面積が違う場合は別施設（大きなグラウンド等）とみなす
            if not area_ok(near_park, osm_area):
                return "", "", near_dist, 0.0
            return near_park, "near", near_dist, 0.0

        if not name:
            return "", "", near_dist, 0.0
        k = norm(name)
        ranked = sorted(
            ((similarity(k, norm(park)), dist, park) for dist, park in dists),
            key=lambda x: (-x[0], x[1]),
        )
        sim, dist, park = ranked[0]
        if sim <= FUZZY_NAME_MIN or dist > FUZZY_MAX_DIST_M:
            return "", "", dist, sim
        # 数字・方位だけの違いは別公園の可能性が高い。重なっている場合のみ採用
        if strip_positional(k) == strip_positional(norm(park)) and dist > NEAR_DIST_M:
            return "", "", dist, sim
        if not area_ok(park, osm_area):
            return "", "", dist, sim
        return park, "fuzzy", dist, sim

    data = fetch(QUERY)
    fetched_keys = {f"{e['type']}/{e['id']}" for e in data["elements"]}
    features = []
    applied = 0
    covered = set()
    new_matches = []

    # leisure=park を先に処理し、公園本体が広場・遊び場より先に枠を取るようにする
    elements = sorted(data["elements"], key=lambda e: (e.get("tags", {}).get("leisure") != "park", e["type"], e["id"]))
    for e in elements:
        tags = e.get("tags", {})
        name = tags.get("name", "")
        key = f"{e['type']}/{e['id']}"

        if e["type"] == "way":
            ring = [pt(p) for p in e.get("geometry", [])]
            if len(ring) < 3:
                continue
            if ring[0] != ring[-1]:
                ring.append(ring[0])
            geom = {"type": "Polygon", "coordinates": [ring]}
        else:
            outers, inners = [], []
            for m in e.get("members", []):
                pts = [pt(p) for p in m.get("geometry", [])]
                if len(pts) < 2:
                    continue
                (inners if m.get("role") == "inner" else outers).append(pts)
            outer_rings = ring_from_ways(outers)
            inner_rings = ring_from_ways(inners)
            if not outer_rings:
                continue
            if len(outer_rings) == 1:
                geom = {"type": "Polygon", "coordinates": [outer_rings[0]] + inner_rings}
            else:
                geom = {"type": "MultiPolygon", "coordinates": [[r] + inner_rings for r in outer_rings]}

        our_name, reason = overrides.get(key, ""), "override" if key in overrides else ""
        dist = sim = None
        if key in overrides:
            applied += 1
        elif not our_name and e["id"] in our_id_names:
            # 「緑地（OSM 1234567）」は ID が名前に入っているので直接結びつける
            our_name, reason = our_id_names[e["id"]], "osm_id"
        elif not our_name and not exact_match(name)[0]:
            our_name, reason, dist, sim = fuzzy_match(name, geom)
        if not our_name:
            our_name, reason = exact_match(name)

        # 公園本体（leisure=park）以外が既に整備済みの公園に重なる場合は、
        # 公園内の遊び場・広場とみなして取り込まない。
        core = tags.get("leisure") == "park"
        if our_name in covered and not (core and reason in ("exact", "alias", "override")):
            continue

        props = {
            "name": name,
            "our_name": our_name,
            "matched": bool(our_name),
            "osm_type": e["type"],
            "osm_id": e["id"],
        }
        if our_name:
            props["match"] = reason
        for k in ["leisure", "landuse", "addr:city", "addr:street", "addr:housenumber", "opening_hours"]:
            if k in tags:
                props[k] = tags[k]

        if our_name:
            covered.add(our_name)
            if key not in previous_matched:
                new_matches.append((our_name, key, reason, tags.get("leisure") or tags.get("landuse", ""), dist, sim))

        features.append({"type": "Feature", "properties": props, "geometry": geom})

    # Overpassが部分的な結果を返すと、override対象のFeatureが抜けることがある。
    # その場合は前回のデータを保持して、公園のポリゴンを失わないようにする。
    for key, city in overrides.items():
        if key in fetched_keys:
            continue
        prev_feat = previous_features.get(key)
        if not prev_feat:
            continue
        feat = json.loads(json.dumps(prev_feat))
        feat["properties"]["our_name"] = city
        feat["properties"]["matched"] = True
        feat["properties"]["match"] = "override"
        features.append(feat)
        applied += 1
        print(
            f"  warning: {key} は今回の取得に含まれないため前回のデータを保持しました",
            file=sys.stderr,
        )

    features.sort(key=lambda f: (f["properties"].get("name", ""), f["properties"]["osm_type"], f["properties"]["osm_id"]))
    out = {"type": "FeatureCollection", "features": features}
    json.dump(out, open(OUT, "w", encoding="utf-8"), ensure_ascii=False, separators=(",", ":"))

    matched = sum(1 for f in features if f["properties"]["matched"])
    named = sum(1 for f in features if f["properties"].get("name"))
    print(f"features {len(features)} (named {named}, matched {matched}) -> {OUT}")
    if overrides:
        print(f"  overrides applied {applied} / {len(overrides)}")
    if new_matches:
        print(f"  newly matched: {len(new_matches)}")
        for our_name, key, reason, tag, dist, sim in sorted(new_matches):
            detail = []
            if sim:
                detail.append(f"類似{sim:.2f}")
            if dist is not None:
                detail.append(f"{dist:.0f}m")
            print(f"    {our_name} <- {key} [{tag}/{reason}] {' '.join(detail)}")

    boundary = build_boundary()
    if boundary:
        boundary_data, ring_count, point_count = boundary
        json.dump(boundary_data, open(BOUNDARY_OUT, "w", encoding="utf-8"), ensure_ascii=False, separators=(",", ":"))
        print(f"boundary rings {ring_count}, points {point_count} -> {BOUNDARY_OUT}")


if __name__ == "__main__":
    main()