path: root/scripts/import-manual-parks.py

#!/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()