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