From a71727538f65ab4f10f5239550be82590f9d47d1 Mon Sep 17 00:00:00 2001 From: Yasutake Yohei <61961825+yasutakeyohei@users.noreply.github.com> Date: Sat, 24 Oct 2026 12:29:24 +0900 Subject: parks: 公園一覧と地図、公共施設の情報を更新 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 近隣公園のクリックを市内公園と同じ挙動にし、Googleマップを公園名で開く - 公共施設のトイレを○にし、ルネこだいらを追加 - 未取得だった施設の座標を追加し、再ジオコーディング時の上書きを防止 - 手入力ポリゴン、近隣公園、PWA、ライセンス表示を追加 --- scripts/update-nearby-parks.py | 378 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 378 insertions(+) create mode 100644 scripts/update-nearby-parks.py (limited to 'scripts/update-nearby-parks.py') 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() -- cgit v1.3.1