985 lines
35 KiB
Python
985 lines
35 KiB
Python
from __future__ import annotations
|
||
|
||
import argparse
|
||
import datetime as dt
|
||
import math
|
||
import sqlite3
|
||
import subprocess
|
||
from pathlib import Path
|
||
|
||
from shapely.affinity import affine_transform
|
||
from shapely.geometry import box
|
||
from shapely.ops import unary_union
|
||
from shapely.prepared import prep
|
||
from shapely.wkt import loads as load_wkt
|
||
|
||
|
||
DB_NAME = "pbf_analysis"
|
||
OUT_DB = Path("out/navsea_kyushu_navigability_test.sqlite")
|
||
OUT_REPORT = Path("report/navsea_kyushu_navigability_test_2026-04-22.md")
|
||
|
||
MYSQL_BASE = [
|
||
"mysql",
|
||
"-h",
|
||
"localhost",
|
||
"-u",
|
||
"root",
|
||
"-p2chi9ks2",
|
||
"-D",
|
||
DB_NAME,
|
||
"-B",
|
||
"-N",
|
||
"--raw",
|
||
"--default-character-set=utf8mb4",
|
||
]
|
||
|
||
WORLD_SPAN_M = 40075016.68557849
|
||
ORIGIN_SHIFT_M = WORLD_SPAN_M / 2.0
|
||
TARGET_Z = 12
|
||
TILE_EXTENT = 1048576.0
|
||
TILE_SPAN_M = WORLD_SPAN_M / (2**TARGET_Z)
|
||
LOCAL_TO_M = TILE_SPAN_M / TILE_EXTENT
|
||
PORT_AOI_RADIUS_M = 1000.0
|
||
COLLISION_BUFFER_M = 80.0
|
||
GROUNDING_BUFFER_M = 40.0
|
||
BOUNDARY_BUFFER_M = 30.0
|
||
ROUTE_BUFFER_M = 25.0
|
||
GRID_200_M = 200.0
|
||
GRID_50_M = 50.0
|
||
KEY_SHIFT = 20
|
||
KEY_MASK = (1 << KEY_SHIFT) - 1
|
||
PORT_TRIGGER_TYPES = (
|
||
"漁港",
|
||
"一般港湾 (港則法区域、重要港湾)",
|
||
"小港湾 (避難港、地方港湾)",
|
||
)
|
||
TILE_COMMIT_INTERVAL = 100
|
||
GRID_PROFILE_PORT_ONLY = "port-tiles"
|
||
GRID_PROFILE_FULL_DOMAIN = "full-domain"
|
||
DEFAULT_COMMIT_INTERVAL = 20
|
||
|
||
|
||
def meta_upsert(conn: sqlite3.Connection, key: str, value: str) -> None:
|
||
conn.execute("INSERT OR REPLACE INTO meta(key, value) VALUES (?, ?)", (key, value))
|
||
|
||
|
||
def meta_get(conn: sqlite3.Connection, key: str, default: str = "") -> str:
|
||
row = conn.execute("SELECT value FROM meta WHERE key = ?", (key,)).fetchone()
|
||
return str(row[0]) if row else default
|
||
|
||
|
||
def table_exists(conn: sqlite3.Connection, table: str) -> bool:
|
||
row = conn.execute(
|
||
"SELECT 1 FROM sqlite_master WHERE type='table' AND name = ? LIMIT 1",
|
||
(table,),
|
||
).fetchone()
|
||
return row is not None
|
||
|
||
|
||
def table_count(conn: sqlite3.Connection, table: str) -> int:
|
||
if not table_exists(conn, table):
|
||
return 0
|
||
return int(conn.execute(f"SELECT COUNT(*) FROM {table}").fetchone()[0])
|
||
|
||
|
||
def mysql_stream(query: str):
|
||
proc = subprocess.Popen(
|
||
MYSQL_BASE + ["-e", query],
|
||
stdout=subprocess.PIPE,
|
||
stderr=subprocess.PIPE,
|
||
text=True,
|
||
encoding="utf-8",
|
||
errors="replace",
|
||
bufsize=1,
|
||
)
|
||
assert proc.stdout is not None
|
||
assert proc.stderr is not None
|
||
try:
|
||
for raw_line in proc.stdout:
|
||
line = raw_line.rstrip("\n")
|
||
if not line:
|
||
continue
|
||
yield line.split("\t")
|
||
finally:
|
||
rc = proc.wait()
|
||
if rc != 0:
|
||
raise RuntimeError(proc.stderr.read().strip())
|
||
|
||
|
||
def encode_key(col: int, row: int) -> int:
|
||
return (col << KEY_SHIFT) | row
|
||
|
||
|
||
def decode_key(key: int) -> tuple[int, int]:
|
||
return key >> KEY_SHIFT, key & KEY_MASK
|
||
|
||
|
||
def mercator_tile_affine(tile_x: int, tile_y: int) -> list[float]:
|
||
xoff = -ORIGIN_SHIFT_M + tile_x * TILE_SPAN_M
|
||
yoff = ORIGIN_SHIFT_M - tile_y * TILE_SPAN_M
|
||
return [LOCAL_TO_M, 0.0, 0.0, -LOCAL_TO_M, xoff, yoff]
|
||
|
||
|
||
def mercator_point_to_tile(mx: float, my: float, z: int = TARGET_Z) -> tuple[int, int]:
|
||
tiles = 2**z
|
||
tx = int(math.floor(((mx + ORIGIN_SHIFT_M) / WORLD_SPAN_M) * tiles))
|
||
ty = int(math.floor(((ORIGIN_SHIFT_M - my) / WORLD_SPAN_M) * tiles))
|
||
return max(0, min(tiles - 1, tx)), max(0, min(tiles - 1, ty))
|
||
|
||
|
||
def cell_index_range(minx: float, miny: float, maxx: float, maxy: float, grid_size: float) -> tuple[int, int, int, int]:
|
||
col_min = int(math.floor((minx + ORIGIN_SHIFT_M) / grid_size))
|
||
col_max = int(math.floor(((maxx + ORIGIN_SHIFT_M) - 1e-9) / grid_size))
|
||
row_min = int(math.floor((ORIGIN_SHIFT_M - maxy) / grid_size))
|
||
row_max = int(math.floor(((ORIGIN_SHIFT_M - miny) - 1e-9) / grid_size))
|
||
return col_min, col_max, row_min, row_max
|
||
|
||
|
||
def cell_bounds(col: int, row: int, grid_size: float) -> tuple[float, float, float, float]:
|
||
minx = -ORIGIN_SHIFT_M + col * grid_size
|
||
maxx = minx + grid_size
|
||
maxy = ORIGIN_SHIFT_M - row * grid_size
|
||
miny = maxy - grid_size
|
||
return minx, miny, maxx, maxy
|
||
|
||
|
||
def make_cell_row(
|
||
grid_level: int,
|
||
col: int,
|
||
row: int,
|
||
state: str,
|
||
blocked: int,
|
||
reason: str | None,
|
||
risk_score: int,
|
||
confidence: str,
|
||
counts: list[int] | None = None,
|
||
) -> tuple:
|
||
counts = counts or [0, 0, 0, 0, 0, 0]
|
||
minx, miny, maxx, maxy = cell_bounds(col, row, float(grid_level))
|
||
parent_col = col // 4 if grid_level == 50 else None
|
||
parent_row = row // 4 if grid_level == 50 else None
|
||
return (
|
||
grid_level,
|
||
col,
|
||
row,
|
||
parent_col,
|
||
parent_row,
|
||
minx,
|
||
miny,
|
||
maxx,
|
||
maxy,
|
||
state,
|
||
blocked,
|
||
reason,
|
||
risk_score,
|
||
confidence,
|
||
counts[0],
|
||
counts[1],
|
||
counts[2],
|
||
counts[3],
|
||
counts[4],
|
||
counts[5],
|
||
)
|
||
|
||
|
||
def classify_cell(counts: list[int]) -> tuple[str, int, str | None, int]:
|
||
land_hits, hole_hits, collision_hits, grounding_hits, route_hits, boundary_hits = counts
|
||
if land_hits > hole_hits and land_hits > 0:
|
||
return "blocked", 1, "land", 100
|
||
if collision_hits > 0:
|
||
return "blocked", 1, "collision", 90
|
||
if grounding_hits > 0:
|
||
return "caution", 0, "grounding", 60
|
||
if boundary_hits > 0:
|
||
return "caution", 0, "boundary", 40
|
||
if route_hits > 0:
|
||
return "preferred", 0, "route", 20
|
||
return "open_water", 0, None, 0
|
||
|
||
|
||
def create_schema(conn: sqlite3.Connection) -> None:
|
||
if table_exists(conn, "meta"):
|
||
return
|
||
conn.executescript(
|
||
"""
|
||
CREATE TABLE meta (
|
||
key TEXT PRIMARY KEY,
|
||
value TEXT NOT NULL
|
||
);
|
||
|
||
CREATE TABLE tile_default (
|
||
z INTEGER NOT NULL,
|
||
x INTEGER NOT NULL,
|
||
y INTEGER NOT NULL,
|
||
default_state TEXT NOT NULL,
|
||
has_port_refine INTEGER NOT NULL DEFAULT 0,
|
||
PRIMARY KEY (z, x, y)
|
||
);
|
||
|
||
CREATE TABLE port_aoi (
|
||
aoi_id INTEGER PRIMARY KEY AUTOINCREMENT,
|
||
feature_id INTEGER NOT NULL,
|
||
fid TEXT,
|
||
canonical_object_type TEXT NOT NULL,
|
||
center_x REAL NOT NULL,
|
||
center_y REAL NOT NULL,
|
||
radius_m REAL NOT NULL
|
||
);
|
||
|
||
CREATE TABLE grid_cell (
|
||
grid_level INTEGER NOT NULL,
|
||
grid_col INTEGER NOT NULL,
|
||
grid_row INTEGER NOT NULL,
|
||
parent_grid_col INTEGER,
|
||
parent_grid_row INTEGER,
|
||
minx REAL NOT NULL,
|
||
miny REAL NOT NULL,
|
||
maxx REAL NOT NULL,
|
||
maxy REAL NOT NULL,
|
||
state TEXT NOT NULL,
|
||
blocked INTEGER NOT NULL DEFAULT 0,
|
||
state_reason TEXT,
|
||
risk_score INTEGER NOT NULL DEFAULT 0,
|
||
confidence TEXT NOT NULL,
|
||
land_hits INTEGER NOT NULL DEFAULT 0,
|
||
hole_hits INTEGER NOT NULL DEFAULT 0,
|
||
collision_hits INTEGER NOT NULL DEFAULT 0,
|
||
grounding_hits INTEGER NOT NULL DEFAULT 0,
|
||
route_hits INTEGER NOT NULL DEFAULT 0,
|
||
boundary_hits INTEGER NOT NULL DEFAULT 0,
|
||
PRIMARY KEY (grid_level, grid_col, grid_row)
|
||
);
|
||
|
||
CREATE INDEX idx_grid_cell_state ON grid_cell(grid_level, state, blocked);
|
||
CREATE INDEX idx_grid_cell_parent ON grid_cell(parent_grid_col, parent_grid_row);
|
||
"""
|
||
)
|
||
conn.executemany(
|
||
"INSERT INTO meta(key, value) VALUES (?, ?)",
|
||
[
|
||
("source_db", DB_NAME),
|
||
("source_scope", "九州"),
|
||
("build_time", dt.datetime.now().isoformat(timespec="seconds")),
|
||
("target_zoom", str(TARGET_Z)),
|
||
("tile_extent", str(int(TILE_EXTENT))),
|
||
("grid_profile", GRID_PROFILE_FULL_DOMAIN),
|
||
("port_aoi_radius_m", str(int(PORT_AOI_RADIUS_M))),
|
||
("collision_buffer_m", str(int(COLLISION_BUFFER_M))),
|
||
("grounding_buffer_m", str(int(GROUNDING_BUFFER_M))),
|
||
("boundary_buffer_m", str(int(BOUNDARY_BUFFER_M))),
|
||
("route_buffer_m", str(int(ROUTE_BUFFER_M))),
|
||
("notes", "测试版；支持九州全域 200m 主层和港口 50m 细化"),
|
||
],
|
||
)
|
||
conn.commit()
|
||
|
||
|
||
def load_tile_defaults(conn: sqlite3.Connection) -> int:
|
||
existing = table_count(conn, "tile_default")
|
||
if existing > 0:
|
||
return existing
|
||
batch: list[tuple[int, int, int, str, int]] = []
|
||
count = 0
|
||
for cols in mysql_stream(
|
||
f"""
|
||
SELECT {TARGET_Z} AS z, x, y
|
||
FROM feature_geojson_lookup
|
||
WHERE z = {TARGET_Z}
|
||
GROUP BY x, y
|
||
ORDER BY x, y
|
||
"""
|
||
):
|
||
z, x, y = map(int, cols)
|
||
batch.append((z, x, y, "open_water", 0))
|
||
count += 1
|
||
if len(batch) >= 5000:
|
||
conn.executemany(
|
||
"""
|
||
INSERT INTO tile_default(z, x, y, default_state, has_port_refine)
|
||
VALUES (?, ?, ?, ?, ?)
|
||
""",
|
||
batch,
|
||
)
|
||
conn.commit()
|
||
batch.clear()
|
||
if batch:
|
||
conn.executemany(
|
||
"""
|
||
INSERT INTO tile_default(z, x, y, default_state, has_port_refine)
|
||
VALUES (?, ?, ?, ?, ?)
|
||
""",
|
||
batch,
|
||
)
|
||
conn.commit()
|
||
return count
|
||
|
||
|
||
def build_port_aoi(conn: sqlite3.Connection) -> tuple[int, int, set[int], set[tuple[int, int]]]:
|
||
existing_features = table_count(conn, "port_aoi")
|
||
existing_cells = conn.execute("SELECT COUNT(*) FROM grid_cell WHERE grid_level = 50").fetchone()[0] if table_exists(conn, "grid_cell") else 0
|
||
existing_tiles = {
|
||
(int(x), int(y))
|
||
for x, y in conn.execute(
|
||
f"SELECT x, y FROM tile_default WHERE z = {TARGET_Z} AND has_port_refine = 1"
|
||
).fetchall()
|
||
}
|
||
if existing_features > 0 and existing_cells > 0 and existing_tiles:
|
||
port_cell_keys = {
|
||
encode_key(int(col), int(row))
|
||
for col, row in conn.execute(
|
||
"SELECT grid_col, grid_row FROM grid_cell WHERE grid_level = 50"
|
||
).fetchall()
|
||
}
|
||
return existing_features, existing_cells, port_cell_keys, existing_tiles
|
||
|
||
placeholders = ",".join(f"'{value}'" for value in PORT_TRIGGER_TYPES)
|
||
query = f"""
|
||
SELECT
|
||
d.feature_id,
|
||
gl.fid,
|
||
d.canonical_object_type,
|
||
d.x,
|
||
d.y,
|
||
g.geometry_wkt
|
||
FROM navsea_detection_objects d
|
||
JOIN feature_geojson_lookup gl
|
||
ON gl.feature_id = d.feature_id
|
||
JOIN feature_geometry_stage g
|
||
ON g.feature_id = d.feature_id
|
||
WHERE d.z = {TARGET_Z}
|
||
AND d.canonical_object_type IN ({placeholders})
|
||
ORDER BY d.x, d.y, d.feature_id
|
||
"""
|
||
inserted_cells = 0
|
||
port_feature_count = 0
|
||
port_cell_keys: set[int] = set()
|
||
port_tiles: set[tuple[int, int]] = set()
|
||
cell_batch: list[tuple] = []
|
||
aoi_batch: list[tuple] = []
|
||
|
||
for cols in mysql_stream(query):
|
||
if len(cols) != 6:
|
||
raise RuntimeError(f"unexpected port trigger column count: {len(cols)}")
|
||
feature_id, fid, canonical_type, tile_x, tile_y, geometry_wkt = cols
|
||
tile_x_i = int(tile_x)
|
||
tile_y_i = int(tile_y)
|
||
raw_geom = load_wkt(geometry_wkt)
|
||
mercator_geom = affine_transform(raw_geom, mercator_tile_affine(tile_x_i, tile_y_i))
|
||
center = mercator_geom.centroid
|
||
center_x = float(center.x)
|
||
center_y = float(center.y)
|
||
aoi_batch.append((int(feature_id), fid or None, canonical_type, center_x, center_y, PORT_AOI_RADIUS_M))
|
||
port_feature_count += 1
|
||
|
||
minx = center_x - PORT_AOI_RADIUS_M
|
||
miny = center_y - PORT_AOI_RADIUS_M
|
||
maxx = center_x + PORT_AOI_RADIUS_M
|
||
maxy = center_y + PORT_AOI_RADIUS_M
|
||
col_min, col_max, row_min, row_max = cell_index_range(minx, miny, maxx, maxy, GRID_50_M)
|
||
|
||
tile_min = mercator_point_to_tile(minx, maxy)
|
||
tile_max = mercator_point_to_tile(maxx, miny)
|
||
for tx in range(tile_min[0], tile_max[0] + 1):
|
||
for ty in range(tile_min[1], tile_max[1] + 1):
|
||
port_tiles.add((tx, ty))
|
||
|
||
for col in range(col_min, col_max + 1):
|
||
for row in range(row_min, row_max + 1):
|
||
key = encode_key(col, row)
|
||
if key in port_cell_keys:
|
||
continue
|
||
cminx, cminy, cmaxx, cmaxy = cell_bounds(col, row, GRID_50_M)
|
||
center_cx = (cminx + cmaxx) / 2.0
|
||
center_cy = (cminy + cmaxy) / 2.0
|
||
if math.hypot(center_cx - center_x, center_cy - center_y) > PORT_AOI_RADIUS_M:
|
||
continue
|
||
port_cell_keys.add(key)
|
||
inserted_cells += 1
|
||
cell_batch.append(
|
||
make_cell_row(
|
||
50,
|
||
col,
|
||
row,
|
||
"open_water",
|
||
0,
|
||
None,
|
||
0,
|
||
"port_aoi_default",
|
||
)
|
||
)
|
||
if len(cell_batch) >= 10000:
|
||
conn.executemany(
|
||
"""
|
||
INSERT OR IGNORE INTO grid_cell(
|
||
grid_level, grid_col, grid_row, parent_grid_col, parent_grid_row,
|
||
minx, miny, maxx, maxy, state, blocked, state_reason, risk_score, confidence,
|
||
land_hits, hole_hits, collision_hits, grounding_hits, route_hits, boundary_hits
|
||
) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)
|
||
""",
|
||
cell_batch,
|
||
)
|
||
conn.commit()
|
||
cell_batch.clear()
|
||
if len(aoi_batch) >= 2000:
|
||
conn.executemany(
|
||
"""
|
||
INSERT INTO port_aoi(
|
||
feature_id, fid, canonical_object_type, center_x, center_y, radius_m
|
||
) VALUES (?, ?, ?, ?, ?, ?)
|
||
""",
|
||
aoi_batch,
|
||
)
|
||
conn.commit()
|
||
aoi_batch.clear()
|
||
|
||
if cell_batch:
|
||
conn.executemany(
|
||
"""
|
||
INSERT OR IGNORE INTO grid_cell(
|
||
grid_level, grid_col, grid_row, parent_grid_col, parent_grid_row,
|
||
minx, miny, maxx, maxy, state, blocked, state_reason, risk_score, confidence,
|
||
land_hits, hole_hits, collision_hits, grounding_hits, route_hits, boundary_hits
|
||
) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)
|
||
""",
|
||
cell_batch,
|
||
)
|
||
conn.commit()
|
||
if aoi_batch:
|
||
conn.executemany(
|
||
"""
|
||
INSERT INTO port_aoi(
|
||
feature_id, fid, canonical_object_type, center_x, center_y, radius_m
|
||
) VALUES (?, ?, ?, ?, ?, ?)
|
||
""",
|
||
aoi_batch,
|
||
)
|
||
conn.commit()
|
||
|
||
if port_tiles:
|
||
conn.executemany(
|
||
f"""
|
||
UPDATE tile_default
|
||
SET has_port_refine = 1
|
||
WHERE z = {TARGET_Z} AND x = ? AND y = ?
|
||
""",
|
||
sorted(port_tiles),
|
||
)
|
||
conn.commit()
|
||
|
||
return port_feature_count, inserted_cells, port_cell_keys, port_tiles
|
||
|
||
|
||
def relevant_feature_query(last_tile: tuple[int, int] | None = None) -> str:
|
||
resume_filter = ""
|
||
if last_tile is not None:
|
||
last_x, last_y = last_tile
|
||
resume_filter = f"\n AND (d.x > {last_x} OR (d.x = {last_x} AND d.y > {last_y}))"
|
||
return f"""
|
||
SELECT
|
||
d.feature_id,
|
||
d.canonical_object_type,
|
||
d.capability,
|
||
d.x,
|
||
d.y,
|
||
d.geom_type,
|
||
g.geometry_wkt
|
||
FROM navsea_detection_objects d
|
||
JOIN feature_geometry_stage g
|
||
ON g.feature_id = d.feature_id
|
||
WHERE d.z = {TARGET_Z}
|
||
AND (
|
||
d.canonical_object_type IN ('P陸域', 'P穴')
|
||
OR d.capability IN ('collision', 'grounding', 'route_reference', 'boundary_reference')
|
||
)
|
||
{resume_filter}
|
||
ORDER BY d.x, d.y, d.feature_id
|
||
"""
|
||
|
||
|
||
def touch_bucket(bucket: dict[int, list[int]], key: int, index: int) -> None:
|
||
counts = bucket.get(key)
|
||
if counts is None:
|
||
counts = [0, 0, 0, 0, 0, 0]
|
||
bucket[key] = counts
|
||
counts[index] += 1
|
||
|
||
|
||
def effective_geometry(raw_geom, canonical_type: str, capability: str | None):
|
||
capability = capability or ""
|
||
if canonical_type == "P陸域":
|
||
return raw_geom, 0
|
||
if canonical_type == "P穴":
|
||
return raw_geom, 1
|
||
if capability == "collision":
|
||
if raw_geom.geom_type in {"Point", "MultiPoint", "LineString", "MultiLineString"}:
|
||
return raw_geom.buffer(COLLISION_BUFFER_M), 2
|
||
return raw_geom.buffer(COLLISION_BUFFER_M / 2.0), 2
|
||
if capability == "grounding":
|
||
if raw_geom.geom_type in {"Point", "MultiPoint", "LineString", "MultiLineString"}:
|
||
return raw_geom.buffer(GROUNDING_BUFFER_M), 3
|
||
return raw_geom, 3
|
||
if capability == "route_reference":
|
||
if raw_geom.geom_type in {"LineString", "MultiLineString"}:
|
||
return raw_geom.buffer(ROUTE_BUFFER_M), 4
|
||
return raw_geom, 4
|
||
if capability == "boundary_reference":
|
||
if raw_geom.geom_type in {"LineString", "MultiLineString"}:
|
||
return raw_geom.buffer(BOUNDARY_BUFFER_M), 5
|
||
return raw_geom, 5
|
||
return None, -1
|
||
|
||
|
||
def apply_feature_to_grid(
|
||
feature_geom,
|
||
grid_size: float,
|
||
target_bucket: dict[int, list[int]],
|
||
bucket_index: int,
|
||
valid_keys: set[int] | None = None,
|
||
) -> None:
|
||
minx, miny, maxx, maxy = feature_geom.bounds
|
||
col_min, col_max, row_min, row_max = cell_index_range(minx, miny, maxx, maxy, grid_size)
|
||
prepared = prep(feature_geom)
|
||
for col in range(col_min, col_max + 1):
|
||
for row in range(row_min, row_max + 1):
|
||
key = encode_key(col, row)
|
||
if valid_keys is not None and key not in valid_keys:
|
||
continue
|
||
cminx, cminy, cmaxx, cmaxy = cell_bounds(col, row, grid_size)
|
||
if prepared.intersects(box(cminx, cminy, cmaxx, cmaxy)):
|
||
touch_bucket(target_bucket, key, bucket_index)
|
||
|
||
|
||
def flush_grid_counts(conn: sqlite3.Connection, grid_level: int, bucket: dict[int, list[int]], confidence: str) -> int:
|
||
if not bucket:
|
||
return 0
|
||
rows: list[tuple] = []
|
||
for key, counts in bucket.items():
|
||
col, row = decode_key(key)
|
||
rows.append(make_cell_row(grid_level, col, row, "open_water", 0, None, 0, confidence, counts))
|
||
conn.executemany(
|
||
"""
|
||
INSERT INTO grid_cell(
|
||
grid_level, grid_col, grid_row, parent_grid_col, parent_grid_row,
|
||
minx, miny, maxx, maxy, state, blocked, state_reason, risk_score, confidence,
|
||
land_hits, hole_hits, collision_hits, grounding_hits, route_hits, boundary_hits
|
||
) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)
|
||
ON CONFLICT(grid_level, grid_col, grid_row) DO UPDATE SET
|
||
confidence=excluded.confidence,
|
||
land_hits=grid_cell.land_hits + excluded.land_hits,
|
||
hole_hits=grid_cell.hole_hits + excluded.hole_hits,
|
||
collision_hits=grid_cell.collision_hits + excluded.collision_hits,
|
||
grounding_hits=grid_cell.grounding_hits + excluded.grounding_hits,
|
||
route_hits=grid_cell.route_hits + excluded.route_hits,
|
||
boundary_hits=grid_cell.boundary_hits + excluded.boundary_hits
|
||
""",
|
||
rows,
|
||
)
|
||
written = len(rows)
|
||
bucket.clear()
|
||
return written
|
||
|
||
|
||
def merge_geometry(geoms: list):
|
||
if not geoms:
|
||
return None
|
||
if len(geoms) == 1:
|
||
return geoms[0]
|
||
return unary_union(geoms)
|
||
|
||
|
||
def flush_tile_features(
|
||
conn: sqlite3.Connection,
|
||
tile_x: int,
|
||
tile_y: int,
|
||
features: list[tuple[str, str | None, object]],
|
||
port_cell_keys: set[int],
|
||
port_tiles: set[tuple[int, int]],
|
||
grid_profile: str,
|
||
) -> tuple[int, int]:
|
||
grouped: list[list] = [[] for _ in range(6)]
|
||
for canonical_type, capability, mercator_geom in features:
|
||
effect_geom, bucket_index = effective_geometry(mercator_geom, canonical_type, capability)
|
||
if effect_geom is None or bucket_index < 0 or effect_geom.is_empty:
|
||
continue
|
||
grouped[bucket_index].append(effect_geom)
|
||
|
||
tile_200: dict[int, list[int]] = {}
|
||
tile_50: dict[int, list[int]] = {}
|
||
in_port_tile = (tile_x, tile_y) in port_tiles
|
||
|
||
for bucket_index, geoms in enumerate(grouped):
|
||
merged = merge_geometry(geoms)
|
||
if merged is None or merged.is_empty:
|
||
continue
|
||
apply_feature_to_grid(merged, GRID_200_M, tile_200, bucket_index)
|
||
if grid_profile in {GRID_PROFILE_PORT_ONLY, GRID_PROFILE_FULL_DOMAIN} and in_port_tile:
|
||
apply_feature_to_grid(merged, GRID_50_M, tile_50, bucket_index, port_cell_keys)
|
||
|
||
confidence_200 = "z12_full_domain_200m" if grid_profile == GRID_PROFILE_FULL_DOMAIN else "z12_sparse_override"
|
||
written_200 = flush_grid_counts(conn, 200, tile_200, confidence_200)
|
||
written_50 = flush_grid_counts(conn, 50, tile_50, "z12_port_refine")
|
||
return written_200, written_50
|
||
|
||
|
||
def write_checkpoint(
|
||
conn: sqlite3.Connection,
|
||
*,
|
||
processed_tiles: int,
|
||
processed_features: int,
|
||
written_200: int,
|
||
written_50: int,
|
||
last_tile: tuple[int, int] | None,
|
||
phase: str,
|
||
) -> None:
|
||
meta_upsert(conn, "run_phase", phase)
|
||
meta_upsert(conn, "processed_tiles", str(processed_tiles))
|
||
meta_upsert(conn, "processed_features", str(processed_features))
|
||
meta_upsert(conn, "written_200", str(written_200))
|
||
meta_upsert(conn, "written_50", str(written_50))
|
||
if last_tile is None:
|
||
meta_upsert(conn, "last_tile_x", "")
|
||
meta_upsert(conn, "last_tile_y", "")
|
||
else:
|
||
meta_upsert(conn, "last_tile_x", str(last_tile[0]))
|
||
meta_upsert(conn, "last_tile_y", str(last_tile[1]))
|
||
meta_upsert(conn, "last_checkpoint_at", dt.datetime.now().isoformat(timespec="seconds"))
|
||
conn.commit()
|
||
|
||
|
||
def read_checkpoint(conn: sqlite3.Connection) -> dict[str, object]:
|
||
last_x = meta_get(conn, "last_tile_x", "")
|
||
last_y = meta_get(conn, "last_tile_y", "")
|
||
last_tile = None
|
||
if last_x != "" and last_y != "":
|
||
last_tile = (int(last_x), int(last_y))
|
||
return {
|
||
"phase": meta_get(conn, "run_phase", "init"),
|
||
"processed_tiles": int(meta_get(conn, "processed_tiles", "0")),
|
||
"processed_features": int(meta_get(conn, "processed_features", "0")),
|
||
"written_200": int(meta_get(conn, "written_200", "0")),
|
||
"written_50": int(meta_get(conn, "written_50", "0")),
|
||
"last_tile": last_tile,
|
||
}
|
||
|
||
|
||
def accumulate_grid_counts(
|
||
conn: sqlite3.Connection,
|
||
port_cell_keys: set[int],
|
||
port_tiles: set[tuple[int, int]],
|
||
*,
|
||
grid_profile: str,
|
||
tile_limit: int | None,
|
||
commit_interval: int,
|
||
checkpoint: dict[str, object] | None = None,
|
||
) -> tuple[int, int, int, int]:
|
||
current_tile: tuple[int, int] | None = None
|
||
current_features: list[tuple[str, str | None, object]] = []
|
||
checkpoint = checkpoint or {}
|
||
processed_features = int(checkpoint.get("processed_features", 0))
|
||
processed_tiles = int(checkpoint.get("processed_tiles", 0))
|
||
written_200 = int(checkpoint.get("written_200", 0))
|
||
written_50 = int(checkpoint.get("written_50", 0))
|
||
last_tile = checkpoint.get("last_tile")
|
||
|
||
for cols in mysql_stream(relevant_feature_query(last_tile)):
|
||
if len(cols) != 7:
|
||
raise RuntimeError(f"unexpected feature column count: {len(cols)}")
|
||
_, canonical_type, capability, tile_x, tile_y, _, geometry_wkt = cols
|
||
tile_x_i = int(tile_x)
|
||
tile_y_i = int(tile_y)
|
||
tile_key = (tile_x_i, tile_y_i)
|
||
if grid_profile == GRID_PROFILE_PORT_ONLY and tile_key not in port_tiles:
|
||
continue
|
||
|
||
if current_tile is not None and tile_key != current_tile:
|
||
w200, w50 = flush_tile_features(
|
||
conn,
|
||
current_tile[0],
|
||
current_tile[1],
|
||
current_features,
|
||
port_cell_keys,
|
||
port_tiles,
|
||
grid_profile,
|
||
)
|
||
written_200 += w200
|
||
written_50 += w50
|
||
processed_tiles += 1
|
||
if processed_tiles % commit_interval == 0:
|
||
write_checkpoint(
|
||
conn,
|
||
processed_tiles=processed_tiles,
|
||
processed_features=processed_features,
|
||
written_200=written_200,
|
||
written_50=written_50,
|
||
last_tile=current_tile,
|
||
phase="running",
|
||
)
|
||
print(
|
||
f"已处理 tile={processed_tiles} / features={processed_features}，"
|
||
f"累计写回 200m={written_200} / 50m={written_50}",
|
||
flush=True,
|
||
)
|
||
current_features = []
|
||
if tile_limit is not None and processed_tiles >= tile_limit:
|
||
current_tile = None
|
||
break
|
||
|
||
raw_geom = load_wkt(geometry_wkt)
|
||
mercator_geom = affine_transform(raw_geom, mercator_tile_affine(tile_x_i, tile_y_i))
|
||
current_features.append((canonical_type, capability, mercator_geom))
|
||
current_tile = tile_key
|
||
processed_features += 1
|
||
|
||
if current_tile is not None:
|
||
w200, w50 = flush_tile_features(
|
||
conn,
|
||
current_tile[0],
|
||
current_tile[1],
|
||
current_features,
|
||
port_cell_keys,
|
||
port_tiles,
|
||
grid_profile,
|
||
)
|
||
written_200 += w200
|
||
written_50 += w50
|
||
processed_tiles += 1
|
||
|
||
write_checkpoint(
|
||
conn,
|
||
processed_tiles=processed_tiles,
|
||
processed_features=processed_features,
|
||
written_200=written_200,
|
||
written_50=written_50,
|
||
last_tile=current_tile,
|
||
phase="grid_complete",
|
||
)
|
||
return processed_features, processed_tiles, written_200, written_50
|
||
|
||
|
||
def reclassify_grid(conn: sqlite3.Connection) -> None:
|
||
conn.execute(
|
||
"""
|
||
UPDATE grid_cell
|
||
SET
|
||
state = CASE
|
||
WHEN land_hits > hole_hits AND land_hits > 0 THEN 'blocked'
|
||
WHEN collision_hits > 0 THEN 'blocked'
|
||
WHEN grounding_hits > 0 THEN 'caution'
|
||
WHEN boundary_hits > 0 THEN 'caution'
|
||
WHEN route_hits > 0 THEN 'preferred'
|
||
ELSE 'open_water'
|
||
END,
|
||
blocked = CASE
|
||
WHEN land_hits > hole_hits AND land_hits > 0 THEN 1
|
||
WHEN collision_hits > 0 THEN 1
|
||
ELSE 0
|
||
END,
|
||
state_reason = CASE
|
||
WHEN land_hits > hole_hits AND land_hits > 0 THEN 'land'
|
||
WHEN collision_hits > 0 THEN 'collision'
|
||
WHEN grounding_hits > 0 THEN 'grounding'
|
||
WHEN boundary_hits > 0 THEN 'boundary'
|
||
WHEN route_hits > 0 THEN 'route'
|
||
ELSE NULL
|
||
END,
|
||
risk_score = CASE
|
||
WHEN land_hits > hole_hits AND land_hits > 0 THEN 100
|
||
WHEN collision_hits > 0 THEN 90
|
||
WHEN grounding_hits > 0 THEN 60
|
||
WHEN boundary_hits > 0 THEN 40
|
||
WHEN route_hits > 0 THEN 20
|
||
ELSE 0
|
||
END
|
||
"""
|
||
)
|
||
conn.commit()
|
||
|
||
|
||
def read_meta(conn: sqlite3.Connection, key: str, default: str = "") -> str:
|
||
row = conn.execute("SELECT value FROM meta WHERE key = ?", (key,)).fetchone()
|
||
return str(row[0]) if row else default
|
||
|
||
|
||
def build_report(report_path: Path, db_path: Path) -> None:
|
||
conn = sqlite3.connect(db_path)
|
||
grid_profile = read_meta(conn, "grid_profile", GRID_PROFILE_PORT_ONLY)
|
||
tile_limit = read_meta(conn, "tile_limit", "")
|
||
stats = {
|
||
"tile_total": conn.execute("SELECT COUNT(*) FROM tile_default").fetchone()[0],
|
||
"tile_port_refine": conn.execute("SELECT COUNT(*) FROM tile_default WHERE has_port_refine = 1").fetchone()[0],
|
||
"port_aoi_total": conn.execute("SELECT COUNT(*) FROM port_aoi").fetchone()[0],
|
||
"grid_200_total": conn.execute("SELECT COUNT(*) FROM grid_cell WHERE grid_level = 200").fetchone()[0],
|
||
"grid_50_total": conn.execute("SELECT COUNT(*) FROM grid_cell WHERE grid_level = 50").fetchone()[0],
|
||
}
|
||
state_rows = conn.execute(
|
||
"""
|
||
SELECT grid_level, state, COUNT(*)
|
||
FROM grid_cell
|
||
GROUP BY grid_level, state
|
||
ORDER BY grid_level, state
|
||
"""
|
||
).fetchall()
|
||
conn.close()
|
||
|
||
lines = [
|
||
"# 九州可航行性栅格测试版",
|
||
"",
|
||
f"- 生成时间：`{dt.datetime.now().isoformat(timespec='seconds')}`",
|
||
f"- 输出 SQLite：`{db_path}`",
|
||
f"- 栅格模式：`{grid_profile}`",
|
||
f"- tile 限制：`{tile_limit or 'none'}`",
|
||
"",
|
||
"## 1. 当前产物结构",
|
||
"",
|
||
"- `tile_default`：所有 `z12` 瓦片的默认状态，当前默认 `open_water`",
|
||
"- `port_aoi`：港口触发点与 50m 细化半径",
|
||
"- `grid_cell`：保存 `200m` 主层和港口 `50m` 细化结果",
|
||
"",
|
||
"## 2. 统计",
|
||
"",
|
||
f"- `z12` tile 总数：`{stats['tile_total']}`",
|
||
f"- 启用港口 50m 细化的 tile 数：`{stats['tile_port_refine']}`",
|
||
f"- 港口 AOI 触发点数：`{stats['port_aoi_total']}`",
|
||
f"- `200m` 覆盖格子数：`{stats['grid_200_total']}`",
|
||
f"- `50m` 细化格子数：`{stats['grid_50_total']}`",
|
||
"",
|
||
"### 2.1 格子状态分布",
|
||
"",
|
||
"| grid_level | state | count |",
|
||
"| --- | --- | ---: |",
|
||
]
|
||
for grid_level, state, count in state_rows:
|
||
lines.append(f"| `{grid_level}` | `{state}` | `{count}` |")
|
||
lines.extend(
|
||
[
|
||
"",
|
||
"## 3. 当前测试口径",
|
||
"",
|
||
"- `P陸域`：直接作为 `blocked=1`",
|
||
"- `P穴`：用于抵消陆地区域覆盖",
|
||
"- `collision`：按 `80m` 缓冲后作为 `blocked=1`",
|
||
"- `grounding`：按 `40m` 缓冲后作为 `caution`",
|
||
"- `boundary_reference`：按 `30m` 缓冲后作为 `caution`",
|
||
"- `route_reference`：按 `25m` 缓冲后作为 `preferred`",
|
||
"",
|
||
"## 4. 当前限制",
|
||
"",
|
||
"- 这是测试版，不是最终船型/吃水感知版",
|
||
"- 当前 `200m` 层以有语义命中的格子为主，不把纯开放水域全量预烘焙进 `grid_cell`",
|
||
"- 港口 `50m` 层是完整 AOI 细化；AOI 外不存完整 `50m` 全覆盖",
|
||
"- 当前仍是 Web Mercator 米制近似，不是分港区本地投影高精版",
|
||
]
|
||
)
|
||
report_path.parent.mkdir(parents=True, exist_ok=True)
|
||
report_path.write_text("\n".join(lines) + "\n", encoding="utf-8")
|
||
|
||
|
||
def open_db(path: Path, *, resume: bool) -> sqlite3.Connection:
|
||
if path.exists() and not resume:
|
||
path.unlink()
|
||
path.parent.mkdir(parents=True, exist_ok=True)
|
||
conn = sqlite3.connect(path)
|
||
conn.execute("PRAGMA journal_mode=WAL")
|
||
conn.execute("PRAGMA synchronous=OFF")
|
||
conn.execute("PRAGMA temp_store=MEMORY")
|
||
conn.execute("PRAGMA cache_size=200000")
|
||
conn.execute("PRAGMA busy_timeout=5000")
|
||
return conn
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser(description="Build Kyushu navigability test grid from z12 NavSea geometries.")
|
||
parser.add_argument("--out-db", type=Path, default=OUT_DB)
|
||
parser.add_argument("--report", type=Path, default=OUT_REPORT)
|
||
parser.add_argument(
|
||
"--grid-profile",
|
||
choices=[GRID_PROFILE_PORT_ONLY, GRID_PROFILE_FULL_DOMAIN],
|
||
default=GRID_PROFILE_FULL_DOMAIN,
|
||
help="200m 主层覆盖范围；full-domain 为九州全域，port-tiles 仅港口触发 tile。",
|
||
)
|
||
parser.add_argument("--tile-limit", type=int, default=None, help="仅处理前 N 个 tile，用于烟测。")
|
||
parser.add_argument("--resume", action="store_true", help="从已有 SQLite checkpoint 继续执行。")
|
||
parser.add_argument(
|
||
"--commit-interval",
|
||
type=int,
|
||
default=DEFAULT_COMMIT_INTERVAL,
|
||
help="每处理多少个 tile 写一次 checkpoint。",
|
||
)
|
||
args = parser.parse_args()
|
||
|
||
conn = open_db(args.out_db, resume=args.resume)
|
||
create_schema(conn)
|
||
if not args.resume:
|
||
conn.executemany(
|
||
"INSERT OR REPLACE INTO meta(key, value) VALUES (?, ?)",
|
||
[
|
||
("grid_profile", args.grid_profile),
|
||
("tile_limit", "" if args.tile_limit is None else str(args.tile_limit)),
|
||
("commit_interval", str(args.commit_interval)),
|
||
("run_phase", "init"),
|
||
("processed_tiles", "0"),
|
||
("processed_features", "0"),
|
||
("written_200", "0"),
|
||
("written_50", "0"),
|
||
("last_tile_x", ""),
|
||
("last_tile_y", ""),
|
||
],
|
||
)
|
||
conn.commit()
|
||
else:
|
||
args.grid_profile = meta_get(conn, "grid_profile", args.grid_profile)
|
||
tile_limit_meta = meta_get(conn, "tile_limit", "")
|
||
args.tile_limit = int(tile_limit_meta) if tile_limit_meta else args.tile_limit
|
||
args.commit_interval = int(meta_get(conn, "commit_interval", str(args.commit_interval)))
|
||
checkpoint = read_checkpoint(conn)
|
||
print(
|
||
f"恢复执行: phase={checkpoint['phase']} "
|
||
f"processed_tiles={checkpoint['processed_tiles']} "
|
||
f"last_tile={checkpoint['last_tile']}",
|
||
flush=True,
|
||
)
|
||
|
||
tile_count = load_tile_defaults(conn)
|
||
meta_upsert(conn, "run_phase", "tile_default_ready")
|
||
conn.commit()
|
||
print(f"已写入 {tile_count} 个 z12 tile 默认状态", flush=True)
|
||
|
||
port_feature_count, port_cell_count, port_cell_keys, port_tiles = build_port_aoi(conn)
|
||
meta_upsert(conn, "run_phase", "port_aoi_ready")
|
||
conn.commit()
|
||
print(
|
||
f"已建立港口 50m AOI: port_features={port_feature_count}, "
|
||
f"port_cells={port_cell_count}, port_tiles={len(port_tiles)}",
|
||
flush=True,
|
||
)
|
||
|
||
processed, processed_tiles, written_200, written_50 = accumulate_grid_counts(
|
||
conn,
|
||
port_cell_keys,
|
||
port_tiles,
|
||
grid_profile=args.grid_profile,
|
||
tile_limit=args.tile_limit,
|
||
commit_interval=args.commit_interval,
|
||
checkpoint=read_checkpoint(conn) if args.resume else None,
|
||
)
|
||
print(
|
||
f"已完成对象落格累计写回: processed_tiles={processed_tiles}, processed_features={processed}, "
|
||
f"grid200_writes={written_200}, grid50_writes={written_50}",
|
||
flush=True,
|
||
)
|
||
reclassify_grid(conn)
|
||
meta_upsert(conn, "run_phase", "done")
|
||
meta_upsert(conn, "finished_at", dt.datetime.now().isoformat(timespec="seconds"))
|
||
conn.execute("ANALYZE")
|
||
conn.commit()
|
||
conn.close()
|
||
|
||
build_report(args.report, args.out_db)
|
||
print(f"完成: 200m rows={written_200}, 50m upserts={written_50}", flush=True)
|
||
print(f"输出 SQLite: {args.out_db.resolve()}", flush=True)
|
||
print(f"输出报告: {args.report.resolve()}", flush=True)
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|