diff options
Diffstat (limited to 'scripts/import-manual-parks.py')
| -rw-r--r-- | scripts/import-manual-parks.py | 445 |
1 files changed, 445 insertions, 0 deletions
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() |