diff options
Diffstat (limited to 'scripts')
| -rw-r--r-- | scripts/geocode-facilities.py | 86 | ||||
| -rw-r--r-- | scripts/import-manual-parks.py | 445 | ||||
| -rw-r--r-- | scripts/make-park-qgis-kit.py | 278 | ||||
| -rw-r--r-- | scripts/update-nearby-parks.py | 378 | ||||
| -rw-r--r-- | scripts/update-osm-parks.py | 266 |
5 files changed, 1431 insertions, 22 deletions
diff --git a/scripts/geocode-facilities.py b/scripts/geocode-facilities.py new file mode 100644 index 0000000..12c727c --- /dev/null +++ b/scripts/geocode-facilities.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python3 +"""小平市の公共施設JSONに緯度・経度を付与する（Nominatim）。 + +使い方: + py -X utf8 scripts/geocode-facilities.py + +- `src/data/kodaira-facilities.json` の各施設に lat / lon / geocode を追加する。 +- 住所で見つからない場合は「小平市 + 施設名」で再試行する。 +- 見つからない場合は lat / lon を null にし、note に理由を追記する。 +- Nominatimの利用ポリシーに従い、1リクエストあたり約1.1秒待つ。 +""" +import json +import time +import urllib.parse +import urllib.request + +PATH = "src/data/kodaira-facilities.json" +UA = "yasutakeyohei.com facilities geocoder (park list)" + +# 住所・名称検索では正しい位置にならない施設。手動の座標を優先する。 +# 値は ((緯度, 経度), geocodeラベル)。 +MANUAL = { + "小平市役所（本庁）": ((35.7279025, 139.4774885), "手動補正"), + "西部出張所": ((35.735038, 139.461827), "小川西町図書館と同じ建物として補正"), + "小川分室": ((35.72604, 139.459319), "小川公民館と同じ建物として補正"), + "東部出張所": ((35.7279272, 139.5140803), "東部市民センターと同じ建物として補正"), + "喜平図書館": ((35.7173213, 139.4908793), "名称"), +} + + +def geocode(query): + url = ( + "https://nominatim.openstreetmap.org/search?format=jsonv2&limit=1&countrycodes=jp&q=" + + urllib.parse.quote(query) + ) + req = urllib.request.Request(url, headers={"User-Agent": UA}) + with urllib.request.urlopen(req, timeout=30) as r: + data = json.loads(r.read().decode("utf-8")) + if not data: + return None + return float(data[0]["lat"]), float(data[0]["lon"]) + + +def main(): + data = json.load(open(PATH, encoding="utf-8")) + for i, f in enumerate(data, 1): + name = f.get("name", "") + address = f.get("address", "") + if name in MANUAL: + (lat, lon), label = MANUAL[name] + f["lat"], f["lon"] = lat, lon + f["geocode"] = label + f["note"] = str(f.get("note", "")).replace("地図の位置は未取得。", "").strip() + print(f"{i}/{len(data)} {name}: {lat:.6f},{lon:.6f} ({label})") + continue + found = None + source = "" + if address: + found = geocode("東京都" + address) + source = "住所" + if not found and name: + time.sleep(1.1) + found = geocode("小平市" + name) + source = "名称" + if found and found[0] is not None: + f["lat"], f["lon"] = found + f["geocode"] = source + f["note"] = str(f.get("note", "")).replace("地図の位置は未取得。", "").strip() + print(f"{i}/{len(data)} {name}: {found[0]:.6f},{found[1]:.6f} ({source})") + else: + f["lat"], f["lon"] = None, None + f["geocode"] = "" + note = str(f.get("note", "")).strip() + if "地図の位置は未取得。" not in note: + f["note"] = (note + " " if note else "") + "地図の位置は未取得。" + print(f"{i}/{len(data)} {name}: NOT FOUND") + time.sleep(1.1) + with open(PATH, "w", encoding="utf-8") as fh: + json.dump(data, fh, ensure_ascii=False, indent=2) + fh.write("\n") + ok = sum(1 for f in data if f.get("lat") is not None) + print(f"geocoded {ok} / {len(data)} -> {PATH}") + + +if __name__ == "__main__": + main() diff --git a/scripts/import-manual-parks.py b/scripts/import-manual-parks.py new file mode 100644 index 0000000..7eb3f1c --- /dev/null +++ b/scripts/import-manual-parks.py @@ -0,0 +1,445 @@ +#!/usr/bin/env python3 +"""QGISで描いた公園ポリゴン(GeoJSON)を kodaira-parks-manual.json に取り込む。 + +使い方: + py -X utf8 scripts/import-manual-parks.py reference/qgis/manual-parks.geojson --dry-run + py -X utf8 scripts/import-manual-parks.py reference/qgis/manual-parks.geojson + +主な仕様: + - 公園名は properties の `our_name` を優先。なければ `name`、それもなければ + 最寄りの市の公園（既定60m以内・曖昧でない場合）から自動判定する + - 1つの公園を複数ポリゴンに分けた場合は、同じ `our_name` で複数Featureを入れる + - 同じ公園名の既存Featureは置き換える（追記したい場合は --append） + - 面積が --min-area 未満のポリゴンは下書きの消し忘れとみなし、取り込まない + - OSM側に既にポリゴンがある公園名は、二重描画を防ぐため既定で拒否（--force で許可） +""" + +import argparse +import difflib +import json +import math +import os +import sys +import unicodedata + +ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) +OURS = os.path.join(ROOT, "src", "data", "kodaira-parks.json") +OSM = os.path.join(ROOT, "src", "data", "kodaira-parks-osm.json") +MANUAL = os.path.join(ROOT, "src", "data", "kodaira-parks-manual.json") +DEFAULT_INPUT = os.path.join(ROOT, "reference", "qgis", "manual-parks.geojson") + +ID_BASE = 900000 + +# 描いた面積と台帳面積の倍率。この範囲を外れたら警告する +AREA_RATIO_WARN_LOW = 0.7 +AREA_RATIO_WARN_HIGH = 1.5 + +# 小平市の範囲（余裕を持たせた bounding box） +LAT_MIN, LAT_MAX = 35.66, 35.79 +LON_MIN, LON_MAX = 139.40, 139.55 + +if hasattr(sys.stdout, "reconfigure"): + sys.stdout.reconfigure(encoding="utf-8") + + +def norm(s): + s = "".join(unicodedata.normalize("NFKC", str(s or "")).split()) + return s.replace("ヶ", "ケ").replace("ツ", "ッ") + + +def load(path): + with open(path, encoding="utf-8") as f: + return json.load(f) + + +def rings(geom): + """外周リングの一覧。空のジオメトリでも落ちないようにする。""" + if not geom: + return [] + coords = geom.get("coordinates") or [] + if geom.get("type") == "Polygon": + return [coords[0]] if coords and coords[0] else [] + if geom.get("type") == "MultiPolygon": + return [poly[0] for poly in coords if poly and poly[0]] + return [] + + +def close_ring(ring): + """最初と最後の点が一致していなければ閉じる。(ring, 閉じたか)""" + if not ring: + return ring, False + if ring[0] != ring[-1]: + return ring + [ring[0]], True + return ring, False + + +def area_m2(geom): + total = 0.0 + for ring in rings(geom): + if len(ring) < 4: + continue + lat0 = sum(p[1] for p in ring) / len(ring) + kx = 111320 * math.cos(math.radians(lat0)) + ky = 111320 + s = 0.0 + for (x1, y1), (x2, y2) in zip(ring, ring[1:]): + s += (x1 * kx) * (y2 * ky) - (x2 * kx) * (y1 * ky) + total += abs(s) / 2 + return total + + +def ring_area_centroid(ring): + a = cx = cy = 0.0 + for (x1, y1), (x2, y2) in zip(ring, ring[1:]): + cross = x1 * y2 - x2 * y1 + a += cross + cx += (x1 + x2) * cross + cy += (y1 + y2) * cross + if a == 0: + return 0.0, ring[0][0], ring[0][1] + a /= 2 + return abs(a), cx / (6 * a), cy / (6 * a) + + +def centroid(geom): + """面積で重み付けした重心。退化している場合は最初の点を返す。""" + parts = [] + for ring in rings(geom): + area, cx, cy = ring_area_centroid(ring) + parts.append((area, cx, cy)) + total = sum(p[0] for p in parts) + if total == 0: + pts = [p for ring in rings(geom) for p in ring] + if not pts: + return None, None + return pts[0][1], pts[0][0] + return sum(a * cy for a, _, cy in parts) / total, sum(a * cx for a, cx, _ in parts) / total + + +def dist_m(lat1, lon1, lat2, lon2): + dy = (lat1 - lat2) * 111320 + dx = (lon1 - lon2) * 111320 * math.cos(math.radians((lat1 + lat2) / 2)) + return math.hypot(dx, dy) + + +def nearest_park(lat, lon, parks, max_distance): + cand = sorted( + ((dist_m(lat, lon, q["lat"], q["lon"]), q["name"]) for q in parks), + key=lambda x: x[0], + ) + if not cand: + return None, "公園一覧に座標がありません" + d1, n1 = cand[0] + if d1 > max_distance: + return None, f"最寄りの公園 {n1} まで {d1:.0f}m（--max-distance {max_distance:.0f}m 超）" + if len(cand) > 1: + d2, n2 = cand[1] + if d2 - d1 < 20 and d2 < d1 * 1.8: + return None, f"{n1} {d1:.0f}m と {n2} {d2:.0f}m が近く、判別できません" + return n1, "" + + +def resolve_name(props, lat, lon, parks, by_norm, max_distance, use_auto): + """(公園名, 判定方法, 補足メッセージ) を返す。判定できなければ (None, ...)。""" + raw = str(props.get("our_name") or "").strip() + if raw: + if norm(raw) in by_norm: + return by_norm[norm(raw)], "our_name", "" + close = difflib.get_close_matches(norm(raw), list(by_norm), n=1, cutoff=0.8) + hint = f"（候補: {by_norm[close[0]]}）" if close else "" + return None, "our_name", f"一覧にない公園名: {raw}{hint}" + + raw = str(props.get("name") or "").strip() + if raw: + if norm(raw) in by_norm: + return by_norm[norm(raw)], "name", "" + close = difflib.get_close_matches(norm(raw), list(by_norm), n=1, cutoff=0.9) + if close: + return by_norm[close[0]], "name(類似)", "" + return None, "name", f"一覧にない公園名: {raw}" + + if use_auto and lat is not None: + name, why = nearest_park(lat, lon, parks, max_distance) + if name: + return name, "最寄り", "" + return None, "最寄り", why + return None, "-", "our_name / name がなく、自動判定も無効" + + +def round_coords(coords, digits=7): + if isinstance(coords, (int, float)): + return round(coords, digits) + return [round_coords(c, digits) for c in coords] + + +def iter_features(obj): + if obj.get("type") == "FeatureCollection": + return obj.get("features") or [] + if obj.get("type") == "Feature": + return [obj] + raise SystemExit("GeoJSON の Feature / FeatureCollection ではありません") + + +def main(): + ap = argparse.ArgumentParser(description="QGISで描いた公園ポリゴンの取り込み") + ap.add_argument("input", nargs="?", default=DEFAULT_INPUT, help="取り込む GeoJSON") + ap.add_argument("--out", default=MANUAL, help="出力先（既定: src/data/kodaira-parks-manual.json）") + ap.add_argument("--dry-run", action="store_true", help="書き込まずに結果だけ表示") + ap.add_argument("--append", action="store_true", help="同名の既存Featureを残して追記する") + ap.add_argument("--allow-partial", action="store_true", help="判定できないFeatureがあっても、判定できた分だけ書き込む") + ap.add_argument("--no-auto-name", action="store_true", help="最寄り公園からの自動判定を使わない") + ap.add_argument("--max-distance", type=float, default=60.0, help="自動判定に使う最寄り公園までの距離(m)") + ap.add_argument("--min-area", type=float, default=20.0, help="これ未満の面積は下書きとみなす(㎡)") + ap.add_argument("--max-area", type=float, default=200000.0, help="これを超える面積は警告(㎡)") + ap.add_argument("--force", action="store_true", help="OSM側に既にポリゴンがある公園も取り込む") + ap.add_argument( + "--force-name", + action="append", + default=[], + metavar="公園名", + help="指定した公園だけ、OSM側にポリゴンがあっても取り込む（複数指定可）", + ) + ap.add_argument("--update-parks", action="store_true", help="kodaira-parks.json の空の緯度・経度を重心で埋める") + args = ap.parse_args() + force_names = set(args.force_name or []) + + if not os.path.exists(args.input): + sys.exit(f"入力が見つかりません: {args.input}") + + ours = load(OURS) + osm = load(OSM) + manual = load(args.out) if os.path.exists(args.out) else {"type": "FeatureCollection", "features": []} + + parks = [] + by_norm = {} + official_areas = {} + for r in ours: + name = str(r.get("公園名") or "").strip() + if not name: + continue + by_norm[norm(name)] = name + try: + official_areas[name] = float(r.get("面積㎡") or "") + except ValueError: + pass + try: + parks.append({"name": name, "lat": float(r["緯度"]), "lon": float(r["経度"])}) + except (KeyError, ValueError): + pass + + osm_names = { + f["properties"].get("our_name") + for f in osm.get("features", []) + if f["properties"].get("matched") and f["properties"].get("our_name") + } + # OSM側のポリゴン面積（手入力と合わせて台帳面積と比べる） + osm_areas = {} + for f in osm.get("features", []): + p = f.get("properties", {}) + name = p.get("our_name") + if p.get("matched") and name: + osm_areas[name] = osm_areas.get(name, 0.0) + area_m2(f["geometry"]) + manual_features = manual.get("features") or [] + manual_names = {f["properties"].get("our_name") for f in manual_features} + + errors = [] + warnings = [] + resolved = [] # dict(name, method, geom, area, lat, lon, props) + + for i, feat in enumerate(iter_features(load(args.input)), start=1): + props = feat.get("properties") or {} + geom = feat.get("geometry") + if not geom or geom.get("type") not in ("Polygon", "MultiPolygon"): + errors.append(f"#{i}: ポリゴンではありません（{geom.get('type') if geom else 'geometry なし'}）") + continue + if str(props.get("draft") or "").lower() in ("true", "1", "yes"): + errors.append(f"#{i}: draft:true のFeatureが残っています（下書きを削除するか属性を消してください）") + continue + if not geom.get("coordinates"): + errors.append( + f"#{i}: ジオメトリが空です（fid {props.get('fid', '?')}, " + f"{props.get('our_name', '名前なし')}）。未完成のポリゴンを削除してください" + ) + continue + + # リングを閉じる + changed = 0 + if geom["type"] == "Polygon": + new_rings = [] + for idx, ring in enumerate(geom["coordinates"]): + fixed, did = close_ring(ring) + changed += did + new_rings.append(fixed) + geom["coordinates"] = new_rings + else: + new_polys = [] + for poly in geom["coordinates"]: + new_rings = [] + for ring in poly: + fixed, did = close_ring(ring) + changed += did + new_rings.append(fixed) + new_polys.append(new_rings) + geom["coordinates"] = new_polys + if changed: + warnings.append(f"#{i}: 閉じていなかったリングを{changed}個補正しました") + + # 座標系のチェック + pts = [p for ring in rings(geom) for p in ring] + if not pts: + errors.append(f"#{i}: 座標が空です") + continue + if any(abs(p[0]) > 180 or abs(p[1]) > 90 for p in pts): + errors.append(f"#{i}: 緯度経度ではない座標です。EPSG:4326 で書き出してください") + continue + lat, lon = centroid(geom) + if not (LAT_MIN <= lat <= LAT_MAX and LON_MIN <= lon <= LON_MAX): + warnings.append( + f"#{i}: 重心が小平市の範囲外です ({lat:.5f}, {lon:.5f})。座標系を確認してください" + ) + + name, method, why = resolve_name( + props, lat, lon, parks, by_norm, args.max_distance, not args.no_auto_name + ) + if not name: + errors.append(f"#{i}: 公園名を判定できません — {why}") + continue + + area = area_m2(geom) + if area < args.min_area: + errors.append( + f"#{i}: {name} の面積が {area:.1f}㎡ です。下書きの消し忘れか、形状を確認してください" + ) + continue + if area > args.max_area: + warnings.append(f"#{i}: {name} の面積が {area:.1f}㎡ と大きすぎます") + forced = args.force or name in force_names + if name in osm_names and name not in manual_names and not forced: + errors.append( + f"#{i}: {name} はOSM側に既にポリゴンがあります" + "（--force または --force-name で強制取り込み）" + ) + continue + + resolved.append( + {"name": name, "method": method, "geom": geom, "area": area, "lat": lat, "lon": lon} + ) + + if errors and not args.allow_partial: + print("取り込みを中止しました。次の点を直してください。\n") + for e in errors: + print(f" - {e}") + if warnings: + print("\n警告:") + for w in warnings: + print(f" - {w}") + sys.exit(1) + if errors and args.allow_partial: + print("次のFeatureは取り込みませんでした:\n") + for e in errors: + print(f" - {e}") + print() + + # 公園名ごとにまとめる（複数ポリゴン可） + grouped = {} + for item in resolved: + grouped.setdefault(item["name"], []).append(item) + + existing_ids = {} + for f in manual_features: + p = f["properties"] + try: + osm_id = int(p.get("osm_id")) + except (TypeError, ValueError): + continue + if p.get("our_name"): + existing_ids[p["our_name"]] = osm_id + next_id = max([ID_BASE, *existing_ids.values()]) + 1 + + new_features = [] + for name, items in grouped.items(): + osm_id = existing_ids.get(name) + if osm_id is None: + osm_id = next_id + next_id += 1 + for item in items: + geom = item["geom"] + geom["coordinates"] = round_coords(geom["coordinates"]) + new_features.append( + { + "type": "Feature", + "properties": { + "name": "", + "our_name": name, + "matched": True, + "osm_type": "manual", + "osm_id": osm_id, + "leisure": "park", + }, + "geometry": geom, + } + ) + + touched = set(grouped) + if args.append: + kept = manual_features + else: + kept = [f for f in manual_features if f["properties"].get("our_name") not in touched] + merged = kept + new_features + merged.sort(key=lambda f: (f["properties"].get("osm_id") or 0)) + + print(f"取り込み: {len(new_features)} Feature / {len(grouped)} 公園") + print(f" 置き換え: {len(touched & manual_names)} 公園 / 新規: {len(touched - manual_names)} 公園") + for name, items in sorted(grouped.items()): + method = items[0]["method"] + area = sum(i["area"] for i in items) + official = official_areas.get(name) + osm_area = osm_areas.get(name, 0.0) + combined = area + osm_area + detail = "" + if official: + ratio = combined / official + extra = f" + OSM {osm_area:.0f}㎡ = {combined:.0f}㎡" if osm_area else "" + detail = f"、台帳 {official:.0f}㎡（{ratio:.2f}倍{extra}）" + if ratio < AREA_RATIO_WARN_LOW or ratio > AREA_RATIO_WARN_HIGH: + warnings.append( + f"{name}: 面積が台帳の {ratio:.2f} 倍です（{combined:.0f}㎡ / {official:.0f}㎡）。描き直しを検討してください" + ) + print(f" - {name}: {len(items)}ポリゴン, {area:.0f}㎡{detail}, 判定={method}") + if warnings: + print("\n警告:") + for w in warnings: + print(f" - {w}") + + if args.dry_run: + print("\n--dry-run のため書き込んでいません。") + return + + with open(args.out, "w", encoding="utf-8") as f: + json.dump({"type": "FeatureCollection", "features": merged}, f, ensure_ascii=False, indent=2) + f.write("\n") + + if args.update_parks and grouped: + updated = 0 + for r in ours: + name = r.get("公園名") + if name not in grouped: + continue + lat = grouped[name][0]["lat"] + lon = grouped[name][0]["lon"] + if not str(r.get("緯度") or "").strip(): + r["緯度"] = f"{lat:.7f}" + updated += 1 + if not str(r.get("経度") or "").strip(): + r["経度"] = f"{lon:.7f}" + updated += 1 + if updated: + with open(OURS, "w", encoding="utf-8") as f: + json.dump(ours, f, ensure_ascii=False, separators=(",", ":")) + print(f"kodaira-parks.json の緯度・経度を {updated} 項目更新しました") + + print(f"\n書き込み: {args.out}") + + +if __name__ == "__main__": + main() diff --git a/scripts/make-park-qgis-kit.py b/scripts/make-park-qgis-kit.py new file mode 100644 index 0000000..0df6d5a --- /dev/null +++ b/scripts/make-park-qgis-kit.py @@ -0,0 +1,278 @@ +#!/usr/bin/env python3 +"""QGISで公園ポリゴンを手描きするための下地データを出力する。 + +使い方: + py -X utf8 scripts/make-park-qgis-kit.py + +出力（reference/ は .gitignore 対象）: + reference/qgis/parks-without-polygon.geojson … ポリゴンが未整備の市の公園（点） + reference/qgis/existing-parks.geojson … 整備済みポリゴン（参考・重複防止） + reference/qgis/README.md … QGISでの手順 + +描いたポリゴンは scripts/import-manual-parks.py で +src/data/kodaira-parks-manual.json に取り込む。 +""" + +import argparse +import json +import math +import os +import re +import sys +import unicodedata + +ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) +OURS = os.path.join(ROOT, "src", "data", "kodaira-parks.json") +OSM = os.path.join(ROOT, "src", "data", "kodaira-parks-osm.json") +MANUAL = os.path.join(ROOT, "src", "data", "kodaira-parks-manual.json") +OUT_DIR = os.path.join(ROOT, "reference", "qgis") +MANUAL_QGIS = os.path.join(OUT_DIR, "manual-parks.geojson") + +if hasattr(sys.stdout, "reconfigure"): + sys.stdout.reconfigure(encoding="utf-8") + + +def load(path): + with open(path, encoding="utf-8") as f: + return json.load(f) + + +def norm(s): + s = "".join(unicodedata.normalize("NFKC", str(s or "")).split()) + return s.replace("ヶ", "ケ").replace("ツ", "ッ") + + +def has_coords(geom): + """空のポリゴン（QGISの未完成Feature）を弾く。""" + if not geom: + return False + coords = geom.get("coordinates") or [] + if geom.get("type") == "Polygon": + return bool(coords and coords[0]) + if geom.get("type") == "MultiPolygon": + return any(poly and poly[0] for poly in coords) + return False + + +def center(geom): + if geom["type"] == "Polygon": + pts = geom["coordinates"][0] + else: + pts = [p for poly in geom["coordinates"] for p in poly[0]] + return ( + (min(p[1] for p in pts) + max(p[1] for p in pts)) / 2, + (min(p[0] for p in pts) + max(p[0] for p in pts)) / 2, + ) + + +def dist_m(a, b): + dy = (a[0] - b[0]) * 111320 + dx = (a[1] - b[1]) * 111320 * math.cos(math.radians((a[0] + b[0]) / 2)) + return math.hypot(dx, dy) + + +def write_geojson(path, features): + with open(path, "w", encoding="utf-8") as f: + json.dump({"type": "FeatureCollection", "features": features}, f, ensure_ascii=False) + f.write("\n") + + +def main(): + ap = argparse.ArgumentParser(description="QGISで公園ポリゴンを手描きするための下地データを出力する") + ap.add_argument("--list", action="store_true", help="ポリゴン未整備の公園一覧も表示する") + args = ap.parse_args() + + ours = load(OURS) + osm = load(OSM) + manual = load(MANUAL) + + # QGISで描いたが、まだ取り込んでいないポリゴンの公園名 + drawn_names = set() + if os.path.exists(MANUAL_QGIS): + try: + for feat in load(MANUAL_QGIS).get("features", []): + p = feat.get("properties") or {} + if not has_coords(feat.get("geometry")): + continue + name = str(p.get("our_name") or "").strip() + if name: + drawn_names.add(norm(name)) + except (ValueError, OSError) as e: + print(f" warning: {MANUAL_QGIS} を読めません: {e}", file=sys.stderr) + + polygons = [] + covered = set() + osm_by_id = {} + osm_unmatched_centers = [] + for feat in osm.get("features", []): + p = feat.get("properties", {}) + osm_by_id[int(p.get("osm_id") or 0)] = p.get("our_name") or "" + if p.get("matched") and p.get("our_name"): + covered.add(p["our_name"]) + polygons.append( + { + "type": "Feature", + "properties": { + "our_name": p["our_name"], + "source": f"OSM {p.get('osm_type', '')}/{p.get('osm_id', '')}", + }, + "geometry": feat["geometry"], + } + ) + else: + osm_unmatched_centers.append((center(feat["geometry"]), p)) + for feat in manual.get("features", []): + p = feat.get("properties", {}) + if p.get("our_name"): + covered.add(p["our_name"]) + polygons.append( + { + "type": "Feature", + "properties": { + "our_name": p["our_name"], + "source": f"manual {p.get('osm_id', '')}", + }, + "geometry": feat["geometry"], + } + ) + + points = [] + no_point = [] + already = [] + for r in ours: + name = r.get("公園名", "") + if not name or name in covered: + continue + # 「緑地（OSM 1234567）」は同じIDのポリゴンが既にある + m = re.match(r"^緑地（OSM (\d+)）$", name) + if m and int(m.group(1)) in osm_by_id: + already.append(name) + continue + try: + lat = float(r.get("緯度") or "") + lon = float(r.get("経度") or "") + except ValueError: + no_point.append(name) + continue + # 名前未設定のOSMポリゴンが同じ場所にある場合も、二重描き防止のため除外 + if any(dist_m((lat, lon), c) <= 20 for c, _ in osm_unmatched_centers): + already.append(name) + continue + points.append( + { + "type": "Feature", + "properties": { + "no": len(points) + 1, + "our_name": name, + "address": r.get("住所", ""), + "kind": r.get("種別", ""), + "official_area_m2": r.get("面積㎡", ""), + "drawn": norm(name) in drawn_names, + "lat": lat, + "lon": lon, + }, + "geometry": {"type": "Point", "coordinates": [lon, lat]}, + } + ) + + counts = { + "covered": len(covered), + "missing": len(points) + len(no_point), + "no_point": len(no_point), + "already": len(already), + "drawn": sum(1 for p in points if p["properties"]["drawn"]), + } + + os.makedirs(OUT_DIR, exist_ok=True) + write_geojson(os.path.join(OUT_DIR, "parks-without-polygon.geojson"), points) + write_geojson(os.path.join(OUT_DIR, "existing-parks.geojson"), polygons) + readme = f"""# QGIS で公園ポリゴンを描く + +`scripts/make-park-qgis-kit.py` の生成物です。再生成すると上書きされます。 + +- ポリゴン整備済み: {counts['covered']} 件 +- OSM側に既にポリゴンあり（未リンク）: {counts['already']} 件 +- ポリゴン未整備: {counts['missing']} 件（うち手入力済み・未取り込み {counts['drawn']} 件、座標なし {counts['no_point']} 件） + +## 1. QGIS に背景地図を追加 + +「レイヤ > レイヤを追加 > XYZ タイルを追加」で次を登録します。 + +| 用途 | URL | 最大ズーム | +|---|---|---| +| 国土地理院 空中写真（推奨） | `https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{{z}}/{{x}}/{{y}}.jpg` | 18 | +| 国土地理院 淡色地図 | `https://cyberjapandata.gsi.go.jp/xyz/pale/{{z}}/{{x}}/{{y}}.png` | 18 | +| OpenStreetMap | `https://tile.openstreetmap.org/{{z}}/{{x}}/{{y}}.png` | 19 | + +## 2. 下地を読み込む + +- `parks-without-polygon.geojson` … まだポリゴンがない公園（点）。`our_name` に公園名が入っています +- `existing-parks.geojson` … 整備済みポリゴン。二重描き防止用 + +## 3. ポリゴン用レイヤを作る + +1. 「レイヤ > レイヤを作成 > 新しい GeoJSON レイヤを作成」 +2. ファイル名は `reference/qgis/manual-parks.geojson`、ジオメトリ型は **Polygon**、CRS は **EPSG:4326** +3. フィールドを1つ追加: 名前 `our_name`、型 **文字列** +4. 編集モードにして「ポリゴンを追加」で公園の輪郭を描く + - 点レイヤのラベルを次の式にすると、公園名と台帳面積が同時に出ます + `"our_name" || '\n' || "official_area_m2" || ' ㎡'` + - スナップを有効にすると描きやすい +5. 属性入力欄に `parks-without-polygon.geojson` の `our_name` をそのまま入力 + （1つの公園が複数ポリゴンに分かれる場合は、同じ `our_name` で複数描く） +6. 編集を保存 + +> `our_name` を空欄にしても、取り込み時に最寄りの公園点から自動で判定します +> （既定 60m 以内・曖昧でない場合のみ）。 + +## 4. 取り込む + +```bash +py -X utf8 scripts/import-manual-parks.py reference/qgis/manual-parks.geojson --dry-run +py -X utf8 scripts/import-manual-parks.py reference/qgis/manual-parks.geojson +``` + +`--dry-run` で公園名・面積・台帳面積との倍率・警告を確認してから本実行してください。 +台帳面積の 0.7〜1.5 倍を外れると警告が出ます。 +取り込み後に `npm run build` で地図と表を確認します。 + +## 進捗の管理 + +`parks-without-polygon.geojson` の `drawn` は、`manual-parks.geojson` に同じ `our_name` のポリゴンがあると `true` になります。 + +- 属性テーブルで `drawn` を確認できます +- スタイルを `drawn = true` で灰色にすると、残りが分かりやすくなります +- 取り込み後に再生成すると、取り込んだ公園はこの一覧から消えます + +> GeoPackage（`.gpkg`）で描いている場合は、`manual-parks.geojson` に書き出してから再生成してください。 + +## 注意 + +- 座標系は EPSG:4326（緯度経度）で保存してください。Web メルカトル（EPSG:3857）のまま書き出すと位置がずれます +- 面積が 20㎡ 未満のポリゴンは、下書きの消し忘れとみなして取り込みを中止します +- `existing-parks.geojson` にある公園は描かないでください（重複します） +""" + with open(os.path.join(OUT_DIR, "README.md"), "w", encoding="utf-8") as f: + f.write(readme) + + print(f"ポリゴン整備済み: {counts['covered']} 件") + print(f"OSM側に既にポリゴンあり（未リンク）: {counts['already']} 件") + print(f"ポリゴン未整備: {counts['missing']} 件（座標なし {counts['no_point']} 件）") + if no_point: + print("座標がないため地図に出せない公園（先に kodaira-parks.json の緯度・経度を埋める）:") + for n in no_point: + print(f" - {n}") + print(f"出力: {OUT_DIR}") + + if args.list: + print("\nポリゴン未整備の公園:") + for p in points: + props = p["properties"] + mark = " [手入力済み]" if props["drawn"] else "" + print(f" {props['no']:3d}. {props['our_name']}（{props['address'] or '住所未登録'}）{mark}") + for name in no_point: + print(f" --- {name}（緯度・経度未登録）") + + +if __name__ == "__main__": + main() 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() diff --git a/scripts/update-osm-parks.py b/scripts/update-osm-parks.py index fffb2f3..fcca883 100644 --- a/scripts/update-osm-parks.py +++ b/scripts/update-osm-parks.py @@ -11,8 +11,11 @@ OpenStreetMapのデータは ODbL 1.0。地図表示時は 「© OpenStreetMap contributors」の帰属表示が必要。 """ +import difflib import json +import math import os +import re import sys import time import unicodedata @@ -30,17 +33,39 @@ OVERRIDES = "src/data/park-name-overrides.json" NORTH_SHIFT_M = 0 LAT_PER_M = 1 / 111320 +# 名前が一致しないポリゴンの自動照合。誤リンクを避けるため保守的に絞る。 +NEAR_DIST_M = 20 # 公園の点がポリゴンからこの距離以内なら同一とみなす +NEAR_NAME_MIN = 0.5 # 近接照合で、名前付きポリゴンに求める名前類似度 +FUZZY_NAME_MIN = 0.8 # 正規化後の名前類似度がこれを超えたら候補（厳密に大） +FUZZY_MAX_DIST_M = 500 # 名前類似だけで照合する場合の距離上限 +# ポリゴン面積と台帳面積の比。極端に違うものは別施設とみなす +AREA_RATIO_MIN = 0.2 +AREA_RATIO_MAX = 5.0 + +# 「緑地（OSM 1234567）」は OSM の ID を名前に含んでいる +OSM_ID_NAME_RE = re.compile(r"^緑地（OSM (\d+)）$") + +# 数字・方位・新旧など「位置の違い」を表す文字。これだけが違う名前は +# 別公園の可能性が高い（例: 七小西第1公園 / 七小西第2公園）。 +POSITIONAL_CHARS = set("0123456789０１２３４５６７８９一二三四五六七八九十第東西南北上中下前後新") + ENDPOINTS = [ "https://overpass-api.de/api/interpreter", "https://overpass.kumi.systems/api/interpreter", "https://overpass.private.coffee/api/interpreter", ] +# leisure=park だけでは取りこぼす児童遊園・広場・緑地も取得する。 +# 照合は下の厳しめの条件（完全一致→近接→名前類似）で行う。 QUERY = f"""[out:json][timeout:90]; area["name"="{CITY}"]["boundary"="administrative"]->.a; ( - way["leisure"="park"](area.a); - relation["leisure"="park"](area.a); + way["leisure"~"^(park|playground|garden|recreation_ground|village_green)$"](area.a); + relation["leisure"~"^(park|playground|garden|recreation_ground|village_green)$"](area.a); + way["landuse"~"^(recreation_ground|village_green)$"](area.a); + relation["landuse"~"^(recreation_ground|village_green)$"](area.a); + way["landuse"="forest"]["name"~"コゲラの森|こげらの森"](area.a); + relation["landuse"="forest"]["name"~"コゲラの森|こげらの森"](area.a); ); out geom; """ @@ -63,10 +88,74 @@ ALIAS = { def norm(s): - s = unicodedata.normalize("NFKC", str(s)).replace(" ", "").replace(" ", "") + s = "".join(unicodedata.normalize("NFKC", str(s or "")).split()) return s.replace("ヶ", "ケ").replace("ツ", "ッ") +def similarity(a, b): + return difflib.SequenceMatcher(None, a, b).ratio() + + +def strip_positional(s): + return "".join(c for c in s if c not in POSITIONAL_CHARS) + + +def rings(geom): + if geom["type"] == "Polygon": + return [geom["coordinates"][0]] + return [poly[0] for poly in geom["coordinates"]] + + +def point_in_ring(lon, lat, ring): + inside = False + for (x1, y1), (x2, y2) in zip(ring, ring[1:]): + if (y1 > lat) != (y2 > lat): + if lon < (x2 - x1) * (lat - y1) / (y2 - y1) + x1: + inside = not inside + return inside + + +def segment_dist_m(lon, lat, a, b): + """点から線分までの距離（メートル）。""" + k = math.cos(math.radians(lat)) + x1, y1 = a[0] * k, a[1] + x2, y2 = b[0] * k, b[1] + px, py = lon * k, lat + dx, dy = x2 - x1, y2 - y1 + if dx == 0 and dy == 0: + t = 0.0 + else: + t = max(0.0, min(1.0, ((px - x1) * dx + (py - y1) * dy) / (dx * dx + dy * dy))) + return math.hypot(px - (x1 + t * dx), py - (y1 + t * dy)) * 111320 + + +def polygon_dist_m(geom, lat, lon): + """公園の代表点からポリゴンまでの距離。内側なら0。""" + rs = rings(geom) + for ring in rs: + if point_in_ring(lon, lat, ring): + return 0.0 + return min( + (segment_dist_m(lon, lat, a, b) for ring in rs for a, b in zip(ring, ring[1:])), + default=float("inf"), + ) + + +def polygon_area_m2(geom): + total = 0.0 + for ring in rings(geom): + if len(ring) < 4: + continue + lat0 = sum(p[1] for p in ring) / len(ring) + kx = 111320 * math.cos(math.radians(lat0)) + ky = 111320 + s = 0.0 + for (x1, y1), (x2, y2) in zip(ring, ring[1:] + ring[:1]): + s += (x1 * kx) * (y2 * ky) - (x2 * kx) * (y1 * ky) + total += abs(s) / 2 + return total + + def fetch(query): last = None for ep in ENDPOINTS: @@ -150,6 +239,41 @@ def main(): our_names = {norm(r["公園名"]): r["公園名"] for r in ours} alias_norm = {norm(k): v for k, v in ALIAS.items()} + # 市の公園の代表点（緯度・経度）と台帳面積 + our_points = [] + our_areas = {} + our_id_names = {} + for r in ours: + name = r["公園名"] + m = OSM_ID_NAME_RE.match(name) + if m: + our_id_names[int(m.group(1))] = name + try: + our_areas[name] = float(r.get("面積㎡") or "") + except ValueError: + pass + try: + lat = float(r.get("緯度") or "") + lon = float(r.get("経度") or "") + except ValueError: + continue + our_points.append((name, lat, lon)) + + # 前回の出力（新規に照合できたもののレポートと、部分取得時の補完に使う） + previous_matched = set() + previous_features = {} + if os.path.exists(OUT): + try: + prev = json.load(open(OUT, encoding="utf-8")) + for f in prev.get("features", []): + p = f.get("properties", {}) + key = f"{p.get('osm_type')}/{p.get('osm_id')}" + previous_features[key] = f + if p.get("matched"): + previous_matched.add(key) + except (ValueError, KeyError): + pass + # OSM ID単位の手動上書き（osm_type/osm_id -> 市の公園名） overrides = {} if os.path.exists(OVERRIDES): @@ -165,38 +289,72 @@ def main(): else: print(f" warning: override city not found: {key} -> {city}", file=sys.stderr) - def match(name): + def exact_match(name): k = norm(name) if k in our_names: - return our_names[k] + return our_names[k], "exact" if k in alias_norm: - return alias_norm[k] + return alias_norm[k], "alias" # 「〇〇児童公園」→「〇〇公園」のような表記差を吸収 k2 = k.replace("児童", "") if k2 in our_names: - return our_names[k2] - return "" + return our_names[k2], "exact" + return "", "" + + def area_ok(park, osm_area): + official = our_areas.get(park) + if not official: + return True + return AREA_RATIO_MIN <= osm_area / official <= AREA_RATIO_MAX + + def fuzzy_match(name, geom): + """ポリゴンが重なる、または名前が十分似ている市の公園を返す。""" + dists = [(polygon_dist_m(geom, q[1], q[2]), q[0]) for q in our_points] + near = sorted(dists) + near_dist, near_park = near[0] + osm_area = polygon_area_m2(geom) + if near_dist <= NEAR_DIST_M: + # 複数の公園が重なる場合はどちらか判定できないので見送る + if len(near) > 1 and near[1][0] <= NEAR_DIST_M: + return "", "", near_dist, 0.0 + # 別名の公園は、ポリゴンが同じ場所にあっても結びつけない + if name and similarity(norm(name), norm(near_park)) < NEAR_NAME_MIN: + return "", "", near_dist, 0.0 + # 極端に面積が違う場合は別施設（大きなグラウンド等）とみなす + if not area_ok(near_park, osm_area): + return "", "", near_dist, 0.0 + return near_park, "near", near_dist, 0.0 + + if not name: + return "", "", near_dist, 0.0 + k = norm(name) + ranked = sorted( + ((similarity(k, norm(park)), dist, park) for dist, park in dists), + key=lambda x: (-x[0], x[1]), + ) + sim, dist, park = ranked[0] + if sim <= FUZZY_NAME_MIN or dist > FUZZY_MAX_DIST_M: + return "", "", dist, sim + # 数字・方位だけの違いは別公園の可能性が高い。重なっている場合のみ採用 + if strip_positional(k) == strip_positional(norm(park)) and dist > NEAR_DIST_M: + return "", "", dist, sim + if not area_ok(park, osm_area): + return "", "", dist, sim + return park, "fuzzy", dist, sim data = fetch(QUERY) + fetched_keys = {f"{e['type']}/{e['id']}" for e in data["elements"]} features = [] applied = 0 - for e in data["elements"]: + covered = set() + new_matches = [] + + # leisure=park を先に処理し、公園本体が広場・遊び場より先に枠を取るようにする + elements = sorted(data["elements"], key=lambda e: (e.get("tags", {}).get("leisure") != "park", e["type"], e["id"])) + for e in elements: tags = e.get("tags", {}) name = tags.get("name", "") key = f"{e['type']}/{e['id']}" - our_name = overrides.get(key) or match(name) - if key in overrides: - applied += 1 - props = { - "name": name, - "our_name": our_name, - "matched": bool(our_name), - "osm_type": e["type"], - "osm_id": e["id"], - } - for k in ["leisure", "addr:city", "addr:street", "addr:housenumber", "opening_hours"]: - if k in tags: - props[k] = tags[k] if e["type"] == "way": ring = [pt(p) for p in e.get("geometry", [])] @@ -221,8 +379,63 @@ def main(): else: geom = {"type": "MultiPolygon", "coordinates": [[r] + inner_rings for r in outer_rings]} + our_name, reason = overrides.get(key, ""), "override" if key in overrides else "" + dist = sim = None + if key in overrides: + applied += 1 + elif not our_name and e["id"] in our_id_names: + # 「緑地（OSM 1234567）」は ID が名前に入っているので直接結びつける + our_name, reason = our_id_names[e["id"]], "osm_id" + elif not our_name and not exact_match(name)[0]: + our_name, reason, dist, sim = fuzzy_match(name, geom) + if not our_name: + our_name, reason = exact_match(name) + + # 公園本体（leisure=park）以外が既に整備済みの公園に重なる場合は、 + # 公園内の遊び場・広場とみなして取り込まない。 + core = tags.get("leisure") == "park" + if our_name in covered and not (core and reason in ("exact", "alias", "override")): + continue + + props = { + "name": name, + "our_name": our_name, + "matched": bool(our_name), + "osm_type": e["type"], + "osm_id": e["id"], + } + if our_name: + props["match"] = reason + for k in ["leisure", "landuse", "addr:city", "addr:street", "addr:housenumber", "opening_hours"]: + if k in tags: + props[k] = tags[k] + + if our_name: + covered.add(our_name) + if key not in previous_matched: + new_matches.append((our_name, key, reason, tags.get("leisure") or tags.get("landuse", ""), dist, sim)) + features.append({"type": "Feature", "properties": props, "geometry": geom}) + # Overpassが部分的な結果を返すと、override対象のFeatureが抜けることがある。 + # その場合は前回のデータを保持して、公園のポリゴンを失わないようにする。 + for key, city in overrides.items(): + if key in fetched_keys: + continue + prev_feat = previous_features.get(key) + if not prev_feat: + continue + feat = json.loads(json.dumps(prev_feat)) + feat["properties"]["our_name"] = city + feat["properties"]["matched"] = True + feat["properties"]["match"] = "override" + features.append(feat) + applied += 1 + print( + f" warning: {key} は今回の取得に含まれないため前回のデータを保持しました", + file=sys.stderr, + ) + features.sort(key=lambda f: (f["properties"].get("name", ""), f["properties"]["osm_type"], f["properties"]["osm_id"])) out = {"type": "FeatureCollection", "features": features} json.dump(out, open(OUT, "w", encoding="utf-8"), ensure_ascii=False, separators=(",", ":")) @@ -232,6 +445,15 @@ def main(): print(f"features {len(features)} (named {named}, matched {matched}) -> {OUT}") if overrides: print(f" overrides applied {applied} / {len(overrides)}") + if new_matches: + print(f" newly matched: {len(new_matches)}") + for our_name, key, reason, tag, dist, sim in sorted(new_matches): + detail = [] + if sim: + detail.append(f"類似{sim:.2f}") + if dist is not None: + detail.append(f"{dist:.0f}m") + print(f" {our_name} <- {key} [{tag}/{reason}] {' '.join(detail)}") boundary = build_boundary() if boundary: |