diff options
Diffstat (limited to 'scripts/update-nearby-parks.py')
| -rw-r--r-- | scripts/update-nearby-parks.py | 243 |
1 files changed, 237 insertions, 6 deletions
diff --git a/scripts/update-nearby-parks.py b/scripts/update-nearby-parks.py index 39268cb..ef30f53 100644 --- a/scripts/update-nearby-parks.py +++ b/scripts/update-nearby-parks.py @@ -20,6 +20,7 @@ OpenStreetMapのデータは ODbL 1.0。地図表示時は import json import math import os +import re import sys import time import urllib.request @@ -29,11 +30,29 @@ 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", @@ -50,6 +69,7 @@ def fetch_overpass(): "(" 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;" ) @@ -141,6 +161,39 @@ def point_in_ring(x, y, ring): 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: @@ -270,6 +323,120 @@ def existing_keys(): } +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] @@ -286,6 +453,25 @@ 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] @@ -297,19 +483,37 @@ def main(): 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"), el.get("id")) in skip: + if el.get("type") == "node": continue - tags = el.get("tags") or {} - name = (tags.get("name") or "").strip() - if not name: + 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 @@ -353,11 +557,38 @@ def main(): "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: @@ -368,10 +599,10 @@ def main(): dists = [f["properties"]["distance_m"] for f in features] print( f"features {len(features)} (distance {min(dists)}-{max(dists)}m, " - f"市内で除外 {skipped_inside}) -> {OUT}" + f"市内で除外 {skipped_inside}, 手入力 {manual_count}) -> {OUT}" ) else: - print(f"features 0 (市内で除外 {skipped_inside}) -> {OUT}") + print(f"features 0 (市内で除外 {skipped_inside}, 手入力 {manual_count}) -> {OUT}") if __name__ == "__main__": |