parks: 公園一覧と地図、公共施設の情報を更新

- 近隣公園のクリックを市内公園と同じ挙動にし、Googleマップを公園名で開く
- 公共施設のトイレを○にし、ルネこだいらを追加
- 未取得だった施設の座標を追加し、再ジオコーディング時の上書きを防止
- 手入力ポリゴン、近隣公園、PWA、ライセンス表示を追加

1 files changed, 378 insertions, 0 deletions

diff --git a/scripts/update-nearby-parks.py b/scripts/update-nearby-parks.py
new file mode 100644
index 0000000..39268cb
--- /dev/null
+++ b/scripts/update-nearby-parks.py
@@ -0,0 +1,378 @@
+#!/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()