550 lines
20 KiB
Python
550 lines
20 KiB
Python
#!/usr/bin/env python3
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
import datetime as dt
|
||
import glob
|
||
import json
|
||
import math
|
||
import sqlite3
|
||
import zipfile
|
||
from dataclasses import dataclass
|
||
from pathlib import Path
|
||
from typing import Iterable
|
||
import xml.etree.ElementTree as ET
|
||
|
||
from shapely.geometry import LineString, MultiLineString, box, mapping
|
||
from shapely.prepared import prep
|
||
|
||
|
||
RADIUS = 6378137.0
|
||
MAX_MERCATOR_LAT = 85.0511287798066
|
||
DEFAULT_INPUT_GLOB = "coastline/C23-06_*_GML.zip"
|
||
DEFAULT_OUTPUT_DIR = "out/coastline/japan_coast_grid"
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class SourcePackage:
|
||
path: Path
|
||
prefecture: str
|
||
code: str
|
||
|
||
|
||
def mercator_x(lon: float) -> float:
|
||
return RADIUS * math.radians(lon)
|
||
|
||
|
||
def mercator_y(lat: float) -> float:
|
||
lat = max(min(lat, MAX_MERCATOR_LAT), -MAX_MERCATOR_LAT)
|
||
return RADIUS * math.log(math.tan(math.pi / 4.0 + math.radians(lat) / 2.0))
|
||
|
||
|
||
def lon_from_mercator(x: float) -> float:
|
||
return math.degrees(x / RADIUS)
|
||
|
||
|
||
def lat_from_mercator(y: float) -> float:
|
||
return math.degrees(2.0 * math.atan(math.exp(y / RADIUS)) - math.pi / 2.0)
|
||
|
||
|
||
def align_floor(value: float, step: float) -> float:
|
||
return math.floor(value / step) * step
|
||
|
||
|
||
def align_ceil(value: float, step: float) -> float:
|
||
return math.ceil(value / step) * step
|
||
|
||
|
||
def parse_text_list(text: str) -> list[float]:
|
||
return [float(part) for part in text.split() if part]
|
||
|
||
|
||
def detect_prefecture(meta_xml: str) -> str:
|
||
title_start = meta_xml.find("<title>")
|
||
if title_start == -1:
|
||
return "unknown"
|
||
title_end = meta_xml.find("</title>", title_start)
|
||
if title_end == -1:
|
||
return "unknown"
|
||
return meta_xml[title_start + 7:title_end].strip()
|
||
|
||
|
||
def load_package(path: Path) -> SourcePackage:
|
||
with zipfile.ZipFile(path) as zf:
|
||
meta_name = next(name for name in zf.namelist() if "META" in name and name.endswith(".xml"))
|
||
meta_xml = zf.read(meta_name).decode("shift_jis", errors="replace")
|
||
prefecture = detect_prefecture(meta_xml)
|
||
code = path.stem.replace("_GML", "")
|
||
return SourcePackage(path=path, prefecture=prefecture, code=code)
|
||
|
||
|
||
def iter_coastline_lines(zip_path: Path) -> Iterable[LineString]:
|
||
with zipfile.ZipFile(zip_path) as zf:
|
||
xml_name = next(name for name in zf.namelist() if name.endswith(".xml") and "META" not in name)
|
||
root = ET.fromstring(zf.read(xml_name))
|
||
|
||
ns = {"gml": "http://www.opengis.net/gml/3.2"}
|
||
for curve in root.findall(".//gml:Curve", ns):
|
||
coords: list[tuple[float, float]] = []
|
||
for pos_list in curve.findall(".//gml:posList", ns):
|
||
if not pos_list.text:
|
||
continue
|
||
values = parse_text_list(pos_list.text)
|
||
if len(values) < 4 or len(values) % 2 != 0:
|
||
continue
|
||
for i in range(0, len(values), 2):
|
||
lat = values[i]
|
||
lon = values[i + 1]
|
||
coords.append((mercator_x(lon), mercator_y(lat)))
|
||
if len(coords) < 2:
|
||
continue
|
||
yield LineString(coords)
|
||
|
||
|
||
def open_sqlite(path: Path) -> sqlite3.Connection:
|
||
if path.exists():
|
||
path.unlink()
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
conn = sqlite3.connect(path)
|
||
conn.execute("PRAGMA journal_mode=OFF")
|
||
conn.execute("PRAGMA synchronous=OFF")
|
||
conn.execute("PRAGMA temp_store=MEMORY")
|
||
conn.execute("PRAGMA cache_size=200000")
|
||
conn.execute("PRAGMA locking_mode=EXCLUSIVE")
|
||
conn.execute("PRAGMA foreign_keys=OFF")
|
||
return conn
|
||
|
||
|
||
def create_schema(conn: sqlite3.Connection, *, store_coastline_feature: bool) -> None:
|
||
statements = [
|
||
"""
|
||
CREATE TABLE meta (
|
||
key TEXT PRIMARY KEY,
|
||
value TEXT NOT NULL
|
||
)
|
||
""",
|
||
"""
|
||
CREATE TABLE package_stat (
|
||
source_code TEXT NOT NULL,
|
||
prefecture TEXT NOT NULL,
|
||
zip_path TEXT NOT NULL,
|
||
line_count INTEGER NOT NULL,
|
||
point_count INTEGER NOT NULL,
|
||
bbox_minx REAL NOT NULL,
|
||
bbox_miny REAL NOT NULL,
|
||
bbox_maxx REAL NOT NULL,
|
||
bbox_maxy REAL NOT NULL,
|
||
grid_minx REAL NOT NULL,
|
||
grid_miny REAL NOT NULL,
|
||
grid_maxx REAL NOT NULL,
|
||
grid_maxy REAL NOT NULL,
|
||
blocked_count INTEGER NOT NULL,
|
||
candidate_count INTEGER NOT NULL
|
||
)
|
||
""",
|
||
"""
|
||
CREATE TABLE grid_cell (
|
||
cell_id TEXT PRIMARY KEY,
|
||
source_code TEXT NOT NULL,
|
||
prefecture TEXT NOT NULL,
|
||
ix INTEGER NOT NULL,
|
||
iy INTEGER NOT NULL,
|
||
minx REAL NOT NULL,
|
||
miny REAL NOT NULL,
|
||
maxx REAL NOT NULL,
|
||
maxy REAL NOT NULL,
|
||
state INTEGER NOT NULL,
|
||
state_name TEXT NOT NULL,
|
||
source_mask INTEGER NOT NULL,
|
||
blocked_permille INTEGER NOT NULL
|
||
)
|
||
""",
|
||
"CREATE INDEX idx_grid_cell_state ON grid_cell(state)",
|
||
"CREATE INDEX idx_grid_cell_source ON grid_cell(source_code, prefecture)",
|
||
]
|
||
if store_coastline_feature:
|
||
statements.insert(
|
||
2,
|
||
"""
|
||
CREATE TABLE coastline_feature (
|
||
feature_id INTEGER PRIMARY KEY AUTOINCREMENT,
|
||
source_code TEXT NOT NULL,
|
||
prefecture TEXT NOT NULL,
|
||
curve_index INTEGER NOT NULL,
|
||
point_count INTEGER NOT NULL,
|
||
minx REAL NOT NULL,
|
||
miny REAL NOT NULL,
|
||
maxx REAL NOT NULL,
|
||
maxy REAL NOT NULL,
|
||
geom_wkt TEXT NOT NULL
|
||
)
|
||
""",
|
||
)
|
||
conn.executescript(";\n".join(statements))
|
||
|
||
|
||
def upsert_cells(conn: sqlite3.Connection, rows: list[tuple]) -> None:
|
||
if not rows:
|
||
return
|
||
conn.executemany(
|
||
"""
|
||
INSERT INTO grid_cell(
|
||
cell_id, source_code, prefecture, ix, iy, minx, miny, maxx, maxy,
|
||
state, state_name, source_mask, blocked_permille
|
||
) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)
|
||
ON CONFLICT(cell_id) DO UPDATE SET
|
||
source_code = CASE WHEN excluded.state > grid_cell.state THEN excluded.source_code ELSE grid_cell.source_code END,
|
||
prefecture = CASE WHEN excluded.state > grid_cell.state THEN excluded.prefecture ELSE grid_cell.prefecture END,
|
||
state = CASE WHEN excluded.state > grid_cell.state THEN excluded.state ELSE grid_cell.state END,
|
||
state_name = CASE WHEN excluded.state > grid_cell.state THEN excluded.state_name ELSE grid_cell.state_name END,
|
||
source_mask = grid_cell.source_mask | excluded.source_mask,
|
||
blocked_permille = CASE
|
||
WHEN excluded.blocked_permille > grid_cell.blocked_permille THEN excluded.blocked_permille
|
||
ELSE grid_cell.blocked_permille
|
||
END
|
||
""",
|
||
rows,
|
||
)
|
||
|
||
|
||
def write_geojson_header(fh) -> None:
|
||
fh.write('{"type":"FeatureCollection","features":[\n')
|
||
|
||
|
||
def write_geojson_feature(fh, feature: dict, first: bool) -> bool:
|
||
if not first:
|
||
fh.write(",\n")
|
||
fh.write(json.dumps(feature, ensure_ascii=False))
|
||
return False
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description="Build Japan coastline 200m grid into SQLite")
|
||
parser.add_argument(
|
||
"--input",
|
||
dest="inputs",
|
||
action="append",
|
||
default=[],
|
||
help="海岸线 GML zip。可重复指定；默认自动扫描 coastline/C23-06_*_GML.zip",
|
||
)
|
||
parser.add_argument(
|
||
"--output-dir",
|
||
default=DEFAULT_OUTPUT_DIR,
|
||
help="输出目录",
|
||
)
|
||
parser.add_argument("--cell-size-m", type=float, default=200.0, help="格网尺寸，默认 200m")
|
||
parser.add_argument("--buffer-m", type=float, default=50.0, help="海岸线缓冲距离，默认 50m")
|
||
parser.add_argument(
|
||
"--margin-cells",
|
||
type=int,
|
||
default=2,
|
||
help="输出范围外扩的格子数，默认 2 个格子",
|
||
)
|
||
parser.add_argument(
|
||
"--write-geojson",
|
||
action="store_true",
|
||
help="同时输出 coastline.geojson 和 grid.geojson",
|
||
)
|
||
parser.add_argument(
|
||
"--store-coastline-feature",
|
||
action="store_true",
|
||
help="把每条海岸线曲线也写入 SQLite，默认关闭以控制库体积",
|
||
)
|
||
args = parser.parse_args()
|
||
|
||
project_root = Path(__file__).resolve().parent.parent
|
||
if args.inputs:
|
||
input_paths: list[Path] = []
|
||
for item in args.inputs:
|
||
if any(ch in item for ch in "*?[]"):
|
||
input_paths.extend(Path(path) for path in sorted(glob.glob(item)))
|
||
else:
|
||
input_paths.append(Path(item))
|
||
else:
|
||
input_paths = [Path(path) for path in sorted(glob.glob(str(project_root / DEFAULT_INPUT_GLOB)))]
|
||
|
||
if not input_paths:
|
||
raise SystemExit("no coastline packages found")
|
||
|
||
packages: list[SourcePackage] = []
|
||
for path in input_paths:
|
||
if not path.is_absolute():
|
||
path = project_root / path
|
||
if not path.exists():
|
||
raise SystemExit(f"missing input: {path}")
|
||
packages.append(load_package(path))
|
||
|
||
out_dir = Path(args.output_dir)
|
||
if not out_dir.is_absolute():
|
||
out_dir = project_root / out_dir
|
||
out_dir.mkdir(parents=True, exist_ok=True)
|
||
|
||
sqlite_path = out_dir / "japan_coast_grid.sqlite"
|
||
conn = open_sqlite(sqlite_path)
|
||
create_schema(conn, store_coastline_feature=args.store_coastline_feature)
|
||
|
||
coastline_geojson_fh = None
|
||
grid_geojson_fh = None
|
||
coastline_geojson_first = True
|
||
grid_geojson_first = True
|
||
if args.write_geojson:
|
||
coastline_geojson_fh = (out_dir / "japan_coastline.geojson").open("w", encoding="utf-8")
|
||
grid_geojson_fh = (out_dir / "japan_grid.geojson").open("w", encoding="utf-8")
|
||
write_geojson_header(coastline_geojson_fh)
|
||
write_geojson_header(grid_geojson_fh)
|
||
|
||
package_stats: list[dict] = []
|
||
total_lines = 0
|
||
total_points = 0
|
||
|
||
try:
|
||
for package in packages:
|
||
print(f"[{package.code}] 开始处理 {package.path.name}")
|
||
lines = list(iter_coastline_lines(package.path))
|
||
if not lines:
|
||
print(f"[{package.code}] 没有可用海岸线，跳过")
|
||
continue
|
||
|
||
coastline_geom = MultiLineString([line.coords[:] for line in lines])
|
||
coastline_buffer = coastline_geom.buffer(args.buffer_m, cap_style=2, join_style=2)
|
||
prepared_buffer = prep(coastline_buffer)
|
||
coast_bounds = coastline_geom.bounds
|
||
|
||
minx = align_floor(coast_bounds[0] - args.buffer_m - args.margin_cells * args.cell_size_m, args.cell_size_m)
|
||
miny = align_floor(coast_bounds[1] - args.buffer_m - args.margin_cells * args.cell_size_m, args.cell_size_m)
|
||
maxx = align_ceil(coast_bounds[2] + args.buffer_m + args.margin_cells * args.cell_size_m, args.cell_size_m)
|
||
maxy = align_ceil(coast_bounds[3] + args.buffer_m + args.margin_cells * args.cell_size_m, args.cell_size_m)
|
||
|
||
line_count = len(lines)
|
||
point_count = sum(len(line.coords) for line in lines)
|
||
total_lines += line_count
|
||
total_points += point_count
|
||
|
||
if args.store_coastline_feature:
|
||
conn.executemany(
|
||
"""
|
||
INSERT INTO coastline_feature(
|
||
source_code, prefecture, curve_index, point_count, minx, miny, maxx, maxy, geom_wkt
|
||
) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?)
|
||
""",
|
||
[
|
||
(
|
||
package.code,
|
||
package.prefecture,
|
||
idx,
|
||
len(line.coords),
|
||
line.bounds[0],
|
||
line.bounds[1],
|
||
line.bounds[2],
|
||
line.bounds[3],
|
||
line.wkt,
|
||
)
|
||
for idx, line in enumerate(lines, start=1)
|
||
],
|
||
)
|
||
|
||
package_blocked = 0
|
||
package_candidate = 0
|
||
package_rows: list[tuple] = []
|
||
|
||
ix0 = int(round(minx / args.cell_size_m))
|
||
iy0 = int(round(miny / args.cell_size_m))
|
||
ix1 = int(round(maxx / args.cell_size_m))
|
||
iy1 = int(round(maxy / args.cell_size_m))
|
||
|
||
for ix in range(ix0, ix1):
|
||
cell_minx = ix * args.cell_size_m
|
||
cell_maxx = cell_minx + args.cell_size_m
|
||
for iy in range(iy0, iy1):
|
||
cell_miny = iy * args.cell_size_m
|
||
cell_maxy = cell_miny + args.cell_size_m
|
||
cell = box(cell_minx, cell_miny, cell_maxx, cell_maxy)
|
||
intersects = prepared_buffer.intersects(cell)
|
||
state = 2 if intersects else 1
|
||
state_name = "HARD_BLOCKED" if intersects else "NAVIGABLE_CANDIDATE"
|
||
blocked_permille = 1000 if intersects else 0
|
||
cell_id = f"{ix}:{iy}"
|
||
|
||
package_rows.append(
|
||
(
|
||
cell_id,
|
||
package.code,
|
||
package.prefecture,
|
||
ix,
|
||
iy,
|
||
cell_minx,
|
||
cell_miny,
|
||
cell_maxx,
|
||
cell_maxy,
|
||
state,
|
||
state_name,
|
||
1,
|
||
blocked_permille,
|
||
)
|
||
)
|
||
if intersects:
|
||
package_blocked += 1
|
||
else:
|
||
package_candidate += 1
|
||
|
||
if args.write_geojson and grid_geojson_fh is not None:
|
||
feature = {
|
||
"type": "Feature",
|
||
"properties": {
|
||
"cell_id": cell_id,
|
||
"source_code": package.code,
|
||
"prefecture": package.prefecture,
|
||
"ix": ix,
|
||
"iy": iy,
|
||
"cell_size_m": args.cell_size_m,
|
||
"state": state,
|
||
"state_name": state_name,
|
||
"source_mask": 1,
|
||
"blocked_permille": blocked_permille,
|
||
},
|
||
"geometry": mapping(cell),
|
||
}
|
||
grid_geojson_first = write_geojson_feature(grid_geojson_fh, feature, grid_geojson_first)
|
||
|
||
if len(package_rows) >= 5000:
|
||
upsert_cells(conn, package_rows)
|
||
conn.commit()
|
||
print(
|
||
f"[{package.code}] 已写入 {len(package_rows)} 行，"
|
||
f"blocked={package_blocked} candidate={package_candidate}"
|
||
)
|
||
package_rows.clear()
|
||
|
||
if package_rows:
|
||
upsert_cells(conn, package_rows)
|
||
conn.commit()
|
||
|
||
package_stats.append(
|
||
{
|
||
"source_code": package.code,
|
||
"prefecture": package.prefecture,
|
||
"zip_path": str(package.path),
|
||
"line_count": line_count,
|
||
"point_count": point_count,
|
||
"bbox_mercator": [coast_bounds[0], coast_bounds[1], coast_bounds[2], coast_bounds[3]],
|
||
"grid_bounds_mercator": [minx, miny, maxx, maxy],
|
||
"blocked_count": package_blocked,
|
||
"candidate_count": package_candidate,
|
||
}
|
||
)
|
||
|
||
conn.execute(
|
||
"""
|
||
INSERT INTO package_stat(
|
||
source_code, prefecture, zip_path, line_count, point_count,
|
||
bbox_minx, bbox_miny, bbox_maxx, bbox_maxy,
|
||
grid_minx, grid_miny, grid_maxx, grid_maxy,
|
||
blocked_count, candidate_count
|
||
) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)
|
||
""",
|
||
(
|
||
package.code,
|
||
package.prefecture,
|
||
str(package.path),
|
||
line_count,
|
||
point_count,
|
||
coast_bounds[0],
|
||
coast_bounds[1],
|
||
coast_bounds[2],
|
||
coast_bounds[3],
|
||
minx,
|
||
miny,
|
||
maxx,
|
||
maxy,
|
||
package_blocked,
|
||
package_candidate,
|
||
),
|
||
)
|
||
conn.commit()
|
||
print(
|
||
f"[{package.code}] 完成 line={line_count} point={point_count} "
|
||
f"blocked={package_blocked} candidate={package_candidate}"
|
||
)
|
||
|
||
conn.execute("ANALYZE")
|
||
conn.commit()
|
||
|
||
final_total = int(conn.execute("SELECT COUNT(*) FROM grid_cell").fetchone()[0])
|
||
blocked_count = int(conn.execute("SELECT COUNT(*) FROM grid_cell WHERE state = 2").fetchone()[0])
|
||
candidate_count = int(conn.execute("SELECT COUNT(*) FROM grid_cell WHERE state = 1").fetchone()[0])
|
||
bbox_row = conn.execute(
|
||
"SELECT MIN(minx), MIN(miny), MAX(maxx), MAX(maxy) FROM grid_cell"
|
||
).fetchone()
|
||
grid_bbox = (
|
||
float(bbox_row[0]) if bbox_row and bbox_row[0] is not None else float("inf"),
|
||
float(bbox_row[1]) if bbox_row and bbox_row[1] is not None else float("inf"),
|
||
float(bbox_row[2]) if bbox_row and bbox_row[2] is not None else float("-inf"),
|
||
float(bbox_row[3]) if bbox_row and bbox_row[3] is not None else float("-inf"),
|
||
)
|
||
|
||
build_time = dt.datetime.now().isoformat(timespec="seconds")
|
||
meta_rows = [
|
||
("build_time", build_time),
|
||
("cell_size_m", str(args.cell_size_m)),
|
||
("buffer_m", str(args.buffer_m)),
|
||
("margin_cells", str(args.margin_cells)),
|
||
("source_scope", "japan"),
|
||
("source_code", "C23-06_*"),
|
||
("source_count", str(len(packages))),
|
||
("source_packages", json.dumps([str(pkg.path) for pkg in packages], ensure_ascii=False)),
|
||
("coastline_line_count", str(total_lines)),
|
||
("coastline_point_count", str(total_points)),
|
||
("grid_cell_count", str(final_total)),
|
||
("blocked_count", str(blocked_count)),
|
||
("candidate_count", str(candidate_count)),
|
||
("grid_bounds_mercator", json.dumps(grid_bbox, ensure_ascii=False)),
|
||
("notes", "全国海岸线 200m 硬阻塞栅格，来源为 coastline/C23-06_*_GML.zip"),
|
||
]
|
||
conn.executemany("INSERT INTO meta(key, value) VALUES (?, ?)", meta_rows)
|
||
conn.commit()
|
||
|
||
meta_path = out_dir / "japan_coast_grid.meta.json"
|
||
meta_payload = {
|
||
"build_time": build_time,
|
||
"cell_size_m": args.cell_size_m,
|
||
"buffer_m": args.buffer_m,
|
||
"margin_cells": args.margin_cells,
|
||
"source_scope": "japan",
|
||
"source_count": len(packages),
|
||
"source_packages": [str(pkg.path) for pkg in packages],
|
||
"coastline_line_count": total_lines,
|
||
"coastline_point_count": total_points,
|
||
"grid_cell_count": final_total,
|
||
"blocked_count": blocked_count,
|
||
"candidate_count": candidate_count,
|
||
"grid_bounds_mercator": list(grid_bbox),
|
||
"package_stats": package_stats,
|
||
}
|
||
meta_path.write_text(json.dumps(meta_payload, ensure_ascii=False, indent=2), encoding="utf-8")
|
||
|
||
if args.write_geojson and coastline_geojson_fh is not None and grid_geojson_fh is not None:
|
||
coastline_geojson_fh.write("\n]}\n")
|
||
grid_geojson_fh.write("\n]}\n")
|
||
|
||
print("完成")
|
||
print(f" SQLite: {sqlite_path}")
|
||
print(f" 元数据: {meta_path}")
|
||
print(f" packages={len(packages)}")
|
||
print(f" coastline_lines={total_lines}")
|
||
print(f" grid_cells={final_total}")
|
||
print(f" blocked={blocked_count}")
|
||
print(f" candidate={candidate_count}")
|
||
if args.write_geojson:
|
||
print(f" 海岸线GeoJSON: {out_dir / 'japan_coastline.geojson'}")
|
||
print(f" 栅格GeoJSON: {out_dir / 'japan_grid.geojson'}")
|
||
finally:
|
||
if coastline_geojson_fh is not None:
|
||
coastline_geojson_fh.close()
|
||
if grid_geojson_fh is not None:
|
||
grid_geojson_fh.close()
|
||
conn.close()
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|