871 lines
34 KiB
Python
871 lines
34 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_v2"
|
||
DEFAULT_CELL_SIZE_M = 200.0
|
||
DEFAULT_BUFFER_M = 50.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(name: str) -> None:
|
||
conn = mysql_connect()
|
||
try:
|
||
with conn.cursor() as cur:
|
||
cur.execute(
|
||
f"CREATE DATABASE IF NOT EXISTS `{name}` "
|
||
"DEFAULT CHARACTER SET utf8mb4 COLLATE utf8mb4_unicode_ci"
|
||
)
|
||
conn.commit()
|
||
finally:
|
||
conn.close()
|
||
|
||
|
||
def column_exists(conn, table_name: str, column_name: str) -> bool:
|
||
with conn.cursor() as cur:
|
||
cur.execute("SELECT DATABASE()")
|
||
database_name = cur.fetchone()[0]
|
||
cur.execute(
|
||
"""
|
||
SELECT COUNT(*)
|
||
FROM information_schema.COLUMNS
|
||
WHERE TABLE_SCHEMA=%s AND TABLE_NAME=%s AND COLUMN_NAME=%s
|
||
""",
|
||
(database_name, table_name, column_name),
|
||
)
|
||
return int(cur.fetchone()[0] or 0) > 0
|
||
|
||
|
||
def ensure_schema(conn) -> None:
|
||
ddl = [
|
||
"""
|
||
CREATE TABLE IF NOT EXISTS navsea_coast_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 IF NOT EXISTS navsea_coast_grid_package_stat (
|
||
layer_name VARCHAR(32) NOT NULL,
|
||
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,
|
||
KEY idx_layer_source (layer_name, source_code)
|
||
) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4
|
||
""",
|
||
"""
|
||
CREATE TABLE IF NOT EXISTS navsea_coast_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_code VARCHAR(64) NOT NULL,
|
||
source_path TEXT 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_source_cell (layer_name, source_code, 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)
|
||
if not column_exists(conn, "navsea_coast_grid_package_stat", "layer_name"):
|
||
with conn.cursor() as cur:
|
||
cur.execute(
|
||
"ALTER TABLE navsea_coast_grid_package_stat "
|
||
"ADD COLUMN layer_name VARCHAR(32) NOT NULL DEFAULT 'coast_200m' FIRST"
|
||
)
|
||
cur.execute("ALTER TABLE navsea_coast_grid_package_stat ADD KEY idx_layer_source (layer_name, source_code)")
|
||
conn.commit()
|
||
|
||
|
||
def replace_layer_rows(conn, layer_name: str) -> None:
|
||
with conn.cursor() as cur:
|
||
cur.execute("DELETE FROM navsea_coast_grid_cell WHERE layer_name=%s", (layer_name,))
|
||
cur.execute("DELETE FROM navsea_coast_grid_layer_meta WHERE layer_name=%s", (layer_name,))
|
||
cur.execute("DELETE FROM navsea_coast_grid_package_stat WHERE layer_name=%s", (layer_name,))
|
||
conn.commit()
|
||
|
||
|
||
def insert_cells(cur, rows: list[tuple]) -> None:
|
||
cur.executemany(
|
||
"""
|
||
INSERT INTO navsea_coast_grid_cell
|
||
(layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, source_code,
|
||
source_path, source_name, min_lon, min_lat, max_lon, max_lat)
|
||
VALUES
|
||
(%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s)
|
||
ON DUPLICATE KEY UPDATE
|
||
row_idx=VALUES(row_idx),
|
||
col_idx=VALUES(col_idx),
|
||
cell_size_m=VALUES(cell_size_m),
|
||
state_name=VALUES(state_name),
|
||
class_name=VALUES(class_name),
|
||
source_code=VALUES(source_code),
|
||
source_path=VALUES(source_path),
|
||
source_name=VALUES(source_name),
|
||
min_lon=VALUES(min_lon),
|
||
min_lat=VALUES(min_lat),
|
||
max_lon=VALUES(max_lon),
|
||
max_lat=VALUES(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_coast_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 default_layer_name(cell_size_m: float) -> str:
|
||
if float(cell_size_m).is_integer():
|
||
return f"coast_{int(cell_size_m)}m"
|
||
return f"coast_{str(cell_size_m).replace('.', '_')}m"
|
||
|
||
|
||
def normalize_source_code(value: str) -> str:
|
||
value = value.strip()
|
||
if value.isdigit():
|
||
return f"C23-06_{int(value):02d}"
|
||
if value.startswith("C23-06_") and value[len("C23-06_"):].isdigit():
|
||
return f"C23-06_{int(value[len('C23-06_'):]):02d}"
|
||
return value
|
||
|
||
|
||
def ts_now() -> str:
|
||
return dt.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
|
||
|
||
|
||
def log(message: str) -> None:
|
||
print(f"[{ts_now()}] {message}", flush=True)
|
||
|
||
|
||
def flush_package_rows(
|
||
conn,
|
||
package_rows: list[tuple],
|
||
*,
|
||
package_code: str,
|
||
package_blocked: int,
|
||
package_candidate: int,
|
||
reason: str,
|
||
) -> None:
|
||
if not package_rows:
|
||
return
|
||
with conn.cursor() as cur:
|
||
insert_cells(cur, package_rows)
|
||
conn.commit()
|
||
log(
|
||
f"[{package_code}] 已写入 {len(package_rows)} 行（{reason}），"
|
||
f"blocked={package_blocked} candidate={package_candidate}"
|
||
)
|
||
package_rows.clear()
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description="Build Japan coastline 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(
|
||
"--cell-size-m",
|
||
type=float,
|
||
default=DEFAULT_CELL_SIZE_M,
|
||
help="海岸线格子大小，默认 200m；例如 50 会自动写入 coast_50m",
|
||
)
|
||
parser.add_argument(
|
||
"--layer-name",
|
||
default=None,
|
||
help="目标 layer_name；默认按 cell size 自动生成，例如 coast_200m / coast_50m",
|
||
)
|
||
parser.add_argument("--buffer-m", type=float, default=DEFAULT_BUFFER_M, 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(
|
||
"--resume",
|
||
action="store_true",
|
||
help="续跑模式：不清理当前 layer，保留已入库格子",
|
||
)
|
||
parser.add_argument(
|
||
"--start-code",
|
||
default=None,
|
||
help="只处理指定 source_code 及之后的包，例如 13 或 C23-06_13",
|
||
)
|
||
parser.add_argument(
|
||
"--progress-every-cells",
|
||
type=int,
|
||
default=50000,
|
||
help="每扫描多少个候选格输出一次进度日志，默认 50000",
|
||
)
|
||
parser.add_argument(
|
||
"--progress-every-seconds",
|
||
type=int,
|
||
default=15,
|
||
help="每隔多少秒至少输出一次进度日志，默认 15",
|
||
)
|
||
parser.add_argument(
|
||
"--coarse-factor",
|
||
type=int,
|
||
default=8,
|
||
help="海岸线粗筛倍率，默认 8；值越大越少细扫，但粗筛更大",
|
||
)
|
||
args = parser.parse_args()
|
||
if args.cell_size_m <= 0:
|
||
raise SystemExit("--cell-size-m must be positive")
|
||
cell_size_m = float(args.cell_size_m)
|
||
layer_name = args.layer_name or default_layer_name(cell_size_m)
|
||
|
||
project_root = Path(__file__).resolve().parent.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")
|
||
if args.progress_every_cells <= 0:
|
||
raise SystemExit("--progress-every-cells must be positive")
|
||
if args.progress_every_seconds < 0:
|
||
raise SystemExit("--progress-every-seconds must be non-negative")
|
||
if args.coarse_factor < 1:
|
||
raise SystemExit("--coarse-factor must be >= 1")
|
||
|
||
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))
|
||
packages.sort(key=lambda package: package.code)
|
||
|
||
if args.start_code:
|
||
start_code = normalize_source_code(args.start_code)
|
||
packages = [package for package in packages if package.code >= start_code]
|
||
if not packages:
|
||
raise SystemExit(f"no coastline packages found from start code: {start_code}")
|
||
|
||
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(args.db_name)
|
||
conn = mysql_connect(args.db_name)
|
||
try:
|
||
log(f"数据库 {args.db_name} 已连接，准备初始化表结构 ...")
|
||
ensure_schema(conn)
|
||
if args.resume:
|
||
first_code = packages[0].code
|
||
last_code = packages[-1].code
|
||
log(
|
||
f"续跑 layer={layer_name} cell={cell_size_m:g}m；"
|
||
f"保留已入库格子，仅处理 {first_code}..{last_code}"
|
||
)
|
||
else:
|
||
log(f"准备重建 layer={layer_name} cell={cell_size_m:g}m；只清理当前 layer，不影响其它尺寸")
|
||
replace_layer_rows(conn, layer_name)
|
||
|
||
total_lines = 0
|
||
total_points = 0
|
||
package_stats: list[dict] = []
|
||
build_start = dt.datetime.now()
|
||
|
||
for package in packages:
|
||
package_start = dt.datetime.now()
|
||
log(f"[{package.code}] 开始处理 {package.path.name}")
|
||
log(f"[{package.code}] 读取原始曲线并构造几何 ...")
|
||
lines = list(iter_coastline_lines(package.path))
|
||
if not lines:
|
||
log(f"[{package.code}] 没有可用海岸线，跳过")
|
||
continue
|
||
|
||
coastline_geom = MultiLineString([list(line.coords) for line in lines])
|
||
coast_bounds = coastline_geom.bounds
|
||
log(
|
||
f"[{package.code}] 几何已构造 lines={len(lines)} points={sum(len(line.coords) for line in lines)} "
|
||
f"bbox={coast_bounds[0]:.2f},{coast_bounds[1]:.2f},{coast_bounds[2]:.2f},{coast_bounds[3]:.2f} "
|
||
f"buffer={args.buffer_m:g}m"
|
||
)
|
||
|
||
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] = []
|
||
package_seen_cells: set[str] = set()
|
||
scanned_cells = 0
|
||
next_progress_at = args.progress_every_cells
|
||
last_progress_at = dt.datetime.now()
|
||
coarse_factor = args.coarse_factor
|
||
coarse_size = cell_size_m * coarse_factor
|
||
|
||
for line_index, line in enumerate(lines, start=1):
|
||
line_start = dt.datetime.now()
|
||
line_buffer = line.buffer(args.buffer_m, cap_style=2, join_style=2)
|
||
if line_buffer.is_empty:
|
||
log(f"[{package.code}] 线段#{line_index} buffer 为空，跳过")
|
||
continue
|
||
prepared_buffer = prep(line_buffer)
|
||
line_bounds = line_buffer.bounds
|
||
log(
|
||
f"[{package.code}] 线段#{line_index} 扫描开始 "
|
||
f"bbox={line_bounds[0]:.2f},{line_bounds[1]:.2f},{line_bounds[2]:.2f},{line_bounds[3]:.2f}"
|
||
)
|
||
|
||
scan_minx = align_floor(line_bounds[0] - args.margin_cells * cell_size_m, coarse_size)
|
||
scan_miny = align_floor(line_bounds[1] - args.margin_cells * cell_size_m, coarse_size)
|
||
scan_maxx = align_ceil(line_bounds[2] + args.margin_cells * cell_size_m, coarse_size)
|
||
scan_maxy = align_ceil(line_bounds[3] + args.margin_cells * cell_size_m, coarse_size)
|
||
|
||
cix0 = int(round(scan_minx / coarse_size))
|
||
ciy0 = int(round(scan_miny / coarse_size))
|
||
cix1 = int(round(scan_maxx / coarse_size))
|
||
ciy1 = int(round(scan_maxy / coarse_size))
|
||
|
||
for cix in range(cix0, cix1):
|
||
coarse_minx = cix * coarse_size
|
||
coarse_maxx = coarse_minx + coarse_size
|
||
for ciy in range(ciy0, ciy1):
|
||
coarse_miny = ciy * coarse_size
|
||
coarse_maxy = coarse_miny + coarse_size
|
||
coarse_cell = box(coarse_minx, coarse_miny, coarse_maxx, coarse_maxy)
|
||
scanned_cells += 1
|
||
|
||
if not prepared_buffer.intersects(coarse_cell):
|
||
package_candidate += coarse_factor * coarse_factor
|
||
now = dt.datetime.now()
|
||
if scanned_cells >= next_progress_at or (
|
||
args.progress_every_seconds
|
||
and (now - last_progress_at).total_seconds() >= args.progress_every_seconds
|
||
):
|
||
elapsed = (now - package_start).total_seconds()
|
||
rate = scanned_cells / elapsed if elapsed > 0 else 0.0
|
||
log(
|
||
f"[{package.code}] 进度 scanned={scanned_cells} blocked={package_blocked} "
|
||
f"candidate={package_candidate} rows={len(package_rows)} "
|
||
f"rate={rate:.1f}/s elapsed={elapsed:.1f}s"
|
||
)
|
||
last_progress_at = now
|
||
next_progress_at = scanned_cells + args.progress_every_cells
|
||
continue
|
||
|
||
fx0 = cix * coarse_factor
|
||
fy0 = ciy * coarse_factor
|
||
fx1 = fx0 + coarse_factor
|
||
fy1 = fy0 + coarse_factor
|
||
|
||
if prepared_buffer.contains(coarse_cell):
|
||
for ix in range(fx0, fx1):
|
||
cell_minx = ix * cell_size_m
|
||
cell_maxx = cell_minx + cell_size_m
|
||
for iy in range(fy0, fy1):
|
||
cell_miny = iy * cell_size_m
|
||
cell_maxy = cell_miny + cell_size_m
|
||
cell_id = f"{ix}:{iy}"
|
||
if cell_id in package_seen_cells:
|
||
continue
|
||
package_seen_cells.add(cell_id)
|
||
package_rows.append(
|
||
(
|
||
layer_name,
|
||
cell_id,
|
||
iy,
|
||
ix,
|
||
cell_size_m,
|
||
"HARD_BLOCKED",
|
||
"HARD_BLOCKED",
|
||
package.code,
|
||
str(package.path),
|
||
package.code,
|
||
cell_minx,
|
||
cell_miny,
|
||
cell_maxx,
|
||
cell_maxy,
|
||
)
|
||
)
|
||
package_blocked += 1
|
||
scanned_cells += coarse_factor * coarse_factor
|
||
now = dt.datetime.now()
|
||
if scanned_cells >= next_progress_at or (
|
||
args.progress_every_seconds
|
||
and (now - last_progress_at).total_seconds() >= args.progress_every_seconds
|
||
):
|
||
elapsed = (now - package_start).total_seconds()
|
||
rate = scanned_cells / elapsed if elapsed > 0 else 0.0
|
||
log(
|
||
f"[{package.code}] 进度 scanned={scanned_cells} blocked={package_blocked} "
|
||
f"candidate={package_candidate} rows={len(package_rows)} "
|
||
f"rate={rate:.1f}/s elapsed={elapsed:.1f}s"
|
||
)
|
||
last_progress_at = now
|
||
next_progress_at = scanned_cells + args.progress_every_cells
|
||
continue
|
||
|
||
for ix in range(fx0, fx1):
|
||
cell_minx = ix * cell_size_m
|
||
cell_maxx = cell_minx + cell_size_m
|
||
for iy in range(fy0, fy1):
|
||
cell_miny = iy * cell_size_m
|
||
cell_maxy = cell_miny + cell_size_m
|
||
cell_id = f"{ix}:{iy}"
|
||
scanned_cells += 1
|
||
cell = box(cell_minx, cell_miny, cell_maxx, cell_maxy)
|
||
if prepared_buffer.intersects(cell):
|
||
if cell_id not in package_seen_cells:
|
||
package_seen_cells.add(cell_id)
|
||
package_rows.append(
|
||
(
|
||
layer_name,
|
||
cell_id,
|
||
iy,
|
||
ix,
|
||
cell_size_m,
|
||
"HARD_BLOCKED",
|
||
"HARD_BLOCKED",
|
||
package.code,
|
||
str(package.path),
|
||
package.code,
|
||
cell_minx,
|
||
cell_miny,
|
||
cell_maxx,
|
||
cell_maxy,
|
||
)
|
||
)
|
||
package_blocked += 1
|
||
else:
|
||
package_candidate += 1
|
||
|
||
now = dt.datetime.now()
|
||
if scanned_cells >= next_progress_at or (
|
||
args.progress_every_seconds
|
||
and (now - last_progress_at).total_seconds() >= args.progress_every_seconds
|
||
):
|
||
elapsed = (now - package_start).total_seconds()
|
||
rate = scanned_cells / elapsed if elapsed > 0 else 0.0
|
||
log(
|
||
f"[{package.code}] 进度 scanned={scanned_cells} blocked={package_blocked} "
|
||
f"candidate={package_candidate} rows={len(package_rows)} "
|
||
f"rate={rate:.1f}/s elapsed={elapsed:.1f}s"
|
||
)
|
||
last_progress_at = now
|
||
next_progress_at = scanned_cells + args.progress_every_cells
|
||
|
||
if package_rows and len(package_rows) >= 100:
|
||
flush_package_rows(
|
||
conn,
|
||
package_rows,
|
||
package_code=package.code,
|
||
package_blocked=package_blocked,
|
||
package_candidate=package_candidate,
|
||
reason=f"line#{line_index}",
|
||
)
|
||
|
||
if package_rows:
|
||
flush_package_rows(
|
||
conn,
|
||
package_rows,
|
||
package_code=package.code,
|
||
package_blocked=package_blocked,
|
||
package_candidate=package_candidate,
|
||
reason=f"line#{line_index}-end",
|
||
)
|
||
|
||
log(
|
||
f"[{package.code}] 线段#{line_index} 完成 "
|
||
f"blocked={package_blocked} candidate={package_candidate} "
|
||
f"elapsed={(dt.datetime.now() - line_start).total_seconds():.1f}s"
|
||
)
|
||
|
||
if package_rows:
|
||
flush_package_rows(
|
||
conn,
|
||
package_rows,
|
||
package_code=package.code,
|
||
package_blocked=package_blocked,
|
||
package_candidate=package_candidate,
|
||
reason="package-end",
|
||
)
|
||
|
||
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(
|
||
"DELETE FROM navsea_coast_grid_package_stat WHERE layer_name=%s AND source_code=%s",
|
||
(layer_name, package.code),
|
||
)
|
||
cur.execute(
|
||
"""
|
||
INSERT INTO navsea_coast_grid_package_stat(
|
||
layer_name, 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, %s
|
||
)
|
||
""",
|
||
(
|
||
layer_name,
|
||
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()
|
||
log(
|
||
f"[{package.code}] 完成 line={line_count} point={point_count} "
|
||
f"blocked={package_blocked} candidate={package_candidate} "
|
||
f"elapsed={(dt.datetime.now() - package_start).total_seconds():.1f}s"
|
||
)
|
||
|
||
with conn.cursor() as cur:
|
||
cur.execute("ANALYZE TABLE navsea_coast_grid_cell")
|
||
cur.execute("ANALYZE TABLE navsea_coast_grid_package_stat")
|
||
cur.execute("ANALYZE TABLE navsea_coast_grid_layer_meta")
|
||
conn.commit()
|
||
|
||
with conn.cursor() as cur:
|
||
cur.execute(f"SELECT COUNT(*) FROM navsea_coast_grid_cell WHERE layer_name=%s", (layer_name,))
|
||
final_total = int(cur.fetchone()[0] or 0)
|
||
cur.execute(
|
||
f"SELECT COUNT(*) FROM navsea_coast_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_coast_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_coast_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,
|
||
"layer_name": layer_name,
|
||
"cell_size_m": cell_size_m,
|
||
"buffer_m": args.buffer_m,
|
||
"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": f"全国海岸线 {cell_size_m:g}m 硬阻塞栅格，直接写入 MySQL",
|
||
}
|
||
manifest_path = out_dir / f"japan_{layer_name}_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=f"全国海岸线 {cell_size_m:g}m 硬阻塞格",
|
||
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()
|
||
|
||
elapsed_all = (dt.datetime.now() - build_start).total_seconds()
|
||
log("完成")
|
||
log(f" MySQL: {args.db_name}")
|
||
log(f" layer={layer_name}")
|
||
log(f" cell_size_m={cell_size_m:g}")
|
||
log(f" packages={len(packages)}")
|
||
log(f" coastline_lines={total_lines}")
|
||
log(f" grid_cells={final_total}")
|
||
log(f" blocked={blocked_count}")
|
||
log(f" candidate={candidate_count}")
|
||
log(f" elapsed={elapsed_all:.1f}s")
|
||
log(f" manifest: {manifest_path}")
|
||
finally:
|
||
conn.close()
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|