path: root/scripts/make-park-qgis-kit.py

#!/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")
NEARBY = os.path.join(ROOT, "src", "data", "kodaira-nearby-parks.json")
OUT_DIR = os.path.join(ROOT, "reference", "qgis")
MANUAL_QGIS = os.path.join(OUT_DIR, "manual-parks.geojson")
NEARBY_MANUAL_QGIS = os.path.join(OUT_DIR, "nearby-parks-manual.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)
    nearby = load(NEARBY)

    # 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)

    nearby_drawn_names = set()
    if os.path.exists(NEARBY_MANUAL_QGIS):
        try:
            for feat in load(NEARBY_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:
                    nearby_drawn_names.add(norm(name))
        except (ValueError, OSError) as e:
            print(f"  warning: {NEARBY_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,
                    "source": "city",
                },
                "geometry": {"type": "Point", "coordinates": [lon, lat]},
            }
        )

    # 市外の近隣公園のうち、まだポリゴンがないもの（点）も描画対象に含める。
    nearby_points = 0
    for feat in nearby.get("features", []):
        geom = feat.get("geometry") or {}
        if geom.get("type") != "Point":
            continue
        p = feat.get("properties") or {}
        name = str(p.get("name") or "").strip()
        if not name:
            continue
        lon, lat = geom["coordinates"][:2]
        points.append(
            {
                "type": "Feature",
                "properties": {
                    "no": len(points) + 1,
                    "our_name": name,
                    "address": f"{p.get('city') or '市外'}（市外）",
                    "kind": "市外の近隣公園",
                    "official_area_m2": p.get("area_m2", ""),
                    "drawn": norm(name) in nearby_drawn_names,
                    "lat": lat,
                    "lon": lon,
                    "source": "nearby",
                    "distance_m": p.get("distance_m", ""),
                },
                "geometry": {"type": "Point", "coordinates": [lon, lat]},
            }
        )
        nearby_points += 1

    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"]),
        "nearby_points": nearby_points,
    }

    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']} 件）
- 市外の近隣公園（点）: {counts['nearby_points']} 件

## 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` に公園名が入っています。`source` が `city` なら小平市内、`nearby` なら市外の近隣公園です
- `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 以内・曖昧でない場合のみ）。

### 市外の近隣公園を描く場合

`parks-without-polygon.geojson` の `source` が `nearby` の点は、小平市外の近隣公園です。
小平市内の公園とは別のファイルに描きます。

1. `reference/qgis/nearby-parks-manual.geojson` を新規作成（Polygon / EPSG:4326、`our_name` フィールド）
2. 公園の輪郭を描き、`our_name` に点の名前を入れる
3. `py -X utf8 scripts/update-nearby-parks.py` で取り込む

## 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 倍を外れると警告が出ます。

OSM側にすでにポリゴンがある公園を手入力で置き換える場合（例: 道路で2つに分かれている公園）は、
同じ `our_name` で複数のポリゴンを描いたうえで、`--force-name` と `--override-name` を付けます。

```bash
py -X utf8 scripts/import-manual-parks.py reference/qgis/manual-parks.geojson --force-name いずみ公園 --override-name いずみ公園 --dry-run
py -X utf8 scripts/import-manual-parks.py reference/qgis/manual-parks.geojson --force-name いずみ公園 --override-name いずみ公園
```

`--override-name` を付けた公園は、OSM由来のポリゴンの代わりに手入力ポリゴンだけを地図に表示します。
市外の近隣公園は `py -X utf8 scripts/update-nearby-parks.py` で取り込みます。
取り込み後に `npm run build` で地図と表を確認します。

## 進捗の管理

`parks-without-polygon.geojson` の `drawn` は、`manual-parks.geojson` または `nearby-parks-manual.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']} 件）")
    print(f"市外の近隣公園（点）: {counts['nearby_points']} 件")
    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()