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

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以内"

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});'
        ");"
        "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 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 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()

    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()
    features = []
    skipped_inside = 0
    for el in elements:
        if (el.get("type"), el.get("id")) in skip:
            continue
        tags = el.get("tags") or {}
        name = (tags.get("name") or "").strip()
        if not name:
            continue

        polygons_ll = build_polygons(el)
        if not polygons_ll:
            continue
        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,
                },
                "geometry": geojson_geometry(polygons_ll),
            }
        )

    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}) -> {OUT}"
        )
    else:
        print(f"features 0 (市内で除外 {skipped_inside}) -> {OUT}")


if __name__ == "__main__":
    main()