diff --git a/AGENTS.md b/AGENTS.md index 8704ecd..43e10e4 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -18,7 +18,7 @@ Geometry Audit: OVERALL: PASS Non-Geometry Audit: OVERALL: PASS ``` -PBF 当前继续冻结。 +Importer 冻结基线继续保持不变；PBF 已在后续下游阶段由用户明确启动，PBF 工作不得反向修改冻结 importer。 --- @@ -343,12 +343,14 @@ cm93-mysql-importer-baseline-v1 当前： ```text -PBF = OUT OF SCOPE / FROZEN +Importer scope: PBF = OUT OF SCOPE +Downstream PBF stage: ACTIVE ``` -除非用户明确启动 PBF 阶段： +用户已明确启动独立 PBF 下游阶段，因此允许在 `src/cm95pbf/`、`tools/`、`config/`、`web/` 和独立 build 输出中继续实现 PBF 能力。 -- 不新增 PBF schema； -- 不新增 PBF exporter； -- 不修改 importer 为 PBF 服务； -- 不以 PBF 需求反向破坏当前 MySQL/raw 数据模型。 +仍然禁止： + +- 为 PBF 修改冻结的 parser / geometry / importer 主流程； +- 为 PBF 改坏当前 MySQL/raw 数据模型； +- 把 PBF 派生结果回写覆盖原始 CM93/CM95 证据。 diff --git a/README.md b/README.md index 49113c0..5d27315 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,7 @@ CM93 source → audit / traceability ``` -PBF 当前继续冻结，不属于本阶段范围。 +本文 importer 冻结基线本身不因下游 PBF 工作而改变；PBF 已在后续独立阶段单独启动，见本文末尾。 --- @@ -464,4 +464,23 @@ tag 应保持不可变。 - 为了下游显示而损失 raw CM93 数据； - 对原始 CM93 文件进行写回。 -PBF 只有在单独明确启动下一阶段后才进入设计和实现。 +该限制针对 importer 冻结阶段；后续 PBF 阶段已由用户明确启动，但 PBF 实现不得反向破坏已冻结的 importer / MySQL / raw 数据模型。 + + +--- + +## 13. 后续独立 PBF 阶段 + +CM95 → NavSea 的 PBF 下游阶段已经单独启动；importer 冻结基线继续保持不变。 + +2026-10-01 新增独立 bathymetry PBF 首批 build： + +```text +config/bathymetry/cm95_bathymetry_north_kyushu_v1.yaml +tools/cm95_bathymetry_build_v1.py +web/bathymetry_test.html +web/bathymetry_test_server.py +docs/CM95_Bathymetry独立PBF_v1.md +``` + +该 bathymetry PBF 是第二个独立 vector source，等深线由 SOUNDG 测深点派生，标记为 `reliability=reference_low`，不与基础 PBF 合并。 diff --git a/config/bathymetry/cm95_bathymetry_north_kyushu_v1.yaml b/config/bathymetry/cm95_bathymetry_north_kyushu_v1.yaml new file mode 100644 index 0000000..446617f --- /dev/null +++ b/config/bathymetry/cm95_bathymetry_north_kyushu_v1.yaml @@ -0,0 +1,24 @@ +schema: navsea.cm95.bathymetry-build/v1 +build_id: cm95-bathymetry-north-kyushu-v1 +dataset_id: 1 +bbox_wgs84: [129.7, 33.1, 130.6, 33.8] +source_layer: bathymetry_contour +source_kind: derived_sounding_tin +reliability: reference_low +native_scale_max: 100000 +water_object_classes: [DEPARE, DRGARE] +contour_levels_m: [2, 5, 10, 20, 30, 50] +max_triangle_edge_m: 4000 +min_segment_length_m: 20 +regional_simplify_m: 2 +minzoom: 8 +maxzoom: 12 +extent: 4096 +tile_buffer_px: 8 +simplify_m_by_zoom: + "8": 80 + "9": 40 + "10": 20 + "11": 10 + "12": 5 +output_dir: build/bathymetry/cm95-bathymetry-north-kyushu-v1 diff --git a/docs/CM95_Bathymetry独立PBF_v1.md b/docs/CM95_Bathymetry独立PBF_v1.md new file mode 100644 index 0000000..7341356 --- /dev/null +++ b/docs/CM95_Bathymetry独立PBF_v1.md @@ -0,0 +1,199 @@ +# CM95 Bathymetry 独立 PBF v1 + +日期：2026-10-01 + +## 1. 目的 + +本阶段新增一套与现有 CM95 基础 PBF 完全独立的 bathymetry PBF。 + +客户端需要显示等深线时增加第二个 vector source；关闭等深线时可以不加载该 source。 +现有基础 PBF 的 schema、tile 和发布生命周期不因本 build 改变。 + +本 PBF 中的等深线属于 NavSea/本项目自行生产的派生参考数据，不视为 CM95 原始等深线、官方水路资料或安全航行依据。 + +## 2. 首批 build + +Build ID： + +```text +cm95-bathymetry-north-kyushu-v1 +``` + +范围： + +```text +W 129.7 +S 33.1 +E 130.6 +N 33.8 +``` + +对应北九州、唐津—福冈附近的首批验证区域。 +生成配置： + +```text +config/bathymetry/cm95_bathymetry_north_kyushu_v1.yaml +``` + +生成器： + +```text +tools/cm95_bathymetry_build_v1.py +``` + +本地生成目录（build/ 被 Git 忽略）： + +```text +build/bathymetry/cm95-bathymetry-north-kyushu-v1/ +``` + +目录内包含： + +```text +build_manifest.json +tilejson.json +contours.geojson +pbf/{z}/{x}/{y}.pbf +``` + +## 3. 输入与生成规则 + +首批 build 只使用 Dataset ID=1 的 SOUNDG sounding points 作为连续水深面的数值输入。 + +为避免低分辨率点拉坏近岸插值，当前只使用： + +```text +native_scale <= 100000 +``` +同坐标存在重复 sounding 时取较浅值。 + +生成流程： + +```text +SOUNDG sounding points + -> Web Mercator + -> Delaunay / TIN + -> 丢弃最大边 > 4000m 的三角形 + -> 三角形内线性等值线 + -> DEPARE + DRGARE 水域 union 裁剪 + -> 区域级合并/简化 + -> 按 XYZ tile 切片 + -> MVT/PBF +``` + +首批 contour levels： + +```text +2m / 5m / 10m / 20m / 30m / 50m +``` + +先生成区域连续线，再切 tile；不得每个 tile 独立插值，以避免 tile seam。 + +## 4. PBF contract + +独立 source-layer： + +```text +bathymetry_contour +``` + +每个 LineString feature 至少包含： + +```text +depth_m +source_kind = derived_sounding_tin +reliability = reference_low +build_id = cm95-bathymetry-north-kyushu-v1 +``` +首批 zoom： + +```text +z8-z12 +extent = 4096 +``` + +粗 zoom 中亚像素短线可能被 MVT encoder 丢弃。生成器不把这种情况静默处理： +实际 drop 数写入 `build_manifest.json`，并以 decode 后实际 feature 数作为 build 真值。 + +## 5. 首批结果 + +输入： + +```text +raw sounding rows 10035 +unique sounding points 10028 +duplicate rows 7 +water area rows 1884 +``` + +TIN： + +```text +triangles total 20028 +accepted 19773 +rejected long edge 255 +``` + +PBF： + +```text +tiles 93 +actual feature instances 15529 +encoder dropped 76 +total bytes 612694 +``` + +其中 encoder drop 主要集中在 z8 的亚像素短线。 +全量解码审计结果： + +```text +tiles 93 / manifest 93 +features 15529 / manifest 15529 +geometry types: LineString 15529 +bad layers: 0 +bad depth values: 0 +bad properties: 0 +AUDIT PASS +``` + +## 6. 本地测试 + +测试页： + +```text +web/bathymetry_test.html +``` + +测试 server： + +```text +web/bathymetry_test_server.py +``` + +启动： + +```powershell +python web\bathymetry_test_server.py --host 0.0.0.0 --port 18890 +``` + +浏览： + +```text +http://127.0.0.1:18890/ +``` + +测试页使用在线 GSI raster 作为普通背景，只加载这一套新的 bathymetry vector source，因此不依赖现有 CM95 大 build。 +## 7. 当前边界 + +本 v1 只验证“独立 bathymetry PBF”生产链可行，不代表全国参数已经冻结。 + +当前参数仍需通过实际海图目视检查决定是否调整，重点包括： + +- 4000m 最大 TIN 边长是否过大或过小； +- 2/5/10/20/30/50m 层级是否合适； +- 岛屿、港池和狭窄水道附近是否存在不合理跨越； +- DEPARE/DRGARE 水域 mask 是否存在局部缺口； +- z8-z12 的简化强度是否合适； +- 是否需要在后续 build 引入额外约束数据。 + +在这些项目通过人工检查前，`reliability=reference_low` 必须保持不变。 diff --git a/tools/cm95_bathymetry_build_v1.py b/tools/cm95_bathymetry_build_v1.py new file mode 100644 index 0000000..390f229 --- /dev/null +++ b/tools/cm95_bathymetry_build_v1.py @@ -0,0 +1,572 @@ +#!/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()) diff --git a/web/bathymetry_test.html b/web/bathymetry_test.html new file mode 100644 index 0000000..c32b9f8 --- /dev/null +++ b/web/bathymetry_test.html @@ -0,0 +1,123 @@ + + + + + +CM95 派生等深线独立 PBF 测试 + + + + +
+
+

独立 bathymetry PBF

+ + +
单独 vector source：bathymetry；source-layer：bathymetry_contour。
+
底图仅用 GSI raster 便于检查，不依赖现有 CM95 大 build。
+
等深线由 SOUNDG 测深点经 TIN 插值、水域裁剪生成，仅作低可靠性参考。
+
2m5m10m20m30m50m
+
加载 build manifest…
+
+ + + + diff --git a/web/bathymetry_test_server.py b/web/bathymetry_test_server.py new file mode 100644 index 0000000..2e371e5 --- /dev/null +++ b/web/bathymetry_test_server.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import json +from http.server import BaseHTTPRequestHandler, ThreadingHTTPServer +from pathlib import Path +from urllib.parse import urlparse + +BATHY_BUILD_ID = "cm95-bathymetry-north-kyushu-v1" + + +def parse_args(): + p = argparse.ArgumentParser(description="CM95 independent bathymetry PBF test server") + p.add_argument("--host", default="0.0.0.0") + p.add_argument("--port", type=int, default=18890) + p.add_argument("--bathymetry-build-id", default=BATHY_BUILD_ID) + return p.parse_args() + + +def main(): + args = parse_args() + web_root = Path(__file__).resolve().parent + project_root = web_root.parent + html_path = web_root / "bathymetry_test.html" + bathy_root = project_root / "build" / "bathymetry" / args.bathymetry_build_id + bathy_pbf_root = bathy_root / "pbf" + bathy_manifest = bathy_root / "build_manifest.json" + + required = [html_path, bathy_pbf_root, bathy_manifest] + missing = [str(path) for path in required if not path.exists()] + if missing: + raise SystemExit("Missing required path(s):\n" + "\n".join(missing)) + manifest_obj = json.loads(bathy_manifest.read_text(encoding="utf-8")) + + class Handler(BaseHTTPRequestHandler): + def send_bytes(self, data: bytes, content_type: str, status: int = 200): + self.send_response(status) + self.send_header("Content-Type", content_type) + self.send_header("Content-Length", str(len(data))) + self.send_header("Cache-Control", "no-store") + self.send_header("Access-Control-Allow-Origin", "*") + self.end_headers() + self.wfile.write(data) + + def do_GET(self): + path = urlparse(self.path).path + if path in ("/", "/index.html"): + self.send_bytes(html_path.read_bytes(), "text/html; charset=utf-8") + return + if path == "/bathymetry-manifest.json": + self.send_bytes( + json.dumps(manifest_obj, ensure_ascii=False, indent=2).encode("utf-8"), + "application/json; charset=utf-8", + ) + return + if path.startswith("/bathymetry/") and path.endswith(".pbf"): + relative = path[len("/bathymetry/"):] + parts = relative.split("/") + if len(parts) != 3: + self.send_error(404) + return + z, x, yfile = parts + if not (z.isdigit() and x.isdigit() and yfile.endswith(".pbf")): + self.send_error(404) + return + y = yfile[:-4] + if not y.isdigit(): + self.send_error(404) + return + tile_path = bathy_pbf_root / z / x / yfile + if not tile_path.is_file(): + self.send_bytes(b"", "application/vnd.mapbox-vector-tile", status=204) + return + self.send_bytes(tile_path.read_bytes(), "application/vnd.mapbox-vector-tile") + return + self.send_error(404) + def log_message(self, fmt, *args): + print(f"{self.address_string()} - {fmt % args}") + + server = ThreadingHTTPServer((args.host, args.port), Handler) + print(f"Bathymetry test: http://{args.host}:{args.port}/") + print(f"Bathymetry PBF: {bathy_pbf_root}") + try: + server.serve_forever() + except KeyboardInterrupt: + pass + finally: + server.server_close() + + +if __name__ == "__main__": + main()