519 lines
17 KiB
Python
519 lines
17 KiB
Python
#!/usr/bin/env python3
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import datetime as dt
|
|
import json
|
|
import math
|
|
import sys
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
import pymysql
|
|
from shapely import get_parts, intersects_xy
|
|
from shapely.geometry import LineString, box
|
|
from shapely.ops import polygonize, unary_union
|
|
from shapely.prepared import prep
|
|
|
|
PROJECT_ROOT = Path(__file__).resolve().parent.parent
|
|
if str(PROJECT_ROOT) not in sys.path:
|
|
sys.path.insert(0, str(PROJECT_ROOT))
|
|
|
|
from coastline.build_fish_port_20m_full_mysql_resume import (
|
|
align_ceil,
|
|
align_floor,
|
|
build_point_lookup,
|
|
cluster_line_indices,
|
|
collect_fish_lines_for_prc,
|
|
discover_target_prcs,
|
|
load_coastlines_for_prc,
|
|
load_coarse_windows_for_bbox,
|
|
load_xml,
|
|
lon_from_mercator,
|
|
lat_from_mercator,
|
|
mercator_x,
|
|
mercator_y,
|
|
rect_lonlat,
|
|
)
|
|
|
|
|
|
DB_NAME = "navsea_japan_coast_grid"
|
|
DB_USER = "root"
|
|
DB_PASSWORD = "2chi9ks2"
|
|
DB_HOST = "localhost"
|
|
DB_SOCKET = "/tmp/mysql.sock"
|
|
|
|
FISH_SRC = "coastline/C09-06.zip"
|
|
DEFAULT_OUT_DIR = "src/pbf/coastline-mysql/fish_port_debug"
|
|
DEFAULT_BUFFER_M = 20.0
|
|
DEFAULT_COARSE_MARGIN_M = 0.0
|
|
COAST_COARSE_LAYER = "coast_200m"
|
|
CELL_SIZE_M = 20.0
|
|
RADIUS = 6378137.0
|
|
|
|
|
|
def mysql_connect(database: str, *, cursorclass=pymysql.cursors.Cursor):
|
|
kwargs = {
|
|
"host": DB_HOST,
|
|
"user": DB_USER,
|
|
"password": DB_PASSWORD,
|
|
"database": database,
|
|
"charset": "utf8mb4",
|
|
"cursorclass": cursorclass,
|
|
}
|
|
if DB_SOCKET:
|
|
kwargs["unix_socket"] = DB_SOCKET
|
|
return pymysql.connect(**kwargs)
|
|
|
|
|
|
def mercator_to_lonlat_bbox(bbox: tuple[float, float, float, float]) -> tuple[float, float, float, float]:
|
|
minx, miny, maxx, maxy = bbox
|
|
return lon_from_mercator(minx), lat_from_mercator(miny), lon_from_mercator(maxx), lat_from_mercator(maxy)
|
|
|
|
|
|
def mercator_geom_to_lonlat_geojson(geom) -> dict:
|
|
geo = geom.__geo_interface__
|
|
|
|
def convert_coords(coords):
|
|
if not isinstance(coords, (list, tuple)):
|
|
return coords
|
|
if coords and isinstance(coords[0], (int, float)):
|
|
x = float(coords[0])
|
|
y = float(coords[1])
|
|
return [float(lon_from_mercator(x)), float(lat_from_mercator(y))]
|
|
return [convert_coords(item) for item in coords]
|
|
|
|
if geo["type"] == "GeometryCollection":
|
|
return {
|
|
"type": "GeometryCollection",
|
|
"geometries": [mercator_geom_to_lonlat_geojson(part) for part in geom.geoms],
|
|
}
|
|
return {
|
|
"type": geo["type"],
|
|
"coordinates": convert_coords(geo["coordinates"]),
|
|
}
|
|
|
|
|
|
def lonlat_line_feature(
|
|
coords: list[tuple[float, float]],
|
|
*,
|
|
kind: str,
|
|
fpc: str,
|
|
prc: str,
|
|
label: str | None = None,
|
|
) -> dict:
|
|
properties = {
|
|
"kind": kind,
|
|
"fpc": fpc,
|
|
"prc": prc,
|
|
}
|
|
if label is not None:
|
|
properties["label"] = label
|
|
return {
|
|
"type": "Feature",
|
|
"properties": properties,
|
|
"geometry": {
|
|
"type": "LineString",
|
|
"coordinates": [[float(lon), float(lat)] for lon, lat in coords],
|
|
},
|
|
}
|
|
|
|
|
|
def mercator_line_feature(
|
|
line: LineString,
|
|
*,
|
|
kind: str,
|
|
fpc: str,
|
|
prc: str,
|
|
label: str | None = None,
|
|
) -> dict:
|
|
coords = [[float(lon_from_mercator(x)), float(lat_from_mercator(y))] for x, y in line.coords]
|
|
return {
|
|
"type": "Feature",
|
|
"properties": {
|
|
"kind": kind,
|
|
"fpc": fpc,
|
|
"prc": prc,
|
|
**({"label": label} if label is not None else {}),
|
|
},
|
|
"geometry": {
|
|
"type": "LineString",
|
|
"coordinates": coords,
|
|
},
|
|
}
|
|
|
|
|
|
def bbox_feature(
|
|
bbox: tuple[float, float, float, float],
|
|
*,
|
|
kind: str,
|
|
fpc: str,
|
|
prc: str,
|
|
label: str,
|
|
) -> dict:
|
|
min_lon, min_lat, max_lon, max_lat = bbox
|
|
return {
|
|
"type": "Feature",
|
|
"properties": {
|
|
"kind": kind,
|
|
"fpc": fpc,
|
|
"prc": prc,
|
|
"label": label,
|
|
},
|
|
"geometry": rect_lonlat(min_lon, min_lat, max_lon, max_lat),
|
|
}
|
|
|
|
|
|
def mercator_bbox_feature(
|
|
bbox: tuple[float, float, float, float],
|
|
*,
|
|
kind: str,
|
|
fpc: str,
|
|
prc: str,
|
|
label: str,
|
|
) -> dict:
|
|
return bbox_feature(
|
|
mercator_to_lonlat_bbox(bbox),
|
|
kind=kind,
|
|
fpc=fpc,
|
|
prc=prc,
|
|
label=label,
|
|
)
|
|
|
|
|
|
def buffer_tag(value: float) -> str:
|
|
if abs(value - round(value)) < 1e-9:
|
|
return str(int(round(value)))
|
|
return str(value).replace(".", "p")
|
|
|
|
|
|
def scan_debug_cells(
|
|
scan_geom,
|
|
coast_prepared,
|
|
*,
|
|
prc: str,
|
|
fpc: str,
|
|
cluster_idx: int,
|
|
part_idx: int,
|
|
window_idx: int,
|
|
max_scan_cells: int | None,
|
|
) -> tuple[list[dict], int, int]:
|
|
fish_prepared = prep(scan_geom)
|
|
minx, miny, maxx, maxy = scan_geom.bounds
|
|
start_x = align_floor(minx, CELL_SIZE_M)
|
|
start_y = align_floor(miny, CELL_SIZE_M)
|
|
end_x = align_ceil(maxx, CELL_SIZE_M)
|
|
end_y = align_ceil(maxy, CELL_SIZE_M)
|
|
estimated_cols = max(0, int(math.ceil((end_x - start_x) / CELL_SIZE_M)))
|
|
estimated_rows = max(0, int(math.ceil((end_y - start_y) / CELL_SIZE_M)))
|
|
estimated_cells = estimated_cols * estimated_rows
|
|
if max_scan_cells is not None and estimated_cells > max_scan_cells:
|
|
raise RuntimeError(
|
|
"scan guard triggered: "
|
|
f"estimated_cells={estimated_cells} max_scan_cells={max_scan_cells} "
|
|
f"grid={estimated_cols}x{estimated_rows} bounds={tuple(round(v, 3) for v in scan_geom.bounds)} "
|
|
f"fpc={fpc} prc={prc} cluster={cluster_idx} part={part_idx} window={window_idx}"
|
|
)
|
|
|
|
half = CELL_SIZE_M / 2.0
|
|
ix_start = int(round(start_x / CELL_SIZE_M))
|
|
ix_end = int(round(end_x / CELL_SIZE_M))
|
|
ix_values = np.arange(ix_start, ix_end, dtype=np.int64)
|
|
x_values = ix_values.astype(float) * CELL_SIZE_M
|
|
x_center_values = x_values + half
|
|
x_right_values = x_values + CELL_SIZE_M
|
|
|
|
features: list[dict] = []
|
|
coast_count = 0
|
|
sea_count = 0
|
|
|
|
y = start_y
|
|
while y < end_y:
|
|
row_strip = box(start_x, y, end_x, y + CELL_SIZE_M)
|
|
if not fish_prepared.intersects(row_strip):
|
|
y += CELL_SIZE_M
|
|
continue
|
|
|
|
mask = (
|
|
intersects_xy(scan_geom, x_center_values, y + half)
|
|
| intersects_xy(scan_geom, x_values, y)
|
|
| intersects_xy(scan_geom, x_right_values, y)
|
|
| intersects_xy(scan_geom, x_values, y + CELL_SIZE_M)
|
|
| intersects_xy(scan_geom, x_right_values, y + CELL_SIZE_M)
|
|
)
|
|
if not np.any(mask):
|
|
y += CELL_SIZE_M
|
|
continue
|
|
|
|
mask = mask | np.r_[False, mask[:-1]] | np.r_[mask[1:], False]
|
|
for idx in np.flatnonzero(mask):
|
|
x = float(x_values[idx])
|
|
cell_geom = box(x, y, x + CELL_SIZE_M, y + CELL_SIZE_M)
|
|
if not fish_prepared.intersects(cell_geom):
|
|
continue
|
|
coast_hit = coast_prepared.intersects(cell_geom)
|
|
min_lon = lon_from_mercator(x)
|
|
min_lat = lat_from_mercator(y)
|
|
max_lon = lon_from_mercator(x + CELL_SIZE_M)
|
|
max_lat = lat_from_mercator(y + CELL_SIZE_M)
|
|
features.append(
|
|
{
|
|
"type": "Feature",
|
|
"properties": {
|
|
"kind": "cell_coast" if coast_hit else "cell_sea",
|
|
"fpc": fpc,
|
|
"prc": prc,
|
|
"cluster_idx": int(cluster_idx),
|
|
"part_idx": int(part_idx),
|
|
"window_idx": int(window_idx),
|
|
"cell_id": f"{int(round(x / CELL_SIZE_M))}:{int(round(y / CELL_SIZE_M))}",
|
|
"row": int(round(y / CELL_SIZE_M)),
|
|
"col": int(round(x / CELL_SIZE_M)),
|
|
"cell_size_m": CELL_SIZE_M,
|
|
"coast_hit": coast_hit,
|
|
},
|
|
"geometry": rect_lonlat(min_lon, min_lat, max_lon, max_lat),
|
|
}
|
|
)
|
|
if coast_hit:
|
|
coast_count += 1
|
|
else:
|
|
sea_count += 1
|
|
y += CELL_SIZE_M
|
|
|
|
return features, coast_count, sea_count
|
|
|
|
|
|
def main() -> None:
|
|
parser = argparse.ArgumentParser(description="Export one fish port debug bundle")
|
|
parser.add_argument("--db-name", default=DB_NAME)
|
|
parser.add_argument("--out-dir", default=DEFAULT_OUT_DIR)
|
|
parser.add_argument("--fpc", required=True, help="港口 FPC 号，例如 4320050")
|
|
parser.add_argument("--coast-buffer-m", type=float, default=DEFAULT_BUFFER_M, help="海岸缓冲距离（米）")
|
|
parser.add_argument(
|
|
"--coast-coarse-margin-m",
|
|
type=float,
|
|
default=DEFAULT_COARSE_MARGIN_M,
|
|
help="coast_200m 粗窗查询外扩距离（米）",
|
|
)
|
|
parser.add_argument("--cluster-gap-m", type=float, default=250.0, help="线聚类间隔（米）")
|
|
parser.add_argument("--max-scan-cells", type=int, default=200000, help="单窗允许扫描的最大 20m 格数")
|
|
args = parser.parse_args()
|
|
|
|
project_root = Path(__file__).resolve().parent.parent
|
|
out_dir = Path(args.out_dir)
|
|
if not out_dir.is_absolute():
|
|
out_dir = project_root / out_dir
|
|
|
|
fpc = str(args.fpc).strip()
|
|
buffer_dir = out_dir / f"fpc_{fpc}" / f"buffer_{buffer_tag(args.coast_buffer_m)}"
|
|
buffer_dir.mkdir(parents=True, exist_ok=True)
|
|
|
|
root = load_xml(project_root / FISH_SRC)
|
|
point_lookup = build_point_lookup(root)
|
|
prcs = discover_target_prcs(root, None, fpc)
|
|
if not prcs:
|
|
raise SystemExit(f"no PRC matched for FPC={fpc}")
|
|
prc = prcs[0]
|
|
|
|
fish_lines = collect_fish_lines_for_prc(root, point_lookup, prc, selected_fpc=fpc)
|
|
if not fish_lines:
|
|
raise SystemExit(f"no fish-port line features selected for FPC={fpc}")
|
|
|
|
fish_bboxes = [line.bbox for line in fish_lines]
|
|
fish_bbox = (
|
|
min(bb[0] for bb in fish_bboxes),
|
|
min(bb[1] for bb in fish_bboxes),
|
|
max(bb[2] for bb in fish_bboxes),
|
|
max(bb[3] for bb in fish_bboxes),
|
|
)
|
|
|
|
coast_lines, coast_bbox = load_coastlines_for_prc(
|
|
project_root,
|
|
prc,
|
|
bbox=(
|
|
mercator_x(fish_bbox[0]),
|
|
mercator_y(fish_bbox[1]),
|
|
mercator_x(fish_bbox[2]),
|
|
mercator_y(fish_bbox[3]),
|
|
),
|
|
search_margin_m=max(args.coast_buffer_m, args.coast_coarse_margin_m),
|
|
)
|
|
coast_geom = unary_union(coast_lines)
|
|
coast_buffer = coast_geom.buffer(args.coast_buffer_m, cap_style=2, join_style=2)
|
|
coast_prepared = prep(coast_buffer)
|
|
|
|
conn = mysql_connect(args.db_name)
|
|
try:
|
|
prc_coarse_windows = load_coarse_windows_for_bbox(
|
|
conn,
|
|
layer_name=COAST_COARSE_LAYER,
|
|
bbox_mercator=(
|
|
mercator_x(fish_bbox[0]),
|
|
mercator_y(fish_bbox[1]),
|
|
mercator_x(fish_bbox[2]),
|
|
mercator_y(fish_bbox[3]),
|
|
),
|
|
margin_m=args.coast_coarse_margin_m,
|
|
)
|
|
|
|
cluster_groups = cluster_line_indices([line.bbox for line in fish_lines], gap_m=args.cluster_gap_m)
|
|
|
|
features: list[dict] = []
|
|
features.append(
|
|
bbox_feature(
|
|
fish_bbox,
|
|
kind="fish_bbox",
|
|
fpc=fpc,
|
|
prc=prc,
|
|
label="fish_bbox",
|
|
)
|
|
)
|
|
features.append(
|
|
mercator_bbox_feature(
|
|
coast_bbox,
|
|
kind="coast_bbox",
|
|
fpc=fpc,
|
|
prc=prc,
|
|
label="coast_bbox",
|
|
)
|
|
)
|
|
for line in fish_lines:
|
|
features.append(
|
|
lonlat_line_feature(
|
|
line.coords,
|
|
kind="fish_line",
|
|
fpc=fpc,
|
|
prc=prc,
|
|
label="fish_line",
|
|
)
|
|
)
|
|
for line in coast_lines:
|
|
features.append(
|
|
mercator_line_feature(
|
|
line,
|
|
kind="coast_line",
|
|
fpc=fpc,
|
|
prc=prc,
|
|
label="coast_line",
|
|
)
|
|
)
|
|
|
|
features.append(
|
|
{
|
|
"type": "Feature",
|
|
"properties": {
|
|
"kind": "coast_buffer",
|
|
"fpc": fpc,
|
|
"prc": prc,
|
|
"buffer_m": args.coast_buffer_m,
|
|
"label": f"buffer={args.coast_buffer_m}m",
|
|
},
|
|
"geometry": mercator_geom_to_lonlat_geojson(coast_buffer),
|
|
}
|
|
)
|
|
|
|
for idx, coarse_window in enumerate(prc_coarse_windows):
|
|
features.append(
|
|
bbox_feature(
|
|
mercator_to_lonlat_bbox(coarse_window),
|
|
kind="coast_200m_window",
|
|
fpc=fpc,
|
|
prc=prc,
|
|
label=f"coast_200m_window_{idx + 1}",
|
|
)
|
|
)
|
|
|
|
scan_window_total = 0
|
|
cell_coast_total = 0
|
|
cell_sea_total = 0
|
|
|
|
for cluster_idx, cluster in enumerate(cluster_groups):
|
|
cluster_lines = [
|
|
LineString([(mercator_x(lon), mercator_y(lat)) for lon, lat in fish_lines[i].coords])
|
|
for i in cluster
|
|
]
|
|
merged = unary_union(cluster_lines)
|
|
polygons = list(polygonize(merged))
|
|
if polygons:
|
|
fish_geom = unary_union(polygons)
|
|
else:
|
|
fish_geom = merged.buffer(CELL_SIZE_M, cap_style=2, join_style=2)
|
|
|
|
fish_parts = [part for part in get_parts(fish_geom) if not part.is_empty]
|
|
for part_idx, fish_part in enumerate(fish_parts):
|
|
part_bbox = fish_part.bounds
|
|
part_coarse_windows = load_coarse_windows_for_bbox(
|
|
conn,
|
|
layer_name=COAST_COARSE_LAYER,
|
|
bbox_mercator=part_bbox,
|
|
margin_m=args.coast_coarse_margin_m,
|
|
)
|
|
for window_idx, coarse_window in enumerate(part_coarse_windows):
|
|
window_geom = box(*coarse_window)
|
|
scan_geom = window_geom
|
|
scan_window_total += 1
|
|
features.append(
|
|
{
|
|
"type": "Feature",
|
|
"properties": {
|
|
"kind": "scan_window",
|
|
"fpc": fpc,
|
|
"prc": prc,
|
|
"cluster_idx": int(cluster_idx),
|
|
"part_idx": int(part_idx),
|
|
"window_idx": int(window_idx),
|
|
"label": f"scan_window_{cluster_idx + 1}_{part_idx + 1}_{window_idx + 1}",
|
|
},
|
|
"geometry": mercator_geom_to_lonlat_geojson(scan_geom),
|
|
}
|
|
)
|
|
cells, coast_count, sea_count = scan_debug_cells(
|
|
scan_geom,
|
|
coast_prepared,
|
|
prc=prc,
|
|
fpc=fpc,
|
|
cluster_idx=cluster_idx,
|
|
part_idx=part_idx,
|
|
window_idx=window_idx,
|
|
max_scan_cells=args.max_scan_cells,
|
|
)
|
|
features.extend(cells)
|
|
cell_coast_total += coast_count
|
|
cell_sea_total += sea_count
|
|
|
|
out_path = buffer_dir / "debug.geojson"
|
|
out_path.write_text(
|
|
json.dumps({"type": "FeatureCollection", "features": features}, ensure_ascii=False, indent=2),
|
|
encoding="utf-8",
|
|
)
|
|
|
|
manifest = {
|
|
"database": args.db_name,
|
|
"generated_at": dt.datetime.now().isoformat(timespec="seconds"),
|
|
"fpc": fpc,
|
|
"prc": prc,
|
|
"buffer_m": args.coast_buffer_m,
|
|
"coast_coarse_margin_m": args.coast_coarse_margin_m,
|
|
"fish_bbox_lonlat": list(fish_bbox),
|
|
"fish_line_count": len(fish_lines),
|
|
"coast_curve_count": len(coast_lines),
|
|
"coast_bbox_mercator": list(coast_bbox),
|
|
"coast_200m_window_count": len(prc_coarse_windows),
|
|
"scan_window_count": scan_window_total,
|
|
"cell_coast_count": cell_coast_total,
|
|
"cell_sea_count": cell_sea_total,
|
|
"feature_count": len(features),
|
|
"export_file": str(out_path.relative_to(project_root)) if project_root in out_path.parents else str(out_path),
|
|
}
|
|
(buffer_dir / "manifest.json").write_text(json.dumps(manifest, ensure_ascii=False, indent=2), encoding="utf-8")
|
|
|
|
print(f"debug: {len(features)} -> {out_path}")
|
|
print(f"manifest: {buffer_dir / 'manifest.json'}")
|
|
finally:
|
|
conn.close()
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|