566 lines
20 KiB
Python
566 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 zipfile
|
||
from dataclasses import dataclass
|
||
from pathlib import Path
|
||
from typing import Iterable
|
||
import xml.etree.ElementTree as ET
|
||
|
||
import pymysql
|
||
from shapely.geometry import LineString, MultiLineString, box, mapping
|
||
from shapely.prepared import prep
|
||
|
||
|
||
DB_NAME = "navsea_japan_coast_grid"
|
||
DB_USER = "root"
|
||
DB_PASSWORD = "2chi9ks2"
|
||
DB_HOST = "localhost"
|
||
DB_SOCKET = "/tmp/mysql.sock"
|
||
|
||
DEFAULT_INPUT_GLOB = "coastline/C23-06_*_GML.zip"
|
||
DEFAULT_OUT_DIR = "out/coastline/japan_coast_grid_mysql"
|
||
LAYER_NAME = "coast_200m"
|
||
CELL_SIZE_M = 200.0
|
||
|
||
RADIUS = 6378137.0
|
||
MAX_MERCATOR_LAT = 85.0511287798066
|
||
NS = {"gml": "http://www.opengis.net/gml/3.2"}
|
||
|
||
|
||
@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 mercator_bbox_to_lonlat(minx: float, miny: float, maxx: float, maxy: float) -> tuple[float, float, float, float]:
|
||
return (
|
||
lon_from_mercator(minx),
|
||
lat_from_mercator(miny),
|
||
lon_from_mercator(maxx),
|
||
lat_from_mercator(maxy),
|
||
)
|
||
|
||
|
||
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))
|
||
|
||
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:
|
||
yield LineString(coords)
|
||
|
||
|
||
def mysql_connect(database: str | None = None):
|
||
kwargs = {
|
||
"host": DB_HOST,
|
||
"user": DB_USER,
|
||
"password": DB_PASSWORD,
|
||
"charset": "utf8mb4",
|
||
"autocommit": False,
|
||
"cursorclass": pymysql.cursors.Cursor,
|
||
}
|
||
if DB_SOCKET:
|
||
kwargs["unix_socket"] = DB_SOCKET
|
||
if database:
|
||
kwargs["database"] = database
|
||
return pymysql.connect(**kwargs)
|
||
|
||
|
||
def ensure_database() -> None:
|
||
conn = mysql_connect()
|
||
try:
|
||
with conn.cursor() as cur:
|
||
cur.execute(
|
||
f"CREATE DATABASE IF NOT EXISTS `{DB_NAME}` "
|
||
"DEFAULT CHARACTER SET utf8mb4 COLLATE utf8mb4_unicode_ci"
|
||
)
|
||
conn.commit()
|
||
finally:
|
||
conn.close()
|
||
|
||
|
||
def ensure_schema(conn) -> None:
|
||
ddl = [
|
||
"DROP TABLE IF EXISTS navsea_grid_cell",
|
||
"DROP TABLE IF EXISTS navsea_grid_package_stat",
|
||
"DROP TABLE IF EXISTS navsea_grid_layer_meta",
|
||
"""
|
||
CREATE TABLE navsea_grid_layer_meta (
|
||
layer_name VARCHAR(32) NOT NULL PRIMARY KEY,
|
||
description VARCHAR(255) NOT NULL,
|
||
source_desc TEXT NOT NULL,
|
||
cell_size_m DOUBLE NOT NULL,
|
||
feature_count BIGINT NOT NULL,
|
||
bbox_min_lon DOUBLE NOT NULL,
|
||
bbox_min_lat DOUBLE NOT NULL,
|
||
bbox_max_lon DOUBLE NOT NULL,
|
||
bbox_max_lat DOUBLE NOT NULL,
|
||
export_file VARCHAR(255) DEFAULT NULL,
|
||
created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP
|
||
) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4
|
||
""",
|
||
"""
|
||
CREATE TABLE navsea_grid_package_stat (
|
||
source_code VARCHAR(64) NOT NULL,
|
||
prefecture VARCHAR(128) NOT NULL,
|
||
zip_path TEXT NOT NULL,
|
||
line_count BIGINT NOT NULL,
|
||
point_count BIGINT NOT NULL,
|
||
bbox_min_lon DOUBLE NOT NULL,
|
||
bbox_min_lat DOUBLE NOT NULL,
|
||
bbox_max_lon DOUBLE NOT NULL,
|
||
bbox_max_lat DOUBLE NOT NULL,
|
||
grid_min_lon DOUBLE NOT NULL,
|
||
grid_min_lat DOUBLE NOT NULL,
|
||
grid_max_lon DOUBLE NOT NULL,
|
||
grid_max_lat DOUBLE NOT NULL,
|
||
blocked_count BIGINT NOT NULL,
|
||
candidate_count BIGINT NOT NULL
|
||
) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4
|
||
""",
|
||
"""
|
||
CREATE TABLE navsea_grid_cell (
|
||
id BIGINT UNSIGNED NOT NULL AUTO_INCREMENT PRIMARY KEY,
|
||
layer_name VARCHAR(32) NOT NULL,
|
||
cell_id VARCHAR(64) NOT NULL,
|
||
row_idx INT NOT NULL,
|
||
col_idx INT NOT NULL,
|
||
cell_size_m DOUBLE NOT NULL,
|
||
state_name VARCHAR(32) NOT NULL,
|
||
class_name VARCHAR(32) NOT NULL,
|
||
source_name VARCHAR(191) NOT NULL,
|
||
min_lon DOUBLE NOT NULL,
|
||
min_lat DOUBLE NOT NULL,
|
||
max_lon DOUBLE NOT NULL,
|
||
max_lat DOUBLE NOT NULL,
|
||
created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP,
|
||
UNIQUE KEY uniq_layer_cell (layer_name, cell_id),
|
||
KEY idx_layer_state (layer_name, state_name),
|
||
KEY idx_layer_rowcol (layer_name, row_idx, col_idx)
|
||
) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4
|
||
""",
|
||
]
|
||
with conn.cursor() as cur:
|
||
for stmt in ddl:
|
||
cur.execute(stmt)
|
||
conn.commit()
|
||
|
||
|
||
def insert_cells(cur, rows: list[tuple]) -> None:
|
||
cur.executemany(
|
||
"""
|
||
INSERT INTO navsea_grid_cell
|
||
(layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, source_name,
|
||
min_lon, min_lat, max_lon, max_lat)
|
||
VALUES
|
||
(%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s)
|
||
ON DUPLICATE KEY UPDATE
|
||
row_idx = row_idx,
|
||
col_idx = col_idx,
|
||
cell_size_m = cell_size_m,
|
||
state_name = IF(VALUES(state_name) = 'HARD_BLOCKED' AND state_name <> 'HARD_BLOCKED', VALUES(state_name), state_name),
|
||
class_name = IF(VALUES(state_name) = 'HARD_BLOCKED' AND state_name <> 'HARD_BLOCKED', VALUES(class_name), class_name),
|
||
source_name = IF(VALUES(state_name) = 'HARD_BLOCKED' AND state_name <> 'HARD_BLOCKED', VALUES(source_name), source_name),
|
||
min_lon = min_lon,
|
||
min_lat = min_lat,
|
||
max_lon = max_lon,
|
||
max_lat = max_lat
|
||
""",
|
||
rows,
|
||
)
|
||
|
||
|
||
def upsert_meta(
|
||
cur,
|
||
*,
|
||
layer_name: str,
|
||
description: str,
|
||
source_desc: str,
|
||
cell_size_m: float,
|
||
feature_count: int,
|
||
bbox: tuple[float, float, float, float],
|
||
export_file: str | None,
|
||
) -> None:
|
||
cur.execute(
|
||
"""
|
||
INSERT INTO navsea_grid_layer_meta
|
||
(layer_name, description, source_desc, cell_size_m, feature_count,
|
||
bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat, export_file)
|
||
VALUES
|
||
(%s, %s, %s, %s, %s, %s, %s, %s, %s, %s)
|
||
ON DUPLICATE KEY UPDATE
|
||
description=VALUES(description),
|
||
source_desc=VALUES(source_desc),
|
||
cell_size_m=VALUES(cell_size_m),
|
||
feature_count=VALUES(feature_count),
|
||
bbox_min_lon=VALUES(bbox_min_lon),
|
||
bbox_min_lat=VALUES(bbox_min_lat),
|
||
bbox_max_lon=VALUES(bbox_max_lon),
|
||
bbox_max_lat=VALUES(bbox_max_lat),
|
||
export_file=VALUES(export_file)
|
||
""",
|
||
(
|
||
layer_name,
|
||
description,
|
||
source_desc,
|
||
cell_size_m,
|
||
feature_count,
|
||
bbox[0],
|
||
bbox[1],
|
||
bbox[2],
|
||
bbox[3],
|
||
export_file,
|
||
),
|
||
)
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description="Build Japan coastline 200m grid directly into MySQL")
|
||
parser.add_argument(
|
||
"--input",
|
||
dest="inputs",
|
||
action="append",
|
||
default=[],
|
||
help="海岸线 GML zip。可重复指定；默认自动扫描 coastline/C23-06_*_GML.zip",
|
||
)
|
||
parser.add_argument(
|
||
"--db-name",
|
||
default=DB_NAME,
|
||
help="MySQL 数据库名",
|
||
)
|
||
parser.add_argument(
|
||
"--out-dir",
|
||
default=DEFAULT_OUT_DIR,
|
||
help="输出目录，用于可选 GeoJSON 和元数据",
|
||
)
|
||
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",
|
||
)
|
||
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.out_dir)
|
||
if not out_dir.is_absolute():
|
||
out_dir = project_root / out_dir
|
||
out_dir.mkdir(parents=True, exist_ok=True)
|
||
|
||
ensure_database()
|
||
conn = mysql_connect(args.db_name)
|
||
try:
|
||
print(f"数据库 {args.db_name} 已连接，准备初始化表结构 ...")
|
||
ensure_schema(conn)
|
||
|
||
total_lines = 0
|
||
total_points = 0
|
||
package_stats: list[dict] = []
|
||
|
||
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([list(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 * CELL_SIZE_M, CELL_SIZE_M)
|
||
miny = align_floor(coast_bounds[1] - args.buffer_m - args.margin_cells * CELL_SIZE_M, CELL_SIZE_M)
|
||
maxx = align_ceil(coast_bounds[2] + args.buffer_m + args.margin_cells * CELL_SIZE_M, CELL_SIZE_M)
|
||
maxy = align_ceil(coast_bounds[3] + args.buffer_m + args.margin_cells * CELL_SIZE_M, 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
|
||
|
||
package_blocked = 0
|
||
package_candidate = 0
|
||
package_rows: list[tuple] = []
|
||
|
||
ix0 = int(round(minx / CELL_SIZE_M))
|
||
iy0 = int(round(miny / CELL_SIZE_M))
|
||
ix1 = int(round(maxx / CELL_SIZE_M))
|
||
iy1 = int(round(maxy / CELL_SIZE_M))
|
||
|
||
for ix in range(ix0, ix1):
|
||
cell_minx = ix * CELL_SIZE_M
|
||
cell_maxx = cell_minx + CELL_SIZE_M
|
||
for iy in range(iy0, iy1):
|
||
cell_miny = iy * CELL_SIZE_M
|
||
cell_maxy = cell_miny + CELL_SIZE_M
|
||
cell = box(cell_minx, cell_miny, cell_maxx, cell_maxy)
|
||
intersects = prepared_buffer.intersects(cell)
|
||
if intersects:
|
||
cell_id = f"{ix}:{iy}"
|
||
package_rows.append(
|
||
(
|
||
LAYER_NAME,
|
||
cell_id,
|
||
iy,
|
||
ix,
|
||
CELL_SIZE_M,
|
||
"HARD_BLOCKED",
|
||
"HARD_BLOCKED",
|
||
package.code,
|
||
cell_minx,
|
||
cell_miny,
|
||
cell_maxx,
|
||
cell_maxy,
|
||
)
|
||
)
|
||
package_blocked += 1
|
||
else:
|
||
package_candidate += 1
|
||
|
||
if len(package_rows) >= 5000:
|
||
with conn.cursor() as cur:
|
||
insert_cells(cur, package_rows)
|
||
conn.commit()
|
||
print(
|
||
f"[{package.code}] 已写入 {len(package_rows)} 行，"
|
||
f"blocked={package_blocked} candidate={package_candidate}"
|
||
)
|
||
package_rows.clear()
|
||
|
||
if package_rows:
|
||
with conn.cursor() as cur:
|
||
insert_cells(cur, 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,
|
||
}
|
||
)
|
||
|
||
with conn.cursor() as cur:
|
||
cur.execute(
|
||
"""
|
||
INSERT INTO navsea_grid_package_stat(
|
||
source_code, prefecture, zip_path, line_count, point_count,
|
||
bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat,
|
||
grid_min_lon, grid_min_lat, grid_max_lon, grid_max_lat,
|
||
blocked_count, candidate_count
|
||
) VALUES (
|
||
%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s
|
||
)
|
||
""",
|
||
(
|
||
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}"
|
||
)
|
||
|
||
with conn.cursor() as cur:
|
||
cur.execute("ANALYZE TABLE navsea_grid_cell")
|
||
cur.execute("ANALYZE TABLE navsea_grid_package_stat")
|
||
cur.execute("ANALYZE TABLE navsea_grid_layer_meta")
|
||
conn.commit()
|
||
|
||
with conn.cursor() as cur:
|
||
cur.execute(f"SELECT COUNT(*) FROM navsea_grid_cell WHERE layer_name=%s", (LAYER_NAME,))
|
||
final_total = int(cur.fetchone()[0] or 0)
|
||
cur.execute(
|
||
f"SELECT COUNT(*) FROM navsea_grid_cell WHERE layer_name=%s AND state_name=%s",
|
||
(LAYER_NAME, "HARD_BLOCKED"),
|
||
)
|
||
blocked_count = int(cur.fetchone()[0] or 0)
|
||
cur.execute(
|
||
f"SELECT COUNT(*) FROM navsea_grid_cell WHERE layer_name=%s AND state_name=%s",
|
||
(LAYER_NAME, "NAVIGABLE_CANDIDATE"),
|
||
)
|
||
candidate_count = int(cur.fetchone()[0] or 0)
|
||
cur.execute(
|
||
"""
|
||
SELECT MIN(min_lon), MIN(min_lat), MAX(max_lon), MAX(max_lat)
|
||
FROM navsea_grid_cell
|
||
WHERE layer_name=%s
|
||
""",
|
||
(LAYER_NAME,),
|
||
)
|
||
bbox_row = cur.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")
|
||
grid_bbox_lonlat = mercator_bbox_to_lonlat(*grid_bbox)
|
||
|
||
manifest = {
|
||
"build_time": build_time,
|
||
"db_name": args.db_name,
|
||
"source_scope": "japan",
|
||
"source_code": "C23-06_*",
|
||
"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_lonlat": list(grid_bbox_lonlat),
|
||
"package_stats": package_stats,
|
||
"notes": "全国海岸线 200m 硬阻塞栅格，直接写入 MySQL",
|
||
}
|
||
manifest_path = out_dir / "japan_coast_grid_mysql.manifest.json"
|
||
manifest_path.write_text(json.dumps(manifest, ensure_ascii=False, indent=2), encoding="utf-8")
|
||
|
||
with conn.cursor() as cur:
|
||
upsert_meta(
|
||
cur,
|
||
layer_name=LAYER_NAME,
|
||
description="全国海岸线 200m 硬阻塞格",
|
||
source_desc="coastline/C23-06_*_GML.zip",
|
||
cell_size_m=CELL_SIZE_M,
|
||
feature_count=final_total,
|
||
bbox=grid_bbox_lonlat,
|
||
export_file=str(manifest_path.relative_to(project_root)),
|
||
)
|
||
conn.commit()
|
||
|
||
print("完成")
|
||
print(f" MySQL: {args.db_name}")
|
||
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}")
|
||
print(f" manifest: {manifest_path}")
|
||
finally:
|
||
conn.close()
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|