573 lines
20 KiB
Python
573 lines
20 KiB
Python
#!/usr/bin/env python3
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import hashlib
|
|
import json
|
|
import math
|
|
import os
|
|
import sys
|
|
from datetime import datetime, timezone
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
ROOT = Path(__file__).resolve().parents[1]
|
|
SRC = ROOT / "src"
|
|
if str(SRC) not in sys.path:
|
|
sys.path.insert(0, str(SRC))
|
|
|
|
import mapbox_vector_tile
|
|
import pymysql
|
|
import yaml
|
|
from dotenv import load_dotenv
|
|
from shapely import delaunay_triangles, line_merge, union_all, wkb
|
|
from shapely.geometry import LineString, MultiPoint, box, mapping
|
|
from shapely.ops import transform
|
|
|
|
from cm95pbf.policy import (
|
|
WEB_MERCATOR_HALF_WORLD,
|
|
lonlat_to_mercator,
|
|
mercator_to_lonlat,
|
|
xyz_tile_bounds,
|
|
)
|
|
|
|
DEFAULT_CONFIG = ROOT / "config" / "bathymetry" / "cm95_bathymetry_north_kyushu_v1.yaml"
|
|
|
|
|
|
def parse_args() -> argparse.Namespace:
|
|
p = argparse.ArgumentParser(
|
|
description="Build an independent low-reliability bathymetry contour PBF from CM95 sounding points."
|
|
)
|
|
p.add_argument("--config", type=Path, default=DEFAULT_CONFIG)
|
|
p.add_argument("--output-dir", type=Path, default=None)
|
|
return p.parse_args()
|
|
|
|
|
|
def load_config(path: Path) -> dict[str, Any]:
|
|
obj = yaml.safe_load(path.read_text(encoding="utf-8"))
|
|
if not isinstance(obj, dict):
|
|
raise RuntimeError(f"invalid config: {path}")
|
|
return obj
|
|
|
|
|
|
def connect_db():
|
|
load_dotenv(ROOT / ".env")
|
|
return pymysql.connect(
|
|
host=os.getenv("CM93_DB_HOST", "127.0.0.1"),
|
|
port=int(os.getenv("CM93_DB_PORT", "3306")),
|
|
user=os.getenv("CM93_DB_USER", "root"),
|
|
password=os.getenv("CM93_DB_PASSWORD"),
|
|
database=os.getenv("CM93_DB_NAME", "cm93_chart"),
|
|
charset=os.getenv("CM93_DB_CHARSET", "utf8mb4"),
|
|
cursorclass=pymysql.cursors.DictCursor,
|
|
autocommit=False,
|
|
connect_timeout=20,
|
|
)
|
|
|
|
def envelope_wkt(bbox_wgs84: list[float]) -> str:
|
|
west, south, east, north = [float(v) for v in bbox_wgs84]
|
|
return (
|
|
"POLYGON(("
|
|
f"{west:.17g} {south:.17g},"
|
|
f"{east:.17g} {south:.17g},"
|
|
f"{east:.17g} {north:.17g},"
|
|
f"{west:.17g} {north:.17g},"
|
|
f"{west:.17g} {south:.17g}"
|
|
"))"
|
|
)
|
|
|
|
|
|
def coord_key(x: float, y: float) -> tuple[float, float]:
|
|
return round(float(x), 3), round(float(y), 3)
|
|
|
|
|
|
def fetch_soundings(conn, cfg: dict[str, Any]) -> tuple[dict[tuple[float, float], float], dict[str, int]]:
|
|
bbox = cfg["bbox_wgs84"]
|
|
native_scale_max = int(cfg["native_scale_max"])
|
|
sql = """
|
|
SELECT ST_X(s.geom) AS lon, ST_Y(s.geom) AS lat, s.depth, f.native_scale
|
|
FROM chart_sounding_point s FORCE INDEX (spx_sounding_geom)
|
|
STRAIGHT_JOIN chart_feature f ON f.id=s.feature_id
|
|
WHERE MBRContains(ST_GeomFromText(%s, 4326), s.geom)
|
|
AND f.dataset_id=%s
|
|
AND f.native_scale <= %s
|
|
"""
|
|
with conn.cursor() as cur:
|
|
cur.execute(sql, (envelope_wkt(bbox), int(cfg["dataset_id"]), native_scale_max))
|
|
rows = cur.fetchall()
|
|
|
|
points: dict[tuple[float, float], float] = {}
|
|
|
|
duplicate_count = 0
|
|
for row in rows:
|
|
depth = float(row["depth"])
|
|
if not math.isfinite(depth) or depth < 0:
|
|
continue
|
|
mx, my = lonlat_to_mercator(float(row["lon"]), float(row["lat"]))
|
|
key = coord_key(mx, my)
|
|
previous = points.get(key)
|
|
if previous is not None:
|
|
duplicate_count += 1
|
|
if depth < previous:
|
|
points[key] = depth
|
|
else:
|
|
points[key] = depth
|
|
|
|
if len(points) < 3:
|
|
raise RuntimeError("fewer than three usable sounding points")
|
|
return points, {
|
|
"raw_sounding_rows": len(rows),
|
|
"unique_sounding_points": len(points),
|
|
"duplicate_sounding_rows": duplicate_count,
|
|
}
|
|
|
|
|
|
def fetch_water_mask(conn, cfg: dict[str, Any]):
|
|
classes = [str(v) for v in cfg.get("water_object_classes", ["DEPARE"])]
|
|
if not classes:
|
|
raise RuntimeError("water_object_classes must not be empty")
|
|
placeholders = ",".join(["%s"] * len(classes))
|
|
sql = f"""
|
|
SELECT ST_AsWKB(g.geom) AS geom_wkb
|
|
FROM chart_area_geometry g FORCE INDEX (spx_area_geom)
|
|
STRAIGHT_JOIN chart_feature f ON f.id=g.feature_id
|
|
WHERE MBRIntersects(g.geom, ST_GeomFromText(%s, 4326))
|
|
AND f.dataset_id=%s
|
|
AND f.native_scale <= %s
|
|
AND f.object_class IN ({placeholders})
|
|
"""
|
|
|
|
params: list[Any] = [
|
|
envelope_wkt(cfg["bbox_wgs84"]),
|
|
int(cfg["dataset_id"]),
|
|
int(cfg["native_scale_max"]),
|
|
*classes,
|
|
]
|
|
with conn.cursor() as cur:
|
|
cur.execute(sql, params)
|
|
rows = cur.fetchall()
|
|
|
|
west, south, east, north = [float(v) for v in cfg["bbox_wgs84"]]
|
|
west_m, south_m = lonlat_to_mercator(west, south)
|
|
east_m, north_m = lonlat_to_mercator(east, north)
|
|
region_box = box(west_m, south_m, east_m, north_m)
|
|
|
|
geoms = []
|
|
for row in rows:
|
|
geom = wkb.loads(bytes(row["geom_wkb"]))
|
|
if geom.is_empty:
|
|
continue
|
|
geom_m = transform(lambda x, y, z=None: lonlat_to_mercator(x, y), geom)
|
|
clipped = geom_m.intersection(region_box)
|
|
if not clipped.is_empty:
|
|
geoms.append(clipped)
|
|
|
|
if not geoms:
|
|
raise RuntimeError("no DEPARE/DRGARE water geometry found in build bbox")
|
|
return union_all(geoms), len(rows)
|
|
|
|
|
|
def edge_length(a: tuple[float, float], b: tuple[float, float]) -> float:
|
|
return math.hypot(a[0] - b[0], a[1] - b[1])
|
|
|
|
def contour_segment(
|
|
coords: list[tuple[float, float]],
|
|
depths: list[float],
|
|
level: float,
|
|
) -> LineString | None:
|
|
eps = 1e-9
|
|
hits: list[tuple[float, float]] = []
|
|
for i, j in ((0, 1), (1, 2), (2, 0)):
|
|
p1, p2 = coords[i], coords[j]
|
|
d1, d2 = depths[i], depths[j]
|
|
if abs(d1 - level) <= eps and abs(d2 - level) <= eps:
|
|
hits.extend((p1, p2))
|
|
continue
|
|
if abs(d1 - d2) <= eps:
|
|
continue
|
|
if (level - d1) * (level - d2) > eps:
|
|
continue
|
|
t = (level - d1) / (d2 - d1)
|
|
if -eps <= t <= 1.0 + eps:
|
|
t = max(0.0, min(1.0, t))
|
|
hits.append((p1[0] + t * (p2[0] - p1[0]), p1[1] + t * (p2[1] - p1[1])))
|
|
|
|
unique: list[tuple[float, float]] = []
|
|
seen: set[tuple[float, float]] = set()
|
|
for hit in hits:
|
|
key = coord_key(*hit)
|
|
if key not in seen:
|
|
seen.add(key)
|
|
unique.append(hit)
|
|
if len(unique) < 2:
|
|
return None
|
|
|
|
best_pair = None
|
|
best_distance = -1.0
|
|
for i in range(len(unique)):
|
|
for j in range(i + 1, len(unique)):
|
|
distance = edge_length(unique[i], unique[j])
|
|
if distance > best_distance:
|
|
|
|
best_distance = distance
|
|
best_pair = (unique[i], unique[j])
|
|
if best_pair is None or best_distance <= 0:
|
|
return None
|
|
return LineString(best_pair)
|
|
|
|
|
|
def line_parts(geom) -> list[LineString]:
|
|
if geom is None or geom.is_empty:
|
|
return []
|
|
if geom.geom_type == "LineString":
|
|
return [geom]
|
|
if geom.geom_type == "MultiLineString":
|
|
return list(geom.geoms)
|
|
out: list[LineString] = []
|
|
for part in getattr(geom, "geoms", []):
|
|
out.extend(line_parts(part))
|
|
return out
|
|
|
|
|
|
def extract_contours(
|
|
points: dict[tuple[float, float], float],
|
|
water_mask,
|
|
cfg: dict[str, Any],
|
|
) -> tuple[dict[float, Any], dict[str, Any]]:
|
|
levels = sorted(float(v) for v in cfg["contour_levels_m"])
|
|
max_edge = float(cfg["max_triangle_edge_m"])
|
|
min_segment = float(cfg["min_segment_length_m"])
|
|
regional_simplify = float(cfg.get("regional_simplify_m", 0))
|
|
depth_by_coord = dict(points)
|
|
|
|
triangles = delaunay_triangles(MultiPoint(list(points.keys())))
|
|
segments: dict[float, list[LineString]] = {level: [] for level in levels}
|
|
triangle_total = 0
|
|
triangle_accepted = 0
|
|
|
|
for triangle in triangles.geoms:
|
|
triangle_total += 1
|
|
coords = [(float(x), float(y)) for x, y, *rest in list(triangle.exterior.coords)[:3]]
|
|
try:
|
|
depths = [depth_by_coord[coord_key(*coord)] for coord in coords]
|
|
except KeyError:
|
|
continue
|
|
longest = max(
|
|
edge_length(coords[0], coords[1]),
|
|
edge_length(coords[1], coords[2]),
|
|
edge_length(coords[2], coords[0]),
|
|
)
|
|
if longest > max_edge:
|
|
continue
|
|
triangle_accepted += 1
|
|
lo, hi = min(depths), max(depths)
|
|
for level in levels:
|
|
if level < lo or level > hi:
|
|
continue
|
|
segment = contour_segment(coords, depths, level)
|
|
if segment is not None and segment.length >= min_segment:
|
|
segments[level].append(segment)
|
|
|
|
contours: dict[float, Any] = {}
|
|
per_level: dict[str, Any] = {}
|
|
for level in levels:
|
|
raw = segments[level]
|
|
if not raw:
|
|
per_level[str(level)] = {"raw_segments": 0, "line_parts": 0, "length_m": 0.0}
|
|
continue
|
|
merged = line_merge(union_all(raw))
|
|
clipped = merged.intersection(water_mask)
|
|
parts = [part for part in line_parts(clipped) if part.length >= min_segment]
|
|
if regional_simplify > 0:
|
|
parts = [
|
|
part.simplify(regional_simplify, preserve_topology=False)
|
|
for part in parts
|
|
if not part.is_empty
|
|
]
|
|
|
|
parts = [part for part in parts if not part.is_empty and part.length >= min_segment]
|
|
if not parts:
|
|
per_level[str(level)] = {
|
|
"raw_segments": len(raw),
|
|
"line_parts": 0,
|
|
"length_m": 0.0,
|
|
}
|
|
continue
|
|
geom = union_all(parts)
|
|
contours[level] = geom
|
|
per_level[str(level)] = {
|
|
"raw_segments": len(raw),
|
|
"line_parts": len(line_parts(geom)),
|
|
"length_m": round(sum(part.length for part in line_parts(geom)), 3),
|
|
}
|
|
|
|
return contours, {
|
|
"triangle_total": triangle_total,
|
|
"triangle_accepted": triangle_accepted,
|
|
"triangle_rejected_long_edge": triangle_total - triangle_accepted,
|
|
"levels": per_level,
|
|
}
|
|
|
|
|
|
def mercator_tile_xy(mx: float, my: float, z: int) -> tuple[int, int]:
|
|
n = 1 << int(z)
|
|
world = WEB_MERCATOR_HALF_WORLD * 2.0
|
|
x = int(math.floor((mx + WEB_MERCATOR_HALF_WORLD) / world * n))
|
|
y = int(math.floor((WEB_MERCATOR_HALF_WORLD - my) / world * n))
|
|
return max(0, min(n - 1, x)), max(0, min(n - 1, y))
|
|
|
|
|
|
def tile_range(bbox_wgs84: list[float], z: int) -> tuple[range, range]:
|
|
west, south, east, north = [float(v) for v in bbox_wgs84]
|
|
west_m, south_m = lonlat_to_mercator(west, south)
|
|
|
|
east_m, north_m = lonlat_to_mercator(east, north)
|
|
x0, y1 = mercator_tile_xy(west_m, south_m, z)
|
|
x1, y0 = mercator_tile_xy(east_m, north_m, z)
|
|
return range(min(x0, x1), max(x0, x1) + 1), range(min(y0, y1), max(y0, y1) + 1)
|
|
|
|
|
|
def line_survives_mvt_quantization(geom: LineString, bounds, extent: int) -> bool:
|
|
minx, miny, maxx, maxy = bounds
|
|
sx = float(extent) / float(maxx - minx)
|
|
sy = float(extent) / float(maxy - miny)
|
|
quantized = {
|
|
(
|
|
int((float(x) - float(minx)) * sx),
|
|
int((float(y) - float(miny)) * sy),
|
|
)
|
|
for x, y, *rest in geom.coords
|
|
}
|
|
return len(quantized) >= 2
|
|
|
|
|
|
def feature_id(build_id: str, z: int, x: int, y: int, depth: float, index: int, geom) -> int:
|
|
payload = f"{build_id}|{z}|{x}|{y}|{depth:.6f}|{index}|{geom.wkb_hex}".encode("ascii")
|
|
value = int.from_bytes(hashlib.blake2b(payload, digest_size=8).digest(), "big") & ((1 << 63) - 1)
|
|
return value or 1
|
|
|
|
|
|
def build_tiles(contours: dict[float, Any], cfg: dict[str, Any], output_dir: Path) -> dict[str, Any]:
|
|
pbf_root = output_dir / "pbf"
|
|
pbf_root.mkdir(parents=True, exist_ok=True)
|
|
source_layer = str(cfg["source_layer"])
|
|
extent = int(cfg["extent"])
|
|
buffer_px = float(cfg.get("tile_buffer_px", 8))
|
|
build_id = str(cfg["build_id"])
|
|
simplify_by_zoom = {int(k): float(v) for k, v in cfg.get("simplify_m_by_zoom", {}).items()}
|
|
|
|
tile_count = 0
|
|
feature_count = 0
|
|
total_bytes = 0
|
|
encoder_dropped = 0
|
|
zoom_stats: dict[str, Any] = {}
|
|
|
|
for z in range(int(cfg["minzoom"]), int(cfg["maxzoom"]) + 1):
|
|
z_tiles = 0
|
|
z_features = 0
|
|
z_bytes = 0
|
|
z_encoder_dropped = 0
|
|
xr, yr = tile_range(cfg["bbox_wgs84"], z)
|
|
for x in xr:
|
|
for y in yr:
|
|
|
|
tile = xyz_tile_bounds(z, x, y)
|
|
pad = tile.width_m * buffer_px / 256.0
|
|
clip_box = box(
|
|
tile.west_m - pad,
|
|
tile.south_m - pad,
|
|
tile.east_m + pad,
|
|
tile.north_m + pad,
|
|
)
|
|
features = []
|
|
simplify_m = simplify_by_zoom.get(z, 0.0)
|
|
local_index = 0
|
|
|
|
for depth, contour_geom in sorted(contours.items()):
|
|
if contour_geom.is_empty or not contour_geom.intersects(clip_box):
|
|
continue
|
|
clipped = contour_geom.intersection(clip_box)
|
|
for part in line_parts(clipped):
|
|
if simplify_m > 0:
|
|
part = part.simplify(simplify_m, preserve_topology=False)
|
|
if part.is_empty or part.length <= 0:
|
|
continue
|
|
if not line_survives_mvt_quantization(part, tile.mercator_tuple, extent):
|
|
continue
|
|
local_index += 1
|
|
features.append(
|
|
{
|
|
"id": feature_id(build_id, z, x, y, depth, local_index, part),
|
|
"geometry": mapping(part),
|
|
"properties": {
|
|
"depth_m": float(depth),
|
|
"source_kind": str(cfg["source_kind"]),
|
|
"reliability": str(cfg["reliability"]),
|
|
"build_id": build_id,
|
|
},
|
|
}
|
|
)
|
|
|
|
if not features:
|
|
continue
|
|
|
|
data = mapbox_vector_tile.encode(
|
|
[{"name": source_layer, "features": features}],
|
|
default_options={
|
|
"quantize_bounds": tile.mercator_tuple,
|
|
"extents": extent,
|
|
},
|
|
)
|
|
decoded = mapbox_vector_tile.decode(data)
|
|
actual = len(decoded.get(source_layer, {}).get("features", []))
|
|
dropped = max(0, len(features) - actual)
|
|
encoder_dropped += dropped
|
|
z_encoder_dropped += dropped
|
|
|
|
tile_path = pbf_root / str(z) / str(x) / f"{y}.pbf"
|
|
tile_path.parent.mkdir(parents=True, exist_ok=True)
|
|
tile_path.write_bytes(data)
|
|
tile_count += 1
|
|
feature_count += actual
|
|
total_bytes += len(data)
|
|
z_tiles += 1
|
|
z_features += actual
|
|
z_bytes += len(data)
|
|
|
|
zoom_stats[str(z)] = {
|
|
"tiles": z_tiles,
|
|
"features": z_features,
|
|
"encoder_dropped_features": z_encoder_dropped,
|
|
"bytes": z_bytes,
|
|
}
|
|
|
|
return {
|
|
"pbf_tiles": tile_count,
|
|
"tile_feature_instances": feature_count,
|
|
"encoder_dropped_features": encoder_dropped,
|
|
"pbf_bytes": total_bytes,
|
|
"zooms": zoom_stats,
|
|
}
|
|
|
|
def write_debug_geojson(contours: dict[float, Any], cfg: dict[str, Any], output_dir: Path) -> int:
|
|
features = []
|
|
for depth, geom in sorted(contours.items()):
|
|
geom_wgs = transform(lambda x, y, z=None: mercator_to_lonlat(x, y), geom)
|
|
features.append(
|
|
{
|
|
"type": "Feature",
|
|
"geometry": mapping(geom_wgs),
|
|
"properties": {
|
|
"depth_m": float(depth),
|
|
"source_kind": str(cfg["source_kind"]),
|
|
"reliability": str(cfg["reliability"]),
|
|
"build_id": str(cfg["build_id"]),
|
|
},
|
|
}
|
|
)
|
|
payload = {"type": "FeatureCollection", "features": features}
|
|
(output_dir / "contours.geojson").write_text(
|
|
json.dumps(payload, ensure_ascii=False, separators=(",", ":")),
|
|
encoding="utf-8",
|
|
)
|
|
return len(features)
|
|
|
|
|
|
def write_tilejson(cfg: dict[str, Any], output_dir: Path) -> None:
|
|
payload = {
|
|
"tilejson": "3.0.0",
|
|
"name": str(cfg["build_id"]),
|
|
"bounds": [float(v) for v in cfg["bbox_wgs84"]],
|
|
"minzoom": int(cfg["minzoom"]),
|
|
"maxzoom": int(cfg["maxzoom"]),
|
|
"tiles": ["pbf/{z}/{x}/{y}.pbf"],
|
|
"vector_layers": [
|
|
{
|
|
"id": str(cfg["source_layer"]),
|
|
|
|
"description": "Derived reference bathymetry contours generated from CM95 SOUNDG points.",
|
|
"minzoom": int(cfg["minzoom"]),
|
|
"maxzoom": int(cfg["maxzoom"]),
|
|
"fields": {
|
|
"depth_m": "Number",
|
|
"source_kind": "String",
|
|
"reliability": "String",
|
|
"build_id": "String",
|
|
},
|
|
}
|
|
],
|
|
}
|
|
(output_dir / "tilejson.json").write_text(
|
|
json.dumps(payload, ensure_ascii=False, indent=2),
|
|
encoding="utf-8",
|
|
)
|
|
|
|
|
|
def main() -> int:
|
|
args = parse_args()
|
|
config_path = args.config.resolve()
|
|
cfg = load_config(config_path)
|
|
output_dir = (
|
|
args.output_dir.resolve()
|
|
if args.output_dir is not None
|
|
else (ROOT / str(cfg["output_dir"])).resolve()
|
|
)
|
|
output_dir.mkdir(parents=True, exist_ok=True)
|
|
|
|
started = datetime.now(timezone.utc)
|
|
print(f"BUILD {cfg['build_id']}")
|
|
print(f"BBOX {cfg['bbox_wgs84']}")
|
|
print(f"OUTPUT {output_dir}")
|
|
|
|
conn = connect_db()
|
|
try:
|
|
points, sounding_stats = fetch_soundings(conn, cfg)
|
|
print("SOUNDINGS", json.dumps(sounding_stats, ensure_ascii=False))
|
|
water_mask, water_rows = fetch_water_mask(conn, cfg)
|
|
print(f"WATER_AREAS {water_rows}")
|
|
finally:
|
|
conn.close()
|
|
|
|
contours, contour_stats = extract_contours(points, water_mask, cfg)
|
|
if not contours:
|
|
raise RuntimeError("no contours generated")
|
|
print("CONTOURS", json.dumps(contour_stats, ensure_ascii=False))
|
|
|
|
debug_feature_count = write_debug_geojson(contours, cfg, output_dir)
|
|
write_tilejson(cfg, output_dir)
|
|
tile_stats = build_tiles(contours, cfg, output_dir)
|
|
print("TILES", json.dumps(tile_stats, ensure_ascii=False))
|
|
|
|
finished = datetime.now(timezone.utc)
|
|
manifest = {
|
|
"schema": "navsea.cm95.bathymetry-build-result/v1",
|
|
"build_id": str(cfg["build_id"]),
|
|
"generated_at": finished.isoformat(),
|
|
"duration_seconds": round((finished - started).total_seconds(), 3),
|
|
"dataset_id": int(cfg["dataset_id"]),
|
|
"bbox_wgs84": [float(v) for v in cfg["bbox_wgs84"]],
|
|
"source_layer": str(cfg["source_layer"]),
|
|
"source_kind": str(cfg["source_kind"]),
|
|
"reliability": str(cfg["reliability"]),
|
|
"contour_levels_m": [float(v) for v in cfg["contour_levels_m"]],
|
|
"native_scale_max": int(cfg["native_scale_max"]),
|
|
"max_triangle_edge_m": float(cfg["max_triangle_edge_m"]),
|
|
"water_object_classes": [str(v) for v in cfg["water_object_classes"]],
|
|
"sounding_stats": sounding_stats,
|
|
"water_area_rows": water_rows,
|
|
"contour_stats": contour_stats,
|
|
"debug_geojson_features": debug_feature_count,
|
|
"tile_stats": tile_stats,
|
|
}
|
|
(output_dir / "build_manifest.json").write_text(
|
|
|
|
json.dumps(manifest, ensure_ascii=False, indent=2),
|
|
encoding="utf-8",
|
|
)
|
|
print("DONE", output_dir / "build_manifest.json")
|
|
return 0
|
|
|
|
|
|
if __name__ == "__main__":
|
|
raise SystemExit(main())
|