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

#!/usr/bin/env python3
"""小平市の境界から500m以内にある近隣市の公園をOpenStreetMapから抽出する。

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

入力:
    src/data/kodaira-boundary.json   小平市の境界（EPSG:4326）
    src/data/kodaira-parks-osm.json  小平市内の公園（重複除外に使う）
出力:
    src/data/kodaira-nearby-parks.json

leisure=park の way / relation を取得し、小平市内のものと既存データに
含まれるものを除いたうえで、重心または頂点が境界から500m以内の
名前付き公園だけを距離順に書き出す。

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

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
DATA_DIR = os.path.join(ROOT, "src", "data")
BOUNDARY_IN = os.path.join(DATA_DIR, "kodaira-boundary.json")
EXISTING_IN = os.path.join(DATA_DIR, "kodaira-parks-osm.json")
OUT = os.path.join(DATA_DIR, "kodaira-nearby-parks.json")
NEARBY_MANUAL = os.path.join(ROOT, "reference", "qgis", "nearby-parks-manual.geojson")

BBOX = (35.68, 139.39, 35.76, 139.56)  # south, west, north, east
MAX_DIST_M = 500.0
M_PER_DEG_LAT = 111320.0
NOTE = "市外・小平市境界から500m以内"
PARK_NAME_RE = re.compile(r"(公園|緑地|広場|遊園|庭園|森)")

# OSMに公園ポリゴンがない近隣公園。位置・面積は自治体の公開情報による。
# 地図では点（範囲未取得）として表示する。
MANUAL_FEATURES = [
    {
        "name": "北町パンダ公園",
        "distance_m": 94,
        "area_m2": 1675,
        "osm_type": "manual",
        "osm_id": "hokucho-panda",
        "city": "国分寺市",
        "lat": 35.721831,
        "lon": 139.447225,
        "note": "市外・小平市境界から500m以内（公園の範囲は未取得。位置は国分寺市の公園情報による）",
    },
]

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


def fetch_overpass():
    """leisure=park を取得する。全エンドポイントが失敗したら例外を投げる。"""
    south, west, north, east = BBOX
    query = (
        "[out:json][timeout:180];"
        "("
        f'way["leisure"="park"]({south},{west},{north},{east});'
        f'relation["leisure"="park"]({south},{west},{north},{east});'
        f'node["name"~"公園|緑地|広場|遊園|庭園|森"]({south},{west},{north},{east});'
        ");"
        "out geom;"
    )
    last_err = None
    for endpoint in ENDPOINTS:
        for attempt in range(2):
            try:
                req = urllib.request.Request(
                    endpoint,
                    data=query.encode("utf-8"),
                    headers={"User-Agent": "kodaira-nearby-parks/1.0 (yasutakeyohei.com)"},
                )
                with urllib.request.urlopen(req, timeout=240) as res:
                    data = json.load(res)
                if "elements" not in data:
                    raise ValueError("レスポンスに elements がありません")
                print(f"  {endpoint} から {len(data['elements'])} 件取得")
                return data["elements"]
            except Exception as exc:  # noqa: BLE001 - 全滅時のみ異常終了する
                last_err = exc
                print(f"  {endpoint} 失敗 ({attempt + 1}/2): {exc}", file=sys.stderr)
                time.sleep(5 * (attempt + 1))
    raise SystemExit(f"Overpass APIの取得に失敗しました（出力は作成していません）: {last_err}")


def projector(lat0, lon0):
    """緯度 lat0・経度 lon0 を原点にした局所平面（メートル）への投影。

    原点を引かないと座標が1e7規模になり、小さい公園の重心計算で
    桁落ちして十数メートルの誤差が出る。
    """
    kx = M_PER_DEG_LAT * math.cos(math.radians(lat0))

    def to_xy(lon, lat):
        return ((lon - lon0) * kx, (lat - lat0) * M_PER_DEG_LAT)

    return to_xy


def same_pt(a, b, eps=1e-9):
    return abs(a[0] - b[0]) <= eps and abs(a[1] - b[1]) <= eps


def is_closed(ring, eps=1e-9):
    return len(ring) >= 4 and same_pt(ring[0], ring[-1], eps)


def stitch_rings(seqs):
    """端点が一致するwayの座標列をつないで閉じたリングにする。"""
    pending = [list(s) for s in seqs if len(s) >= 2]
    rings = []
    while pending:
        cur = pending.pop(0)
        growing = True
        while not is_closed(cur) and growing:
            growing = False
            for i, seq in enumerate(pending):
                if same_pt(cur[-1], seq[0]):
                    cur = cur + seq[1:]
                elif same_pt(cur[-1], seq[-1]):
                    cur = cur + seq[-2::-1]
                elif same_pt(cur[0], seq[-1]):
                    cur = seq[:-1] + cur
                elif same_pt(cur[0], seq[0]):
                    cur = seq[:0:-1] + cur
                else:
                    continue
                del pending[i]
                growing = True
                break
        if len(cur) >= 3:
            if not is_closed(cur):
                cur = cur + [cur[0]]
            rings.append(cur)
    return rings


def point_in_ring(x, y, ring):
    """レイキャスティングによる点の内外判定（平面座標でも経緯度でも可）。"""
    inside = False
    j = len(ring) - 1
    for i in range(len(ring)):
        xi, yi = ring[i]
        xj, yj = ring[j]
        if (yi > y) != (yj > y):
            if x < (xj - xi) * (y - yi) / (yj - yi) + xi:
                inside = not inside
        j = i
    return inside


def polygon_contains(polygon, x, y):
    if not point_in_ring(x, y, polygon[0]):
        return False
    return not any(point_in_ring(x, y, hole) for hole in polygon[1:])


def name_from_named_node(polygons_ll, named_nodes, used_names):
    """名前のない公園ポリゴンの中にある名前付きノードから公園名を推測する。

    公園ポリゴンの名称がノード側にしかない場合（例: 窪東公園）に使う。
    既に別の公園で使われている名前は候補から除く。
    複数ある場合は重心に最も近いものを選ぶ。
    """
    best_name = ""
    best_dist = float("inf")
    for polygon in polygons_ll:
        ring = polygon[0]
        xs = [pt[0] for pt in ring]
        ys = [pt[1] for pt in ring]
        cx = sum(xs) / len(xs)
        cy = sum(ys) / len(ys)
        for name, lon, lat in named_nodes:
            if name in used_names:
                continue
            if not polygon_contains(polygon, lon, lat):
                continue
            dist = (lon - cx) ** 2 + (lat - cy) ** 2
            if dist < best_dist:
                best_dist = dist
                best_name = name
    return best_name


def way_polygons(el):
    pts = [(p["lon"], p["lat"]) for p in el.get("geometry") or [] if p and "lon" in p]
    if len(pts) < 3:
        return []
    if not same_pt(pts[0], pts[-1]):
        pts.append(pts[0])
    return [[pts]]


def relation_polygons(el):
    outer, inner = [], []
    for member in el.get("members") or []:
        if member.get("type") != "way":
            continue
        pts = [(p["lon"], p["lat"]) for p in member.get("geometry") or [] if p and "lon" in p]
        if len(pts) >= 2:
            (inner if member.get("role") == "inner" else outer).append(pts)
    inner_rings = stitch_rings(inner)
    polygons = []
    for ring in stitch_rings(outer):
        holes = [r for r in inner_rings if point_in_ring(r[0][0], r[0][1], ring)]
        polygons.append([ring] + holes)
    return polygons


def build_polygons(el):
    if el.get("type") == "relation":
        return relation_polygons(el)
    return way_polygons(el)


def ring_area(ring):
    total = 0.0
    for (x1, y1), (x2, y2) in zip(ring, ring[1:] + ring[:1]):
        total += x1 * y2 - x2 * y1
    return total / 2.0


def ring_centroid(ring):
    area = ring_area(ring)
    if abs(area) < 1e-9:
        xs = [p[0] for p in ring]
        ys = [p[1] for p in ring]
        return (sum(xs) / len(xs), sum(ys) / len(ys)), 0.0
    cx = cy = 0.0
    for (x1, y1), (x2, y2) in zip(ring, ring[1:] + ring[:1]):
        cross = x1 * y2 - x2 * y1
        cx += (x1 + x2) * cross
        cy += (y1 + y2) * cross
    return (cx / (6.0 * area), cy / (6.0 * area)), abs(area)


def polygon_centroid(polygon):
    """外周と穴を考慮した多角形の重心と面積。"""
    (ox, oy), outer_area = ring_centroid(polygon[0])
    sx, sy, net = ox * outer_area, oy * outer_area, outer_area
    for ring in polygon[1:]:
        (cx, cy), area = ring_centroid(ring)
        sx -= cx * area
        sy -= cy * area
        net -= area
    if net <= 1e-9:
        return (ox, oy), outer_area
    return (sx / net, sy / net), net


class BoundaryIndex:
    """小平市境界の内外判定と、任意の点からの最短距離（メートル）。"""

    def __init__(self, polygons):
        self.polygons = polygons
        self.segments = []
        for polygon in polygons:
            for ring in polygon:
                for (x1, y1), (x2, y2) in zip(ring, ring[1:] + ring[:1]):
                    dx, dy = x2 - x1, y2 - y1
                    self.segments.append((x1, y1, dx, dy, dx * dx + dy * dy))

    def contains(self, x, y):
        for polygon in self.polygons:
            if point_in_ring(x, y, polygon[0]) and not any(
                point_in_ring(x, y, hole) for hole in polygon[1:]
            ):
                return True
        return False

    def min_distance(self, points):
        """points（メートル座標の配列）から境界線までの最短距離。"""
        best = float("inf")
        for px, py in points:
            for ax, ay, dx, dy, len2 in self.segments:
                if len2 > 0.0:
                    t = ((px - ax) * dx + (py - ay) * dy) / len2
                    if t < 0.0:
                        t = 0.0
                    elif t > 1.0:
                        t = 1.0
                else:
                    t = 0.0
                dist = math.hypot(px - ax - t * dx, py - ay - t * dy)
                if dist < best:
                    best = dist
                    if best == 0.0:
                        return 0.0
        return best


def load_boundary():
    data = json.load(open(BOUNDARY_IN, encoding="utf-8"))
    polygons = []
    for feature in data["features"]:
        geom = feature["geometry"]
        if geom["type"] == "Polygon":
            polygons.append(geom["coordinates"])
        elif geom["type"] == "MultiPolygon":
            polygons.extend(geom["coordinates"])
    if not polygons:
        raise SystemExit("小平市の境界ポリゴンが見つかりません")
    return polygons


def existing_keys():
    data = json.load(open(EXISTING_IN, encoding="utf-8"))
    return {
        (f["properties"].get("osm_type"), f["properties"].get("osm_id"))
        for f in data["features"]
    }


def exclusion_polygons():
    """近隣公園から除外する大きな公園のポリゴン（経緯度）。

    都立小金井公園の内側にある広場・公園は、独立した近隣公園として
    重複表示しない。小金井公園本体は小平市の一覧側に残す。
    """
    data = json.load(open(EXISTING_IN, encoding="utf-8"))
    polygons = []
    for feature in data["features"]:
        name = feature["properties"].get("our_name") or ""
        if name != "都立小金井公園":
            continue
        geom = feature["geometry"]
        if geom["type"] == "Polygon":
            polygons.append(geom["coordinates"])
        elif geom["type"] == "MultiPolygon":
            polygons.extend(geom["coordinates"])
    return polygons


def manual_polygons_ll(geom):
    if geom.get("type") == "Polygon":
        return [geom.get("coordinates") or []]
    if geom.get("type") == "MultiPolygon":
        return geom.get("coordinates") or []
    return []


def merge_nearby_manual(features, boundary, to_xy):
    """reference/qgis/nearby-parks-manual.geojson の手入力ポリゴンを反映する。

    `our_name` が既存の近隣公園と一致すれば、点やOSMポリゴンを
    手入力ポリゴンで置き換える。一致しなければ新しい公園として追加する。
    """
    if not os.path.exists(NEARBY_MANUAL):
        return 0
    try:
        data = json.load(open(NEARBY_MANUAL, encoding="utf-8"))
    except (OSError, ValueError) as exc:
        print(f"  warning: {NEARBY_MANUAL} を読めません: {exc}", file=sys.stderr)
        return 0

    by_name = {f["properties"]["name"]: f for f in features}
    count = 0
    for index, feat in enumerate(data.get("features", []), 1):
        props = feat.get("properties") or {}
        name = (props.get("our_name") or props.get("name") or "").strip()
        geom = feat.get("geometry") or {}
        polygons_ll = manual_polygons_ll(geom)
        if not name or not polygons_ll:
            continue

        polygons = []
        for polygon in polygons_ll:
            rings = []
            for ring in polygon:
                pts = [(pt[0], pt[1]) for pt in ring if len(pt) >= 2]
                if len(pts) < 3:
                    continue
                if not same_pt(pts[0], pts[-1]):
                    pts.append(pts[0])
                rings.append(pts)
            if rings:
                polygons.append(rings)
        if not polygons:
            continue

        polygons_xy = [
            [[to_xy(lon, lat) for lon, lat in ring] for ring in polygon]
            for polygon in polygons
        ]
        area = 0.0
        cx = cy = 0.0
        for polygon in polygons_xy:
            (px, py), poly_area = polygon_centroid(polygon)
            cx += px * poly_area
            cy += py * poly_area
            area += poly_area
        if area > 1e-9:
            cx /= area
            cy /= area
        else:
            vertices = [pt for polygon in polygons_xy for ring in polygon for pt in ring]
            cx = sum(pt[0] for pt in vertices) / len(vertices)
            cy = sum(pt[1] for pt in vertices) / len(vertices)

        points = [(cx, cy)]
        for polygon in polygons_xy:
            for ring in polygon:
                points.extend(ring)
        distance = boundary.min_distance(points)

        new_props = {
            "name": name,
            "distance_m": int(round(distance)),
            "area_m2": int(round(area)),
            "osm_type": "manual",
            "osm_id": f"nearby-{index}",
            "city": str(props.get("city") or "").strip(),
            "note": "市外・小平市境界から500m以内（手入力ポリゴン）",
            "name_source": "manual",
        }
        existing = by_name.get(name)
        if existing:
            existing["properties"].update(new_props)
            existing["geometry"] = geom
        else:
            feature = {"type": "Feature", "properties": new_props, "geometry": geom}
            features.append(feature)
            by_name[name] = feature
        count += 1
    return count


def geojson_geometry(polygons):
    def out_ring(ring):
        return [[round(x, 7), round(y, 7)] for x, y in ring]

    if len(polygons) == 1:
        return {"type": "Polygon", "coordinates": [out_ring(r) for r in polygons[0]]}
    return {
        "type": "MultiPolygon",
        "coordinates": [[out_ring(r) for r in polygon] for polygon in polygons],
    }


def main():
    print("Overpass API から取得中...")
    elements = fetch_overpass()

    # 名前のない公園ポリゴンに名前を付けるための、公園らしい名前のノード。
    named_nodes = []
    park_feature_names = set()
    for el in elements:
        if el.get("type") == "node":
            tags = el.get("tags") or {}
            name = (tags.get("name") or "").strip()
            lat = el.get("lat")
            lon = el.get("lon")
            if len(name) < 3 or not PARK_NAME_RE.search(name):
                continue
            if lat is None or lon is None:
                continue
            named_nodes.append((name, lon, lat))
        else:
            name = ((el.get("tags") or {}).get("name") or "").strip()
            if name:
                park_feature_names.add(name)

    boundary_ll = load_boundary()
    lons = [lon for polygon in boundary_ll for ring in polygon for lon, _ in ring]
    lats = [lat for polygon in boundary_ll for ring in polygon for _, lat in ring]
    to_xy = projector(sum(lats) / len(lats), sum(lons) / len(lons))
    boundary_xy = [
        [[to_xy(lon, lat) for lon, lat in ring] for ring in polygon]
        for polygon in boundary_ll
    ]
    boundary = BoundaryIndex(boundary_xy)

    skip = existing_keys()
    exclusions = exclusion_polygons()
    features = []
    skipped_inside = 0
    used_names = set(park_feature_names)
    for el in elements:
        if el.get("type") == "node":
            continue
        if (el.get("type"), el.get("id")) in skip:
            continue

        polygons_ll = build_polygons(el)
        if not polygons_ll:
            continue

        # 小金井公園の内側にある広場・公園は、独立した近隣公園にしない。
        if exclusions:
            vertices = [pt for polygon in polygons_ll for pt in polygon[0]]
            cx_ll = sum(pt[0] for pt in vertices) / len(vertices)
            cy_ll = sum(pt[1] for pt in vertices) / len(vertices)
            if any(polygon_contains(poly, cx_ll, cy_ll) for poly in exclusions):
                continue

        tags = el.get("tags") or {}
        name = (tags.get("name") or "").strip()
        name_source = "osm"
        if not name:
            name = name_from_named_node(polygons_ll, named_nodes, used_names)
            name_source = "nearby_node"
        if not name:
            continue
        used_names.add(name)
        polygons = [
            [[to_xy(lon, lat) for lon, lat in ring] for ring in polygon]
            for polygon in polygons_ll
        ]

        area = 0.0
        cx = cy = 0.0
        for polygon in polygons:
            (px, py), poly_area = polygon_centroid(polygon)
            cx += px * poly_area
            cy += py * poly_area
            area += poly_area
        if area > 1e-9:
            cx /= area
            cy /= area
        else:
            vertices = [pt for polygon in polygons for ring in polygon for pt in ring]
            cx = sum(pt[0] for pt in vertices) / len(vertices)
            cy = sum(pt[1] for pt in vertices) / len(vertices)

        if boundary.contains(cx, cy):
            skipped_inside += 1
            continue

        points = [(cx, cy)]
        for polygon in polygons:
            for ring in polygon:
                points.extend(ring)
        distance = boundary.min_distance(points)
        if distance > MAX_DIST_M:
            continue

        features.append(
            {
                "type": "Feature",
                "properties": {
                    "name": name,
                    "distance_m": int(round(distance)),
                    "area_m2": int(round(area)),
                    "osm_type": el.get("type"),
                    "osm_id": el.get("id"),
                    "city": (tags.get("addr:city") or "").strip(),
                    "note": NOTE,
                    "name_source": name_source,
                },
                "geometry": geojson_geometry(polygons_ll),
            }
        )

    existing_names = {f["properties"]["name"] for f in features}
    for manual in MANUAL_FEATURES:
        if manual["name"] in existing_names:
            continue
        features.append(
            {
                "type": "Feature",
                "properties": {
                    "name": manual["name"],
                    "distance_m": manual["distance_m"],
                    "area_m2": manual["area_m2"],
                    "osm_type": manual["osm_type"],
                    "osm_id": manual["osm_id"],
                    "city": manual["city"],
                    "note": manual["note"],
                    "name_source": "manual",
                },
                "geometry": {
                    "type": "Point",
                    "coordinates": [manual["lon"], manual["lat"]],
                },
            }
        )

    manual_count = merge_nearby_manual(features, boundary, to_xy)

    features.sort(key=lambda f: (f["properties"]["distance_m"], f["properties"]["name"]))
    out = {"type": "FeatureCollection", "features": features}
    with open(OUT, "w", encoding="utf-8") as fp:
        json.dump(out, fp, ensure_ascii=False, separators=(",", ":"))
        fp.write("\n")

    if features:
        dists = [f["properties"]["distance_m"] for f in features]
        print(
            f"features {len(features)} (distance {min(dists)}-{max(dists)}m, "
            f"市内で除外 {skipped_inside}, 手入力 {manual_count}) -> {OUT}"
        )
    else:
        print(f"features 0 (市内で除外 {skipped_inside}, 手入力 {manual_count}) -> {OUT}")


if __name__ == "__main__":
    main()