From 0456c443deea73f942a5b3efcead97a0b407de8d Mon Sep 17 00:00:00 2001 From: OpenAI Codex Date: Sat, 17 Oct 2026 22:27:23 +0800 Subject: [PATCH] feat: checkpoint current CM95 builder and bathymetry work --- .env.linux.example | 9 + .gitignore | 5 +- README.md | 2 + .../cm95_bathymetry_north_kyushu_v1.yaml | 24 + ...m95_bathymetry_north_kyushu_v2_smooth.yaml | 30 + ...cm95_bathymetry_north_kyushu_v3_dense.yaml | 34 + ..._bathymetry_north_kyushu_v4_tin_dense.yaml | 29 + ...5_bathymetry_north_kyushu_v5_adaptive.yaml | 37 + docs/CM95_Bathymetry独立PBF_v1.md | 199 ++++++ ...FAIL根因归并与修复批次规划_Codex任务_v1.md | 417 +++++++++++ ...与Japan_PBF_Canonical归一化缺口修复上下文_v2.md | 558 +++++++++++++++ ...划暂停与Japan_PBF低Zoom航标修复上下文_v1.md | 655 ++++++++++++++++++ scripts/setup_linux.sh | 51 ++ src/cm93db/audit.py | 15 +- src/cm93db/cli.py | 10 +- src/cm93db/database.py | 68 +- src/cm93db/importer.py | 22 +- src/cm95pbf/builder.py | 332 ++++++++- src/cm95pbf/production.py | 5 +- tools/cm93_mysql_final_geometry_audit.py | 35 +- tools/cm95_bathymetry_build_national_v5.py | 528 ++++++++++++++ tools/cm95_bathymetry_build_v1.py | 603 ++++++++++++++++ tools/cm95_bathymetry_build_v2_smooth.py | 502 ++++++++++++++ tools/cm95_bathymetry_build_v3_dense.py | 515 ++++++++++++++ tools/cm95_bathymetry_build_v4_tin_dense.py | 399 +++++++++++ tools/cm95_bathymetry_build_v5_adaptive.py | 588 ++++++++++++++++ tools/cm95_linux_full_rebuild_preflight.py | 110 +++ tools/cm95_v1_4_6_query_equivalence.py | 27 +- web/bathymetry_test.html | 191 +++++ web/bathymetry_test_server.py | 93 +++ 30 files changed, 6033 insertions(+), 60 deletions(-) create mode 100644 .env.linux.example create mode 100644 config/bathymetry/cm95_bathymetry_north_kyushu_v1.yaml create mode 100644 config/bathymetry/cm95_bathymetry_north_kyushu_v2_smooth.yaml create mode 100644 config/bathymetry/cm95_bathymetry_north_kyushu_v3_dense.yaml create mode 100644 config/bathymetry/cm95_bathymetry_north_kyushu_v4_tin_dense.yaml create mode 100644 config/bathymetry/cm95_bathymetry_north_kyushu_v5_adaptive.yaml create mode 100644 docs/CM95_Bathymetry独立PBF_v1.md create mode 100644 docs/CM95_S52_120个FAIL根因归并与修复批次规划_Codex任务_v1.md create mode 100644 docs/CM95_全球计划暂停与Japan_PBF_Canonical归一化缺口修复上下文_v2.md create mode 100644 docs/CM95_全球计划暂停与Japan_PBF低Zoom航标修复上下文_v1.md create mode 100755 scripts/setup_linux.sh create mode 100644 tools/cm95_bathymetry_build_national_v5.py create mode 100644 tools/cm95_bathymetry_build_v1.py create mode 100644 tools/cm95_bathymetry_build_v2_smooth.py create mode 100644 tools/cm95_bathymetry_build_v3_dense.py create mode 100644 tools/cm95_bathymetry_build_v4_tin_dense.py create mode 100644 tools/cm95_bathymetry_build_v5_adaptive.py create mode 100755 tools/cm95_linux_full_rebuild_preflight.py create mode 100644 web/bathymetry_test.html create mode 100644 web/bathymetry_test_server.py diff --git a/.env.linux.example b/.env.linux.example new file mode 100644 index 0000000..8cc27a1 --- /dev/null +++ b/.env.linux.example @@ -0,0 +1,9 @@ +CM93_DB_HOST=192.168.200.230 +CM93_DB_PORT=3306 +CM93_DB_NAME=cm95_japan +CM93_DB_USER=root +CM93_DB_PASSWORD=change_me +CM93_DB_CHARSET=utf8mb4 + +# Linux CM95 日本源数据目录。 +CM93_SOURCE_DIR=/root/sourceserver/cm95-data/CM93JAPAN diff --git a/.gitignore b/.gitignore index fb465e6..bdd1d69 100644 --- a/.gitignore +++ b/.gitignore @@ -43,4 +43,7 @@ cm93_spatial_audit_*/ .vscode/ .idea/ .DS_Store -Thumbs.db +Thumbs.db + +!.env.linux.example +.codex-build-logs/ diff --git a/README.md b/README.md index 49113c0..12a8925 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,8 @@ CM95 当前仓库的已验证基线，是一个将 CM93 数据解析、标准化并保存到 MySQL 5.7 的 importer。 +> Linux 迁移说明（2026-09-15）：历史 production baseline 仍为 MySQL 5.7；当前 Linux 环境已增加 MySQL 8.0.x 显式兼容层。MySQL 8 对 EPSG:4326 默认采用 SRS 轴顺序，因此所有 4326 WKB/WKT 构造与 WKB 读回必须通过项目 helper 强制 `axis-order=long-lat`，禁止在新代码中绕过 helper 直接拼空间 SQL。`192.168.200.230` 的 MySQL 8.0.45 已通过 schema 初始化、单 Cell 导入、64 点精度审计（0° 误差）、较大 Cell 精度审计及单 z12 PBF tile 构建；Japan 1492 Cell 全量 MySQL 8 production 验收尚未完成。 + 当前基线： ```text 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/config/bathymetry/cm95_bathymetry_north_kyushu_v2_smooth.yaml b/config/bathymetry/cm95_bathymetry_north_kyushu_v2_smooth.yaml new file mode 100644 index 0000000..a1e498c --- /dev/null +++ b/config/bathymetry/cm95_bathymetry_north_kyushu_v2_smooth.yaml @@ -0,0 +1,30 @@ +schema: navsea.cm95.bathymetry-build/v2 +build_id: cm95-bathymetry-north-kyushu-v2-smooth +dataset_id: 1 +bbox_wgs84: [129.7, 33.1, 130.6, 33.8] +source_layer: bathymetry_contour +source_kind: derived_sounding_idw_grid +reliability: reference_low +native_scale_max: 100000 +water_object_classes: [DEPARE, DRGARE] +contour_levels_m: [2, 5, 10, 20, 30, 50] +grid_spacing_m: 300 +idw_search_radius_m: 3000 +idw_max_nearest_m: 1500 +idw_min_neighbors: 3 +idw_max_neighbors: 16 +idw_power: 2.0 +chaikin_iterations: 2 +min_segment_length_m: 30 +regional_simplify_m: 3 +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-v2-smooth diff --git a/config/bathymetry/cm95_bathymetry_north_kyushu_v3_dense.yaml b/config/bathymetry/cm95_bathymetry_north_kyushu_v3_dense.yaml new file mode 100644 index 0000000..0cc5d0a --- /dev/null +++ b/config/bathymetry/cm95_bathymetry_north_kyushu_v3_dense.yaml @@ -0,0 +1,34 @@ +schema: navsea.cm95.bathymetry-build/v3 +build_id: cm95-bathymetry-north-kyushu-v3-dense +dataset_id: 1 +bbox_wgs84: [129.7, 33.1, 130.6, 33.8] +source_layer: bathymetry_contour +source_kind: derived_sounding_idw_grid_v3 +reliability: reference_low +native_scale_max: 100000 +water_object_classes: [DEPARE, DRGARE] +contour_levels_m: [2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 35, 40, 45, 50] +major_contours_m: [2, 5, 10, 20, 30, 50] +grid_spacing_m: 250 +idw_search_radius_m: 3000 +idw_max_nearest_m: 1800 +idw_min_neighbors: 3 +idw_max_neighbors: 16 +idw_power: 2.0 +gap_fill_max_component_nodes: 49 +gap_fill_radius_cells: 4 +gap_fill_min_neighbors: 4 +chaikin_iterations: 2 +min_segment_length_m: 25 +regional_simplify_m: 2 +minzoom: 8 +maxzoom: 12 +extent: 4096 +tile_buffer_px: 8 +simplify_m_by_zoom: + "8": 70 + "9": 35 + "10": 18 + "11": 8 + "12": 4 +output_dir: build/bathymetry/cm95-bathymetry-north-kyushu-v3-dense diff --git a/config/bathymetry/cm95_bathymetry_north_kyushu_v4_tin_dense.yaml b/config/bathymetry/cm95_bathymetry_north_kyushu_v4_tin_dense.yaml new file mode 100644 index 0000000..8676dbe --- /dev/null +++ b/config/bathymetry/cm95_bathymetry_north_kyushu_v4_tin_dense.yaml @@ -0,0 +1,29 @@ +schema: navsea.cm95.bathymetry-build/v4 +build_id: cm95-bathymetry-north-kyushu-v4-tin-dense +dataset_id: 1 +bbox_wgs84: [129.7, 33.1, 130.6, 33.8] +source_layer: bathymetry_contour +source_kind: derived_sounding_tin_surface_v4 +reliability: reference_low +native_scale_max: 100000 +water_object_classes: [DEPARE, DRGARE] +contour_levels_m: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50] +major_contours_m: [2, 5, 10, 20, 30, 50] +grid_spacing_m: 200 +max_triangle_edge_m: 4000 +surface_smooth_iterations: 1 +surface_smooth_strength: 0.35 +chaikin_iterations: 1 +min_segment_length_m: 20 +regional_simplify_m: 1.5 +minzoom: 8 +maxzoom: 12 +extent: 4096 +tile_buffer_px: 8 +simplify_m_by_zoom: + "8": 60 + "9": 30 + "10": 15 + "11": 7 + "12": 3 +output_dir: build/bathymetry/cm95-bathymetry-north-kyushu-v4-tin-dense diff --git a/config/bathymetry/cm95_bathymetry_north_kyushu_v5_adaptive.yaml b/config/bathymetry/cm95_bathymetry_north_kyushu_v5_adaptive.yaml new file mode 100644 index 0000000..3c0f0dc --- /dev/null +++ b/config/bathymetry/cm95_bathymetry_north_kyushu_v5_adaptive.yaml @@ -0,0 +1,37 @@ +schema: navsea.cm95.bathymetry-build/v5 +build_id: cm95-bathymetry-north-kyushu-v5-adaptive +dataset_id: 1 +bbox_wgs84: [129.7, 33.1, 130.6, 33.8] +source_layer: bathymetry_contour +source_kind: derived_sounding_adaptive_tin_v5 +reliability: reference_low +native_scale_max: 100000 +water_object_classes: [DEPARE, DRGARE] +contour_levels_m: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50] +major_contours_m: [2, 5, 10, 20, 30, 50] +grid_spacing_m: 150 +adaptive_edge_factor: 8.0 +adaptive_edge_min_m: 600 +adaptive_edge_max_m: 4000 +smooth_spacing_low_m: 100 +smooth_spacing_high_m: 700 +smooth_strength_low: 0.05 +smooth_strength_high: 0.30 +surface_smooth_iterations: 1 +gap_fill_max_component_nodes: 49 +gap_fill_radius_cells: 4 +gap_fill_min_neighbors: 4 +chaikin_iterations: 1 +min_segment_length_m: 18 +regional_simplify_m: 1.0 +minzoom: 8 +maxzoom: 12 +extent: 4096 +tile_buffer_px: 8 +simplify_m_by_zoom: + "8": 55 + "9": 28 + "10": 12 + "11": 5 + "12": 2 +output_dir: build/bathymetry/cm95-bathymetry-north-kyushu-v5-adaptive diff --git a/docs/CM95_Bathymetry独立PBF_v1.md b/docs/CM95_Bathymetry独立PBF_v1.md new file mode 100644 index 0000000..3e904b6 --- /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 192.168.200.184 --port 18890 +``` + +浏览： + +```text +http://192.168.200.184: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/docs/CM95_S52_120个FAIL根因归并与修复批次规划_Codex任务_v1.md b/docs/CM95_S52_120个FAIL根因归并与修复批次规划_Codex任务_v1.md new file mode 100644 index 0000000..f597b7b --- /dev/null +++ b/docs/CM95_S52_120个FAIL根因归并与修复批次规划_Codex任务_v1.md @@ -0,0 +1,417 @@ +# CM95 S-52 120 个 FAIL 根因归并与修复批次规划 Codex 任务 v1 + +日期：2026-09-20 + +## 1. 任务目标 + +本轮只做“120 个 FAIL 的根因归并与修复批次规划”，不修改任何产品功能代码。 + +唯一输入基线： + +- CM95 HEAD：`716acba74dec2c5bbc76e2cd62ca522618d810e2` +- 分支：`main` +- OpenCPN 固定基线：`f7211b4e1874c453e37a7f190aa191c55a25c32e` +- v2 审计报告：`docs/报告/CM95_S52_OpenCPN渲染流程对照审计报告_v2.md` +- v2 审计标准：`docs/CM95_S52_语义守恒与运行时能力审计标准_v2.md` +- 机器审计：`build/audits/cm95_s52_opencpn_code_first_equivalence_v2.json` + +当前机器审计结论固定为： + +- branch records：206 +- FAIL：120 +- NOT_PROVEN：86 +- 正式 PBF：5064 +- delivery signatures：1976 +- overall Target Portrayal：FAIL + +本轮不得把 86 个 NOT_PROVEN 当成已确认 Bug，也不得把 120 FAIL 直接当成 120 个独立代码缺陷。## 2. 核心问题 + +需要回答的不是“有哪些 120 个失败分支”，而是： + +1. 这 120 个 FAIL 最少可以归并为多少个独立根因？ +2. 每个根因实际发生在哪一层？ +3. 哪些 FAIL 只是同一根因的不同 branch 展开？ +4. 哪些看似相似的 FAIL 实际必须分开修？ +5. 哪些根因可以局部修复？ +6. 哪些根因说明当前架构缺少统一 Rule IR / executor？ +7. 哪些修复之间存在依赖顺序？ +8. 如何分批修复，尽量让一次修改消除一组 FAIL？ +9. 每个修复批次完成后应重跑哪些 Gate / branch / witness？ +10. 哪些 NOT_PROVEN 需要补证据，而不是写代码？ + +## 3. 本轮禁止事项 + +本轮禁止： + +- 修改 `src/`、`web/`、`tools/`、runtime registry、builder、production、style compiler 等产品实现； +- 修改 PBF schema； +- 修改 release gate； +- 修改 OpenCPN 固定参考源码； +- 新增或修改 Jest； +- 为了减少 FAIL 数量而改变 v2 审计判定； +- 把多个不同根因为了“数量更少”强行合并； +- 重新扫描 5064 个 PBF； +- 重新做完整 OpenCPN 源码审计； +- 从截图推断根因。 + +允许读取产品代码来定位根因，但不能编辑产品代码。## 4. 根因定义 + +一个 Root Cause 必须满足： + +- 有明确的错误机制或缺失机制； +- 能指向具体决策层、数据层或 renderer 层； +- 可以说明为什么会同时导致一组 branch FAIL； +- 修复该机制后，理论上可以消除其映射的全部或大部分 FAIL； +- 与其他 Root Cause 有清晰边界。 + +不要使用过宽的根因，例如： + +- “CM95 和 OpenCPN 不一致” +- “renderer 不完整” +- “CSP 有问题” + +应使用可执行粒度，例如： + +- “hazard resolver 没有把 QUAPNT01 suffix 拼入 ordered instruction chain” +- “presentation scalarization 只保留每种 primitive 的首个值，破坏 ordered multi-opcode sequence” +- “Web SOUNDG runtime 用普通 text 代替 SNDFRM02 glyph sequence executor” +- “RESTRN 集合分类把 14/15 错误纳入 info set” + +## 5. Root Cause ID + +所有根因使用稳定 ID： + +`RC-001`、`RC-002`、…… + +每个 RC 必须记录： + +- title +- layer +- mechanism +- exact evidence +- affected branch_ids +- affected FAIL count +- affected CSP / opcode +- current files/functions +- required change type +- architecture impact +- risk +- dependency +- proposed batch +- verification target +- confidence## 6. 归并原则 + +### 6.1 可以合并 + +当多个 FAIL： + +- 源于同一个函数/同一决策错误； +- 或源于同一个中间表示丢失； +- 或源于同一个 runtime consumer 缺失； +- 且一次修复可以共同消除， + +则应归并为一个 RC。 + +### 6.2 不应合并 + +即使表现相似，如果： + +- 修复位置不同； +- 条件语义不同； +- 一个是 Python resolver 错误，一个是 runtime executor 缺失； +- 一个是 lookup 错误，一个是 CSP 错误； +- 一个是 OpenCPN upstream anomaly / NOT_PROVEN； + +则不得为减少数量而合并。 + +### 6.3 多根因 branch + +一个 branch FAIL 可以同时受多个 RC 影响。 + +因此必须支持 many-to-many： + +`branch_id <-> root_cause_id` + +不能要求 120 个 FAIL 每个只能对应一个 RC。## 7. 层级分类 + +每个 RC 至少归入以下一类： + +- `NORMALIZATION_LOOKUP` +- `CSP_DECISION` +- `NESTED_CSP_CHAIN` +- `INSTRUCTION_IR` +- `PRESENTATION_MATERIALIZATION` +- `PBF_RUNTIME_FACTS` +- `STYLE_COMPILER` +- `RUNTIME_EXECUTOR` +- `RESOURCE_CONTRACT` +- `MARINER_STATE` +- `RENDER_ORDER_PRIORITY` +- `EVIDENCE_GAP` + +如果一个 RC 跨层，必须指定 primary layer 和 secondary layers。 + +## 8. 优先检查的已知根因簇 + +以下只作为 seed，不是结论，Codex 必须用 v2 JSON 验证后再决定是否拆分/合并： + +1. OBSTRN/WRECKS nested suffix 丢失： + - QUAPNT01 + - SOUNDG03 + - OBJNAM TX + - WRECKS Wk/WRECKS07 +2. ordered Rule Chain 被 scalar hint 压缩。 +3. Web/runtime 缺通用 multi-opcode executor。 +4. SOUNDG glyph chain 被普通文字替代。 +5. multi-SY / multi-LS / multi-TX 不守恒。 +6. RESTRN 集合语义错误。 +7. RESARE 集合语义 + OP side effect 不等价。 +8. SLCONS LOWACC / attribute presence 差异。 +9. LIGHTS05/06 颜色、半径、文本/描述/单位差异。 +10. DEPCNT SAFECD/SAFECN、quality、side effect 差异。 +11. DEPARE bad range / two_shades / DRGARE 差异。 +12. Lookup matching / fallback / $AREAS II25 等 Gate A 差异。## 9. 必须建立 FAIL → RC 覆盖矩阵 + +机器审计中的 120 个 FAIL branch 必须全部进入覆盖矩阵。 + +每条 FAIL 至少记录： + +| 字段 | 要求 | +|---|---| +| branch_id | v2 JSON 原始 branch ID | +| procedure | CSP/flow | +| fail_reason | v2 已确认原因 | +| root_cause_ids | 一个或多个 RC | +| primary_rc | 主要根因 | +| secondary_rc | 可选 | +| fix_batch | 修复批次 | +| verification | 修复后验证方式 | + +硬要求： + +- 120/120 FAIL 都必须有 root cause 映射； +- 不允许出现未解释的 FAIL； +- 如确实无法归因，单独建 `RC-XXX UNKNOWN_ROOT_CAUSE`，并明确为何无法继续，不得静默遗漏； +- 映射总数可以大于 120，因为允许一条 FAIL 对应多个 RC。 + +## 10. RC → FAIL 反向矩阵 + +每个 RC 必须列出： + +- affected_fail_count； +- affected_branch_ids； +- affected procedures； +- 是否存在 real-data signature； +- 是否只存在理论分支； +- 修复后预计可直接关闭的 FAIL； +- 修复后仍可能因其他 RC 保持 FAIL 的 branch。 + +必须区分： + +- `directly_resolved` +- `partially_unblocked` +- `still_blocked_by_other_rc`## 11. 独立根因数量 + +报告必须给出： + +- 根因总数； +- 按 layer 分类数量； +- 按 risk 分类数量； +- 按 proposed batch 分类数量。 + +同时给出两个指标： + +1. `raw_fail_count = 120` +2. `deduplicated_root_cause_count = N` + +不要人为追求 N 越小越好。 + +## 12. 修复批次规划 + +建立稳定 Batch ID： + +`BATCH-01`、`BATCH-02`、…… + +每个批次必须说明： + +- 目标； +- 包含哪些 RC； +- 为什么应一起修； +- 涉及的预计文件/函数； +- 是否改变公共数据契约； +- 是否改变 PBF； +- 是否改变 runtime； +- 是否需要重新 build release； +- 是否需要重新跑 5064 PBF audit； +- 是否需要视觉验收； +- 风险等级； +- 前置依赖； +- 完成判定； +- 回归范围。 + +## 13. 批次排序原则 + +默认排序： + +1. 明确且局部的 semantic bug； +2. 上游数据/IR 守恒问题； +3. 通用 Rule Chain / opcode executor； +4. renderer 专项； +5. LIGHTS 等复杂专项； +6. render order / priority 等高风险结构项； +7. NOT_PROVEN 补证据任务。 + +但最终排序必须依据实际 RC dependency graph，不得机械照抄。## 14. 架构分叉判断 + +每个 RC 必须判断属于： + +- `LOCAL_FIX` +- `SHARED_HELPER_FIX` +- `IR_CONTRACT_CHANGE` +- `GENERIC_EXECUTOR_REQUIRED` +- `SPECIALIZED_RENDERER_REQUIRED` +- `EVIDENCE_ONLY` + +特别回答： + +1. 是否继续逐 CSP 手工复制 OpenCPN 条件会制造新差异？ +2. 哪些 RC 可以靠修 `csp.py` 解决？ +3. 哪些 RC 即使修 `csp.py`，Web runtime 仍会 FAIL？ +4. 哪些 RC 证明需要 ordered Runtime Rule IR？ +5. 哪些 opcode 应由通用 executor 执行？ +6. SOUNDG、CA 等哪些应保留专用 renderer？ +7. 是否需要让 PBF 携带完整 ordered instruction IR，还是只需要 runtime facts + decision IR？ + +本轮只做设计判断，不实施。 + +## 15. NOT_PROVEN 单独处理 + +86 个 NOT_PROVEN 不是本轮根因归并主体。 + +但报告必须总结它们按“缺什么证据”归类，例如： + +- source attributes unavailable； +- association identity unavailable； +- upstream resource snapshot incomplete； +- global render scheduling unavailable； +- runtime state snapshot unavailable。 + +为这些项建立 `EVIDENCE-TASK-01...`，不要建立代码 RC，除非已有明确代码缺陷证据。 + +## 16. 本轮正式输出 + +本轮只允许新增/更新以下正式产物： + +1. `docs/报告/CM95_S52_120个FAIL根因归并与修复批次规划报告_v1.md` + - 完整列出 deduplicated Root Cause； + - 包含 120/120 FAIL → RC 覆盖矩阵； + - 包含 RC → FAIL 反向矩阵； + - 包含 RC dependency graph； + - 包含 BATCH-01... 修复批次规划； + - 包含 86 个 NOT_PROVEN 的 evidence-task 分类； + - 给出 raw_fail_count=120 与 deduplicated_root_cause_count=N； + - 明确哪些是局部修复、共享 helper、IR contract、通用 executor、专用 renderer、纯证据任务。 + +2. `build/audits/cm95_s52_120_fail_root_cause_plan_v1.json` + - 机器可读； + - 至少包含： + - cm95_head； + - opencpn_commit； + - source_audit_artifact； + - raw_fail_count； + - not_proven_count； + - deduplicated_root_cause_count； + - root_causes[]； + - fail_to_root_causes[]； + - batches[]； + - evidence_tasks[]； + - dependency_edges[]； + - coverage_summary； + - validation_summary。 + - 必须能机器确认 120/120 FAIL 全覆盖； + - 每个 FAIL 至少有一个 primary_rc； + - 多根因 branch 可列 secondary RC； + - 不得把 NOT_PROVEN 伪装成代码 RC。 + +不得创建第二份内容重复的“总结/摘要”文档。 + +## 17. 必做验证 + +完成前必须执行结构化自检，至少确认： + +- [ ] 输入审计 JSON 的 overall status 仍为 FAIL； +- [ ] 输入 branch records=206； +- [ ] 输入 FAIL=120； +- [ ] 输入 NOT_PROVEN=86； +- [ ] 120/120 FAIL 全部进入 fail_to_root_causes； +- [ ] 每个 FAIL 有且只有一个 primary_rc； +- [ ] secondary RC 可为 0..N； +- [ ] 所有引用的 RC ID 都存在； +- [ ] 所有 RC 的 affected_fail_count 与实际映射一致； +- [ ] deduplicated_root_cause_count 与 root_causes 数量一致； +- [ ] 所有 batch 引用的 RC 均存在； +- [ ] dependency graph 不包含未知 RC/Batch； +- [ ] NOT_PROVEN evidence tasks 与代码 RC 分离； +- [ ] 报告中的数量与 JSON 完全一致； +- [ ] 没有修改产品实现、Jest、PBF schema、release gate； +- [ ] 没有重新扫描 5064 个 PBF； +- [ ] 原工作区已有无关修改保持原样。 + +如果一个 FAIL 只能暂时归入 UNKNOWN_ROOT_CAUSE，必须在报告和 JSON 中明确列出，不得为了 100% 覆盖而虚构根因。 + +## 18. 修复批次输出要求 + +每个 BATCH 必须明确区分： + +- `expected_direct_fail_reduction`：若本批次单独完成，理论上可直接关闭的 FAIL 数； +- `expected_unblocked_fail_count`：本批次完成后可继续验证，但仍受其他 RC 阻塞的 FAIL 数； +- `requires_release_rebuild`； +- `requires_full_pbf_rescan`； +- `requires_runtime_validation`； +- `requires_visual_validation`； +- `risk`； +- `dependencies`。 + +不要把“预计可关闭”写成实际已关闭。本轮没有改代码，因此实际 FAIL 数仍是 120。 + +修复批次应尽量减少重复构建和重复全量验收；只有数据契约/PBF 变化时才安排全 release rebuild / 5064 PBF 重扫。纯 Python 语义修复、runtime executor 修复、renderer 专项应根据实际影响决定最小必要回归范围。 + +## 19. Git 收尾 + +当前仓库存在与本任务无关的未提交修改，必须严格隔离。 + +本轮只允许提交： + +- `docs/CM95_S52_120个FAIL根因归并与修复批次规划_Codex任务_v1.md` +- `docs/报告/CM95_S52_120个FAIL根因归并与修复批次规划报告_v1.md` +- `build/audits/cm95_s52_120_fail_root_cause_plan_v1.json` + +不得使用 `git add .`。 + +完成后： + +1. 显式 stage 上述三个路径； +2. 检查 `git diff --cached --name-only` 与 `--stat`； +3. 确认没有产品代码/Jest/原无关修改进入 staged； +4. commit 到当前 `main`； +5. commit message：`docs: plan fixes for OpenCPN S-52 audit failures`； +6. push `origin/main`； +7. 最终完成消息记录 commit hash、push 结果、deduplicated_root_cause_count、batch 数量及最高优先级批次。 + +push 失败时必须明确报告，不得声称完成。 + +## 20. 最终判定原则 + +本轮不是修复任务，不能因为规划完成就改变 v2 审计状态。 + +结束时必须保持： + +- raw FAIL = 120； +- raw NOT_PROVEN = 86； +- Target Portrayal = FAIL。 + +本轮成功标准是： + +`120 FAIL → 可解释的 Root Cause 集合 → 有依赖关系的修复 Batch → 明确回归目标` + +而不是“把失败数量写少”。 diff --git a/docs/CM95_全球计划暂停与Japan_PBF_Canonical归一化缺口修复上下文_v2.md b/docs/CM95_全球计划暂停与Japan_PBF_Canonical归一化缺口修复上下文_v2.md new file mode 100644 index 0000000..413ff9c --- /dev/null +++ b/docs/CM95_全球计划暂停与Japan_PBF_Canonical归一化缺口修复上下文_v2.md @@ -0,0 +1,558 @@ +# CM95 全球计划暂停与 Japan PBF Canonical 归一化缺口修复上下文 v2 + +> 本文取代 v1 中关于“低 Zoom retention 是主因”的判断。 +> +> 后续同 Zoom、同 Tile、同 Feature 审计已经证明：相当一部分 S-52 页面可见、NavSea canonical 页面不可见的对象，实际上仍存在于同一份 PBF 中，只是停留在 S-52 fallback source-layer，没有进入 NavSea canonical Layer Contract。 + +## 1. 当前决定 + +Global CM95 继续暂停： + +~~~ +Global CM95 Production = PAUSED +~~~ + +暂停决定继续有效，但原因修正为： + +- Japan 正式 PBF 中仍存在大量 S-52 fallback 对象； +- 完整 S-52 Style 可以正常显示这些对象； +- NavSea canonical Style 不消费这些 fallback layer； +- 因此 CM95 → NavSea canonical normalization 尚未完整； +- 在该缺口修复并验证前，不应把当前规则扩展到全球。 + +当前 Japan 优先级： + +~~~ +确认 S-52 fallback → canonical 缺口 +→ 完善 canonical mapping / composition +→ 同 Zoom S-52 vs canonical 对照 +→ 多区域验证 +→ NavSea 真机回归 +→ 再判断是否仍有独立 Zoom retention 问题 +~~~ + +## 2. 关键反证 + +用户提供的 Japan S-52 截图显示： + +- 较大 Zoom 可见大量航标； +- 较小 Zoom 仍然可见大量航标、浮标、顶标、灯标和其他 S-52 符号。 + +对应 X220 compare server： + +~~~ +web/japan_compare_server.py +~~~ + +实际使用： + +~~~ +Release: +cm95-d1-62a85bd40d3a-ae61b8866e80 + +PBF: +build/releases-s52-v2/cm95-d1-62a85bd40d3a-ae61b8866e80/pbf/delivery + +S-52 Style: +build/s52_style_compiler_japan_compare_v1/style.s52.json +~~~ + +因此不能再简单用“NavSea 低 Zoom 看不到”推导为“PBF 里没有”。 + +## 3. 同 Zoom / 同 Tile 审计 + +审计对象： + +~~~ +Tile: +11/1765/820 +~~~ + +同一个 z11 tile 中，带 S-52 symbol code 或 light symbol code 的 feature 共： + +~~~ +374 +~~~ + +其中： + +~~~ +navigation_marks: +77 +~~~ + +这 77 个对象： + +~~~ +chart_icon_image: 77/77 +canonical_object_type: 77/77 +~~~ + +说明已经成功 canonical 化的 navigation_marks 本身基本正常。 + +## 4. 真正的缺口：s52_points + +同一个 z11 tile 中： + +~~~ +s52_points 中带 S-52 symbol 的对象: +130 +~~~ + +而 NavSea canonical Style 对： + +~~~ +source-layer = s52_points +~~~ + +的 Layer 数量为： + +~~~ +0 +~~~ + +所以这些对象： + +- S-52 Style 可以显示； +- NavSea canonical Style 根本不会显示。 + +这是“同一个 PBF、同一个 Zoom，S-52 页面对象更多”的直接原因之一。 + +## 5. s52_points 中的重要对象 + +该 z11 tile 的 130 个 s52_points symbol 中： + +~~~ +fishing_facility 72 +s52_boylat 19 +s52_topmar 18 +s52_radrfl 14 +s52_marcul 3 +s52_boycar 2 +s52_cgusta 1 +s52_mstcon 1 +~~~ + +常见 S-52 symbol code： + +~~~ +FSHHAV01 +RADRFL03 +TOPMAR33 +TOPMAR22 +TOPMAR13 +TOPMAR05 +TOPMAR01 +BOYCON60 +BOYCON61 +BOYPIL60 +BOYPIL61 +BOYPIL68 +BOYSPR60 +BOYSPR61 +MARCUL02 +CGUSTA02 +MSTCON04 +~~~ + +特别是 BOYLAT、TOPMAR、BOYCAR、RADRFL，本身就是视觉上容易被理解为航标/浮标/顶标相关对象的类别。 + +## 6. Trace 已确认：不是被 Zoom 删除 + +正式 trace： + +~~~ +build/releases-s52-v2/ +cm95-d1-62a85bd40d3a-ae61b8866e80/ +reports/delivery/11/1765/820.trace.jsonl +~~~ + +BOYLAT 示例： + +~~~ +output_layer: +s52_points + +native_scale: +50000 + +mapping_rule_id: +S52-FALLBACK-BOYLAT + +mapping_status: +S52_FALLBACK +~~~ + +TOPMAR 示例： + +~~~ +output_layer: +s52_points + +native_scale: +50000 + +mapping_rule_id: +S52-FALLBACK-TOPMAR + +mapping_status: +S52_FALLBACK +~~~ + +RADRFL 示例： + +~~~ +output_layer: +s52_points + +native_scale: +50000 + +mapping_rule_id: +S52-FALLBACK-RADRFL + +mapping_status: +S52_FALLBACK +~~~ + +BOYCAR 示例： + +~~~ +output_layer: +s52_points + +mapping_rule_id: +S52-FALLBACK-BOYCAR + +mapping_status: +S52_FALLBACK +~~~ + +这些对象在 z11 中真实存在，而且 native_scale=50000。 + +所以至少这批对象与“z11 不读取 20k 数据”无关。 + +## 7. 正式 mapping 中没有对应 canonical 规则 + +当前： + +~~~ +config/mappings/cm95_to_navsea_v1.yaml +~~~ + +中没有以下正式 canonical mapping： + +~~~ +BOYLAT +TOPMAR +RADRFL +BOYCAR +MARCUL +CGUSTA +MSTCON +~~~ + +所以 Builder 在 mapping.resolve 找不到规则后进入 S-52 fallback 路径。 + +## 8. Builder 当前 fallback 行为 + +当前 Builder 中的 S-52 fallback 会把无 canonical mapping、但有 S-52 lookup 的对象写入： + +~~~ +POINT → s52_points +LINE → s52_lines +AREA → s52_areas +~~~ + +并设置： + +~~~ +canonical_object_type = s52_ +mapping_rule_id = S52-FALLBACK- +mapping_status = S52_FALLBACK +~~~ + +例如： + +~~~ +BOYLAT → s52_points / s52_boylat +TOPMAR → s52_points / s52_topmar +RADRFL → s52_points / s52_radrfl +~~~ + +这保证完整 S-52 Style 不丢对象，但 NavSea canonical Style 不消费这些 layer，因此看不到。 + +## 9. 另一个例子：FSHHAV + +正式 mapping 中 CM95-023： + +~~~ +CM95 class: +FSHHAV + +target_source_layer: +s52_points + +canonical_object_type: +fishing_facility + +mapping_status: +MAPPED + +style_fields: +空 +~~~ + +同一个 z11 tile 中这类对象有 72 个。 + +它们通过： + +~~~ +s52_symbol_code = FSHHAV01 +~~~ + +在 S-52 Style 中显示，但没有 chart_icon_image，因此 canonical Style 不显示。 + +所以问题不是 BOYLAT 的单一特例，而是更广泛的： + +~~~ +S-52 presentation layer +vs +NavSea canonical layer contract +~~~ + +边界不完整。 + +## 10. 对 v1 Zoom retention 判断的修正 + +v1 曾基于 z12→z11 stable FID / 空间匹配推断低 Zoom 有对象被裁掉。 + +现在确认，这种统计不能直接解释为“同一个现实对象被 Builder 删除”。 + +原因包括： + +1. 不同 native scale 可能使用不同 feature/FID 表达同一现实对象； +2. S-52 fallback 对象可能存在于另一个 source-layer； +3. canonical Style 没有消费 fallback layer； +4. 类型和位置可能因不同 cell / native scale 发生变化。 + +因此： + +~~~ +跨 Zoom stable FID retention +只能作为审计指标 +不能单独证明对象被删除 +~~~ + +当前不要据此修改 Japan z11 scale ladder。 + +## 11. 当前最可能主根因 + +当前证据最强的主链路： + +~~~ +CM95 PBF 中存在正式 S-52 可显示对象 + ↓ +部分对象只有 S-52 presentation semantics + ↓ +Builder 输出到 s52_points / s52_lines / s52_areas + ↓ +NavSea canonical Style 不消费这些 Layer + ↓ +NavSea 看起来“丢对象” +~~~ + +所以当前问题首先归类为： + +~~~ +CM95 → NavSea canonical normalization completeness +~~~ + +而不是： + +~~~ +PBF Zoom generalization +~~~ + +## 12. CM95 项目下一步 + +### 12.1 暂时不要先改 Zoom policy + +当前不要因为 v1 的判断直接修改： + +~~~ +z11 scale ladder +20k retention +navigation aid retention ladder +~~~ + +必须先完成 canonical completeness audit。 + +### 12.2 建立同 Zoom completeness audit + +对同一个 PBF tile： + +~~~ +完整 S-52 Style 可显示 feature +vs +NavSea canonical Style 可显示 feature +~~~ + +逐 feature 分类： + +~~~ +A. 已 canonical，能显示 +B. 已 canonical，但 Style 漏匹配 +C. S52_FALLBACK，无 canonical mapping +D. MAPPED，但仍只保留在 s52_points/s52_lines/s52_areas +E. 需要 composite/relation 后才应进入 canonical +F. 明确不应该进入 NavSea canonical +~~~ + +### 12.3 优先审计类别 + +至少： + +~~~ +BOYLAT +TOPMAR +BOYCAR +RADRFL +MARCUL +CGUSTA +MSTCON +FSHHAV +~~~ + +特别是： + +~~~ +BOYLAT +TOPMAR +BOYCAR +~~~ + +应检查是否应成为 canonical navigation_marks 的组成部分或独立标准对象。 + +不能简单把全部 s52_points 直接显示到 NavSea。 + +### 12.4 必须保持 composite 语义 + +当前 navigation_marks 中已经存在大量： + +~~~ +mapping_rule_id = CM95-NAV-COMPOSITE +mapping_status = COMPOSITE +~~~ + +所以 BOYLAT / TOPMAR 等是否独立输出，必须考虑： + +- 是否已经被 composite 消费； +- relation20； +- 浮标 + 顶标组合； +- 灯 + 标体组合； +- 是否会重复显示； +- stable FID； +- S-52 presentation order。 + +## 13. 历史线索 + +仓库历史/实验 trace 中存在某些： + +~~~ +S52-FALLBACK-BOYLAT +S52-FALLBACK-TOPMAR +~~~ + +曾输出到： + +~~~ +navigation_marks +~~~ + +而当前正式 Japan release 同类对象输出到： + +~~~ +s52_points +~~~ + +这说明项目历史上可能已经尝试过相关整合逻辑。 + +CM95 项目可以调查这些历史实现及其放弃/调整原因，但不要直接恢复旧代码。 + +## 14. Global 暂停原因更新 + +Global CM95 继续暂停。 + +正式原因更新为： + +~~~ +Japan 正式 Release 已证明： +S-52 presentation coverage +大于 +NavSea canonical coverage +~~~ + +也就是 canonical normalization 尚未完整。 + +恢复 Global 之前至少需要： + +1. Japan canonical completeness audit 完成； +2. 安全相关 S52_FALLBACK 分类完成； +3. canonical mapping/composite 规则冻结； +4. 同 Zoom S-52 vs canonical 重要对象覆盖通过； +5. 不产生 duplicate navigation mark； +6. stable FID 不破坏； +7. Pickup / ChartDanger 所需字段完整； +8. Japan 多区域抽样通过； +9. NavSea 真机回归通过； +10. 再单独判断是否仍有真正 Zoom retention 问题。 + +## 15. NavSea 边界保持不变 + +NavSea 项目不要为了 CM95 修改： + +~~~ +1/2/3/4 挡 +Pickup +ChartDangerAlarmSource +CrossCursor +PbfLayerGroupConfig +MapLibre query logic +~~~ + +正确方向： + +~~~ +CM95 项目完善 canonical PBF contract +→ NavSea 原逻辑直接消费 +~~~ + +## 16. 当前结论摘要 + +~~~ +Global CM95: +PAUSED + +Japan PBF: +对象并非简单低 Zoom 消失 + +正式 z11 PBF: +存在大量 S-52 可显示对象 + +navigation_marks canonical: +已映射对象基本正常 + +主要缺口: +s52_points / s52_lines / s52_areas 中的 +S52_FALLBACK / presentation-only 对象 +没有进入 NavSea canonical contract + +最优先任务: +同 Zoom canonical completeness audit + +当前不要先做: +z11 scale ladder 修改 +20k 全局下放 +navigation aid retention policy +~~~ diff --git a/docs/CM95_全球计划暂停与Japan_PBF低Zoom航标修复上下文_v1.md b/docs/CM95_全球计划暂停与Japan_PBF低Zoom航标修复上下文_v1.md new file mode 100644 index 0000000..8c5b86c --- /dev/null +++ b/docs/CM95_全球计划暂停与Japan_PBF低Zoom航标修复上下文_v1.md @@ -0,0 +1,655 @@ +# CM95 全球计划暂停与 Japan PBF 低 Zoom 航标修复上下文 v1 + +## 1. 文档目的 + +本文用于把 NavSea 侧已经确认的事实、CM95 Japan PBF 暴露的问题，以及后续应由 CM95 项目承担的工作边界正式交接到 CM95 项目。 + +本文是上下文与状态交接文档，不是实现报告，也不代表 Japan PBF 修复已经完成。 + +后续 CM95/PBF 工作以 X220 上的实际仓库为唯一开发和生产事实源： + +~~~ +机器：X220 +仓库：/root/sourceserver/cm95 +~~~ + +CARBON7i 上曾做过的 PBF 实验只能作为调查线索，不作为正式实现、生产基线或最终结论。 + +--- + +## 2. 当前正式决定 + +### 2.1 全球 CM95 计划暂停 + +从现在开始： + +~~~ +Global CM95 Production = PAUSED +~~~ + +暂停原因： + +Japan 区域已经确认存在与 Zoom 有关的助航对象保留缺陷。若在修复前继续扩展全球 PBF，很可能把同一生成缺陷复制到全球数据，之后产生大规模返工。 + +因此，在 Japan 修复、验证和冻结完成之前： + +- 不继续扩大全球 CM95 生成范围； +- 不把当前 Japan 的 PBF scale/generalization 策略直接推广到全球； +- 不将当前 Japan 的低 Zoom 行为视为全球可接受基线。 + +### 2.2 Japan PBF 优先修复 + +当前 CM95 项目优先级改为： + +~~~ +Japan PBF correctness +→ Japan multi-region validation +→ Japan production freeze +→ 再评估 Global CM95 resume +~~~ + +--- + +## 3. NavSea 侧当前架构结论 + +NavSea 端已经重新收敛到以下架构： + +~~~ +CM95 PBF + ↓ +NavSea canonical source-layer / property contract + ↓ +NavSea canonical Style + ↓ +现有 1/2/3/4 挡 +现有 Pickup +现有 ChartDangerAlarmSource +~~~ + +禁止继续走下面的方向： + +~~~ +CM95 raw S-52 Style + ↓ +让 NavSea 为 CM95 特别修改挡位 / Pickup / 报警 +~~~ + +结论： + +CM95 Provider 应负责把自己的数据生产为 NavSea 已定义好的 canonical Layer / property contract；NavSea 业务模块不应该理解 CM95 特有 Layer。 + +--- + +## 4. NavSea 真机已经验证通过的部分 + +测试机： + +~~~ +192.168.200.145 +~~~ + +### 4.1 1/2/3/4 挡 + +已经实际验证： + +~~~ +3 → 2 → 1 → 4 +~~~ + +均正常切换。 + +同时观察： + +~~~ +MapLibre resource cancel count = 0 +ParseStyle error = 0 +resourceUrl error = 0 +~~~ + +说明当前挡位已经回到原 NavSea existing-layer visibility 热切换，不再通过重建整张 Style/MapView 实现。 + +### 4.2 Pickup + +连续点击多个 CM95 航标，现有链路可以直接返回 NavSea 标准 Pickup 信息。 + +现有正式链路： + +~~~ +runNavSeaPickupJson +→ pickupInterpreter +→ CrossCursor +~~~ + +实际返回过： + +- 小灯； +- 航路標識灯； +- red / green / white； +- omni； +- 距离与方位。 + +因此 Pickup 不需要增加 CM95 专用逻辑。 + +### 4.3 海图报警 + +现有 ChartDangerAlarmSource 仍使用 NavSea canonical render layer 和 canonical properties。 + +已经确认 CM95 canonical PBF 中存在相关： + +- navigation_hazard_point； +- navigation_marks； +- baseline_area； +- facility / land / coast / structure 等标准 source-layer。 + +动态浅水报警 danger-depth-overlay 继续通过当前 Provider 的 analysisSourceId + baseline_area 生成。 + +CM95 的 baseline_area 已观察到标准深度分区，例如： + +~~~ +depth_zone_0-2m +depth_zone_0-5m +depth_zone_0-10m +depth_zone_2-5m +depth_zone_5-10m +depth_zone_10-20m +... +~~~ + +因此 ChartDangerAlarmSource 不应为了 CM95 改业务逻辑。 + +--- + +## 5. Canonical Style 已经排除为“丢灯”根因 + +在福冈湾真实 CM95 PBF 上做过逐对象检查。 + +对 PBF 中已经存在的灯类 navigation_marks feature，逐个执行 NavSea canonical Style 的灯标 filters。 + +结果： + +~~~ +unmatched = 0 +~~~ + +也就是说： + +只要 CM95 PBF 中已经存在该灯标，NavSea canonical Style 可以匹配并显示。 + +所以剩余的“低 Zoom 丢灯”问题发生在 Style 之前： + +~~~ +CM95 source +→ PBF Builder / scale selection / generalization / retention +→ PBF +→ NavSea canonical Style +~~~ + +当前不应继续通过修改 NavSea Style filter 来补偿。 + +--- + +## 6. 已确认的低 Zoom 现象 + +在固定福冈湾区域做过跨 Zoom 审计。 + +CM95 灯类对象在 Zoom 降低时存在明显减少。 + +特别是 z12 → z11： + +- z12 中存在的一部分灯标在 z11 消失； +- 主要缺失对象中包含 minor_light； +- 使用空间匹配方式时，CM95 z12→z11 的保留率显著低于 Current NavSea。 + +曾观察到大致： + +~~~ +CM95 z12 → z11 +50m 内保留率约 85% +100m 内保留率约 87% +~~~ + +而同类 Current NavSea 对照接近： + +~~~ +z12 → z11 ≈ 100% +~~~ + +因此 z11 已经可以明确视为 Japan PBF 的一个重点修复 Zoom。 + +z10 以下也存在更激进的数据缩减，但不应简单假设所有 z12 对象都必须无条件保留到全部低 Zoom；需要由 CM95 项目建立明确的安全对象 retention policy。 + +--- + +## 7. Current 与 CM95 原始数据并非完全一致 + +Current 与 CM95 是不同数据源。 + +已观察到类似差异： + +~~~ +Current: breakwater_lighthouse +CM95 附近对象: light_beacon +~~~ + +以及： + +~~~ +Current: harbor_lighthouse +CM95: coastal_lighthouse_over_15m +~~~ + +还存在几十米乃至更大的坐标差异。 + +因此禁止用以下方式修复： + +~~~ +距离很近 +→ 认为是同一个对象 +→ 自动把 CM95 类型改成 Current 类型 +~~~ + +航标 canonical 分类必须依据 CM93/CM95 自身语义事实，例如： + +- object class； +- relation/group； +- native scale； +- S-52 lookup； +- 原始 attributes； +- 已冻结的 mapping contract。 + +空间邻近只能用于审计和发现异常，不能作为生产语义改判依据。 + +--- + +## 8. CARBON7i 曾做过的调查线索 + +以下内容是在确认“X220 才是实际 CM95/PBF 工作机器”之前，在 CARBON7i 副本上做的实验。 + +这些内容只能作为 X220 调查线索，必须在 X220 当前真实代码、配置和数据上重新验证。不得直接复制实验代码。 + +### 8.1 一个强候选根因 + +CARBON7i 副本中的 Japan policy 曾显示类似： + +~~~ +z11: +50k +100k +200k +... + +z12: +7.5k +20k +50k +100k +... +~~~ + +同时福冈区域数据库中大量 LIGHTS 存在于： + +~~~ +20k +50k +100k +... +~~~ + +由此得到一个需要在 X220 正式重新验证的强候选： + +~~~ +z12 会读取 20k LIGHTS +z11 从 50k 开始 +→ 20k 独有灯标在 z11 被整体排除 +~~~ + +### 8.2 不应把整个 z11 提升为 20k full detail + +CARBON7i 的实验还显示，如果简单把整个 z11 scale ladder 改为包含 20k： + +代表 tile 曾出现： + +~~~ +旧 z11: +约 3906 features + +整个 z11 加入 20k: +约 6851 features +~~~ + +约增加 75%。 + +原因是 20k 中不仅有灯标，还包含大量： + +- seabed text； +- area； +- line； +- label； +- 港区细节； +- 其他 feature。 + +因此不建议把整个 z11 直接升级为 20k/full detail。 + +### 8.3 值得在 X220 重新验证的设计方向 + +CARBON7i 实验曾验证一个更合理的方向： + +为安全助航对象建立独立的 scale retention ladder，例如： + +~~~ +普通 z11 海图内容: +50k → 100k → 200k → ... + +navigation_marks z11: +7.5k → 20k → 50k → 100k → ... +~~~ + +只针对助航组件，例如： + +~~~ +LIGHTS +BCNLAT +BCNSPP +BOYSPP +RTPBCN +TOWERS +_extgn +~~~ + +使用细尺度 quilt。 + +其他 source-layer 仍按正常 z11 ladder 生成。 + +CARBON7i 实验曾得到过以下结果： + +单代表 tile： + +~~~ +z12 → z11 navigation_marks stable FID retention: +约 16.8% → 100% + +总 feature: +3906 → 3933 + +PBF 体积: +约 +1% + +非 navigation_marks layer count: +0 差异 +~~~ + +福冈 12 张 z11 tile 实验曾得到： + +~~~ +z12 navigation_marks FID retention: +约 35.8% → 100% + +非 navigation layer count 差异: +0 + +总 PBF 体积: +约 +1.1% + +hard gate: +0 + +encoder drop: +0 +~~~ + +这些只是实验线索。 + +正式 CM95 项目必须在 X220 当前实际 builder、policy、mapping、manifest、release 上重新实现和验证。 + +--- + +## 9. X220 当前事实源 + +后续开发机器： + +~~~ +X220 +~~~ + +正式仓库： + +~~~ +/root/sourceserver/cm95 +~~~ + +本文建立时观察到的分支： + +~~~ +main +~~~ + +本文建立前调查时观察到的 HEAD： + +~~~ +ed53864598f3059f524841fbe25dcd488fc1301a +~~~ + +注意：X220 当前 worktree 自身已经存在未提交修改。 + +因此后续任务必须： + +1. 先读当前实际 diff； +2. 不覆盖现有未提交工作； +3. 只做本任务相关最小修改； +4. 不把无关文件一起提交； +5. 以 X220 当前实际文件为准，不使用 CARBON7i 副本替代。 + +X220 中已经确认存在： + +~~~ +/root/sourceserver/cm95/config/pbf/ +/root/sourceserver/cm95/src/cm95pbf/ +/root/sourceserver/cm95/build/production/ +/root/sourceserver/cm95/build/releases-s52-v2/ +/root/sourceserver/cm95/scripts/ +/root/sourceserver/cm95/tools/ +~~~ + +并存在 Japan production plan： + +~~~ +build/production/cm95_s52_japan_production_plan_v1.json +~~~ + +--- + +## 10. CM95 项目后续正式工作 + +### 10.1 第一阶段：X220 重新确认根因 + +必须基于 X220 当前实际版本重新检查： + +~~~ +config/pbf/cm95_tile_policy_japan_full_detail_v1.yaml +src/cm95pbf/builder.py +src/cm95pbf/production.py +build/production/cm95_s52_japan_production_plan_v1.json +~~~ + +至少回答： + +- z11 实际 scale ladder； +- z12 实际 scale ladder； +- LIGHTS / beacon / buoy 等助航对象的 native scale 分布； +- z12 中存在但 z11 消失的对象主要来自哪些 native scale； +- 20k LIGHTS 是否确实是主要来源； +- 当前 builder 是如何做 scale-region quilt 和 finer-scale coverage suppression 的。 + +### 10.2 不修改 importer + +当前问题属于 PBF production 阶段。 + +除非出现新的直接证据，否则不要修改： + +~~~ +CM93 parser +geometry normalization +MySQL importer +raw source records +~~~ + +### 10.3 优先研究 navigation aid 独立 retention + +优先方向： + +~~~ +navigation_marks 使用独立 scale ladder +~~~ + +而不是： + +~~~ +整个 z11 直接启用 20k/full detail +~~~ + +实现时应保持： + +- stable FID 不变； +- composite relation 语义不变； +- 多 LIGHTS 合成逻辑不变； +- S-52 presentation 字段不退化； +- trace / report / hard gate 继续准确； +- 其他 source-layer 不被细尺度数据无意放大。 + +### 10.4 单 Tile A/B 验证 + +至少比较： + +~~~ +旧 z11 +候选修复 z11 +对应 z12 四个 child tiles +~~~ + +指标： + +- navigation_marks stable FID retention； +- canonical_object_type 分布； +- 非 navigation layer feature count； +- PBF byte size； +- duplicate FID； +- hard gate； +- encoder drop； +- build report； +- trace consistency。 + +### 10.5 Japan 多区域抽样 + +正式冻结前至少抽样： + +- 函馆； +- 仙台； +- 东京湾； +- 大阪湾； +- 广岛湾； +- 福冈； +- 唐津； +- 长崎； +- 那霸。 + +不能只凭福冈一片区域通过就宣布 Japan 修复完成。 + +### 10.6 Japan 全量重建 + +不要边全国构建边修改 policy。 + +应先冻结： + +~~~ +builder contract +Japan tile policy +mapping +audit tool +acceptance thresholds +~~~ + +然后再生成正式 Japan production plan 并全量构建。 + +--- + +## 11. Global CM95 恢复条件 + +Global CM95 不因为 Japan 局部修复通过就立即恢复。 + +建议至少满足： + +1. Japan z11 navigation aid retention 根因确认； +2. Japan 修复实现完成； +3. stable FID 不破坏； +4. composite navigation aids 无重复； +5. 非 navigation layer 无非预期变化； +6. hard gate PASS； +7. encoder drop = 0； +8. PBF 体积增长在可接受范围； +9. 多区域抽样通过； +10. Japan 全量 production audit 通过； +11. Japan release 完成； +12. NavSea 真机低 Zoom 显示回归通过。 + +达到以上条件后，才允许讨论： + +~~~ +Global CM95 Production = RESUME +~~~ + +即使恢复全球计划，也应先把 Japan 修复抽象成 provider-neutral / production-safe contract，再逐区域扩大，不应直接把 Japan 特例硬编码到全球路径。 + +--- + +## 12. NavSea 项目边界 + +CM95 修复期间，NavSea 项目不要再修改以下模块来补偿 CM95 数据： + +~~~ +1/2/3/4 挡 +Pickup +ChartDangerAlarmSource +CrossCursor +PbfLayerGroupConfig +MapLibre query compatibility +~~~ + +CM95 最终交付的 PBF 应满足现有 NavSea canonical Layer / property contract。 + +NavSea 只消费，不承担 CM95 PBF 数据补全和 zoom retention。 + +--- + +## 13. 当前状态摘要 + +~~~ +NavSea CM95 canonical 接入： + 基本通过 + +NavSea 1/2/3/4 挡： + 真机通过 + +NavSea Pickup： + 真机通过 + +NavSea ChartDanger contract： + 已核对 + +Canonical Style filter 漏灯： + 已排除 + +CM95 Japan z11 低 Zoom 航标保留： + 已确认存在问题 + +Global CM95： + PAUSED + +Japan PBF： + 下一优先修复目标 + +实际工作机器： + X220 + +正式仓库： + /root/sourceserver/cm95 +~~~ diff --git a/scripts/setup_linux.sh b/scripts/setup_linux.sh new file mode 100755 index 0000000..5961824 --- /dev/null +++ b/scripts/setup_linux.sh @@ -0,0 +1,51 @@ +#!/usr/bin/env bash +set -euo pipefail + +PROJECT_ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +SOURCE_ROOT="${CM93_SOURCE_DIR:-/root/sourceserver/cm95-data/CM93JAPAN}" +VENV="$PROJECT_ROOT/.venv" + +cd "$PROJECT_ROOT" + +if [[ ! -d "$SOURCE_ROOT" ]]; then + echo "CM93 source not found: $SOURCE_ROOT" >&2 + exit 2 +fi + +python3 - <<'PY' +import sys +if sys.version_info[:2] != (3, 12): + raise SystemExit(f"Python 3.12 required, got {sys.version.split()[0]}") +print(f"Python {sys.version.split()[0]}: PASS") +PY + +if [[ ! -x "$VENV/bin/python" ]]; then + python3 -m venv "$VENV" +fi + +"$VENV/bin/python" -m pip install --upgrade pip setuptools wheel +"$VENV/bin/pip" install -e . + +if [[ ! -f .env ]]; then + cp .env.linux.example .env + chmod 600 .env + echo "Created .env from .env.linux.example" +else + chmod 600 .env + echo "Existing .env kept" +fi + +"$VENV/bin/python" tools/cm95_linux_full_rebuild_preflight.py +"$VENV/bin/cm93db" scan "$SOURCE_ROOT" +"$VENV/bin/cm93db" inspect-dicts "$SOURCE_ROOT" + +if grep -Eq '^CM93_DB_PASSWORD=(|change_me)$' .env; then + echo + echo "Linux source/code preflight: PASS" + echo "Database credential: PENDING (CM93_DB_PASSWORD is empty)" + echo "Edit $PROJECT_ROOT/.env before DB/full-rebuild operations." +else + echo + echo "Linux source/code preflight: PASS" + echo "Database password is configured; run cm93db env-check before rebuild." +fi diff --git a/src/cm93db/audit.py b/src/cm93db/audit.py index a884ac8..ae7ee11 100644 --- a/src/cm93db/audit.py +++ b/src/cm93db/audit.py @@ -10,7 +10,13 @@ from shapely import from_wkb from .cm93_binary import parse_cell from .config import DbConfig -from .database import connect, assert_mysql57, server_version +from .database import ( + connect, + assert_mysql_compatible, + server_version, + mysql_compat_mode_for_connection, + as_wkb_sql, +) from .dictionaries import load_object_dictionary, load_attribute_dictionary from .geometry import build_coverages, build_feature from .scanner import scan_cm93, CellFile @@ -173,13 +179,14 @@ def _dataset_id(conn, root: Path) -> int | None: def _precision_db_geom(conn, cell_id: int, sample: PrecisionSample): + db_mode = mysql_compat_mode_for_connection(conn) if sample.geometry_type == "SOUNDING": - sql=("SELECT ST_AsWKB(s.geom) AS wkb FROM chart_feature f JOIN chart_sounding_point s ON s.feature_id=f.id " + sql=(f"SELECT {as_wkb_sql(db_mode, 's.geom')} AS wkb FROM chart_feature f JOIN chart_sounding_point s ON s.feature_id=f.id " "WHERE f.cell_id=%s AND f.source_record_index=%s AND s.point_index=%s") params=(cell_id,sample.record_index,sample.point_index or 0) else: table={"POINT":"chart_point_geometry","LINE":"chart_line_geometry","AREA":"chart_area_geometry"}[sample.geometry_type] - sql=(f"SELECT ST_AsWKB(g.geom) AS wkb FROM chart_feature f JOIN {table} g ON g.feature_id=f.id " + sql=(f"SELECT {as_wkb_sql(db_mode, 'g.geom')} AS wkb FROM chart_feature f JOIN {table} g ON g.feature_id=f.id " "WHERE f.cell_id=%s AND f.source_record_index=%s") params=(cell_id,sample.record_index) with conn.cursor() as cur: @@ -207,7 +214,7 @@ def run_audit(cfg: DbConfig, root: Path, limit: int | None = None, scales: set[s conn=connect(cfg) try: - assert_mysql57(server_version(conn)) + assert_mysql_compatible(server_version(conn)) did=_dataset_id(conn,root) if did is None: result.db_missing_cells=len(stats) diff --git a/src/cm93db/cli.py b/src/cm93db/cli.py index a864f03..ea0c4bc 100644 --- a/src/cm93db/cli.py +++ b/src/cm93db/cli.py @@ -7,7 +7,7 @@ from datetime import datetime from pathlib import Path from .config import load_config, load_source_dir -from .database import connect, server_version, assert_mysql57, create_database, apply_schema +from .database import connect, server_version, assert_mysql_compatible, mysql_compat_mode, create_database, apply_schema from .scanner import scan_cm93 from .dictionaries import load_object_dictionary, load_attribute_dictionary from .schema_contract import REQUIRED_TABLES @@ -46,7 +46,7 @@ def cmd_env_check(_args): try: v = server_version(conn) print(f"MySQL: {v}") - assert_mysql57(v) + assert_mysql_compatible(v) with conn.cursor() as cur: cur.execute("SHOW VARIABLES LIKE 'max_allowed_packet'") row=cur.fetchone() @@ -54,7 +54,7 @@ def cmd_env_check(_args): val=int(row.get("Value",0)); print(f"max_allowed_packet: {val} bytes") if val < 4*1024*1024: print("WARNING: max_allowed_packet 较小；建议至少 4 MiB。raw table 已按 512 KiB 分块。") - print("MySQL 5.7 compatibility: OK") + print(f"MySQL compatibility: {mysql_compat_mode(v)} OK") finally: conn.close() @@ -74,9 +74,9 @@ def cmd_db_check(_args): conn = connect(cfg) try: v = server_version(conn) - assert_mysql57(v) + assert_mysql_compatible(v) with conn.cursor() as cur: - cur.execute("SELECT table_name FROM information_schema.tables WHERE table_schema=%s", (cfg.database,)) + cur.execute("SELECT TABLE_NAME AS table_name FROM information_schema.tables WHERE table_schema=%s", (cfg.database,)) tables = {row["table_name"] for row in cur.fetchall()} missing = sorted(REQUIRED_TABLES - tables) print(f"MySQL: {v}") diff --git a/src/cm93db/database.py b/src/cm93db/database.py index bf3aa71..f4de612 100644 --- a/src/cm93db/database.py +++ b/src/cm93db/database.py @@ -27,23 +27,73 @@ def server_version(conn) -> str: return str(cur.fetchone()["v"]) -def assert_mysql57(version: str) -> None: +def mysql_version_tuple(version: str) -> tuple[int, int, int]: m = re.match(r"^(\d+)\.(\d+)\.(\d+)", version) if not m: raise RuntimeError(f"无法识别 MySQL 版本: {version}") - major, minor, _patch = map(int, m.groups()) - if (major, minor) != (5, 7): - raise RuntimeError( - f"项目数据库基线是 MySQL 5.7.34；当前服务器为 {version}。" - "不要直接继续初始化，先确认兼容策略。" - ) + return tuple(map(int, m.groups())) + + +def mysql_compat_mode(version: str) -> str: + major, minor, _patch = mysql_version_tuple(version) + if (major, minor) == (5, 7): + return "mysql57" + if (major, minor) == (8, 0): + return "mysql80" + raise RuntimeError( + f"CM95 当前仅验证 MySQL 5.7.x / 8.0.x；当前服务器为 {version}。" + ) + + +def assert_mysql_compatible(version: str) -> None: + mysql_compat_mode(version) + + +_MYSQL_COMPAT_MODE_CACHE: dict[tuple[str, int, str], str] = {} + + +def mysql_compat_mode_for_connection(conn) -> str: + key = ( + str(getattr(conn, "host", "")), + int(getattr(conn, "port", 0) or 0), + str(getattr(conn, "user", "")), + ) + cached = _MYSQL_COMPAT_MODE_CACHE.get(key) + if cached is not None: + return cached + mode = mysql_compat_mode(server_version(conn)) + _MYSQL_COMPAT_MODE_CACHE[key] = mode + return mode + + +def assert_mysql57(version: str) -> None: + # Backward-compatible public name retained for older tools/scripts. + assert_mysql_compatible(version) + + +def geom_from_wkb_sql(mode: str, value_sql: str = "%s", srid: int = 4326) -> str: + if mode == "mysql80": + return f"ST_GeomFromWKB({value_sql},{srid},'axis-order=long-lat')" + return f"ST_GeomFromWKB({value_sql},{srid})" + + +def geom_from_text_sql(mode: str, value_sql: str = "%s", srid: int = 4326) -> str: + if mode == "mysql80": + return f"ST_GeomFromText({value_sql},{srid},'axis-order=long-lat')" + return f"ST_GeomFromText({value_sql},{srid})" + + +def as_wkb_sql(mode: str, geom_sql: str) -> str: + if mode == "mysql80": + return f"ST_AsWKB({geom_sql},'axis-order=long-lat')" + return f"ST_AsWKB({geom_sql})" def create_database(cfg: DbConfig) -> None: safe = cfg.database.replace("`", "``") conn = connect(cfg, with_database=False) try: - assert_mysql57(server_version(conn)) + assert_mysql_compatible(server_version(conn)) with conn.cursor() as cur: cur.execute( f"CREATE DATABASE IF NOT EXISTS `{safe}` " @@ -90,7 +140,7 @@ def apply_schema(cfg: DbConfig, schema_path: Path) -> None: sql = schema_path.read_text(encoding="utf-8") conn = connect(cfg) try: - assert_mysql57(server_version(conn)) + assert_mysql_compatible(server_version(conn)) with conn.cursor() as cur: for stmt in _split_sql(sql): cur.execute(stmt) conn.commit() diff --git a/src/cm93db/importer.py b/src/cm93db/importer.py index c12fe2f..96eab70 100644 --- a/src/cm93db/importer.py +++ b/src/cm93db/importer.py @@ -12,7 +12,13 @@ from shapely.geometry.base import BaseGeometry from . import __version__ from .cm93_binary import ParsedCell, ParsedAttribute, parse_cell, CM93ParseError from .config import DbConfig -from .database import connect, assert_mysql57, server_version +from .database import ( + connect, + assert_mysql_compatible, + server_version, + mysql_compat_mode_for_connection, + geom_from_wkb_sql, +) from .dictionaries import load_object_dictionary, load_attribute_dictionary, ObjectDef, AttributeDef from .geometry import build_coverages, build_feature, BuiltFeature from .scanner import scan_cm93, CellFile @@ -489,6 +495,8 @@ def _insert_parsed_cell( ) -> tuple[int, int, int, int, int]: _clear_cell_payload(conn, cell_id) h = cell.header + db_mode = mysql_compat_mode_for_connection(conn) + geom_from_wkb = geom_from_wkb_sql(db_mode) with conn.cursor() as cur: cur.execute( @@ -593,7 +601,7 @@ def _insert_parsed_cell( "wgsox,wgsoy,attributes_json,geom," "geometry_repaired,geometry_issue) " "VALUES(%s,%s,%s,%s,%s,%s,%s,%s,%s,%s," - "ST_GeomFromWKB(%s,4326),%s,%s)", + f"{geom_from_wkb},%s,%s)", cov_rows, ) @@ -770,19 +778,19 @@ def _insert_parsed_cell( point_rows, "INSERT INTO chart_point_geometry" "(feature_id,geom) " - "VALUES(%s,ST_GeomFromWKB(%s,4326))", + f"VALUES(%s,{geom_from_wkb})", ), ( line_rows, "INSERT INTO chart_line_geometry" "(feature_id,geom) " - "VALUES(%s,ST_GeomFromWKB(%s,4326))", + f"VALUES(%s,{geom_from_wkb})", ), ( area_rows, "INSERT INTO chart_area_geometry" "(feature_id,geom) " - "VALUES(%s,ST_GeomFromWKB(%s,4326))", + f"VALUES(%s,{geom_from_wkb})", ), ]: for i in range(0, len(rows), 500): @@ -796,7 +804,7 @@ def _insert_parsed_cell( "INSERT INTO chart_sounding_point" "(feature_id,point_index,raw_x,raw_y,raw_z,depth,geom) " "VALUES(%s,%s,%s,%s,%s,%s," - "ST_GeomFromWKB(%s,4326))", + f"{geom_from_wkb})", sounding_rows[i:i + 1000], ) @@ -881,7 +889,7 @@ def run_import( conn = connect(cfg) try: - assert_mysql57( + assert_mysql_compatible( server_version(conn) ) dataset_id = _dataset( diff --git a/src/cm95pbf/builder.py b/src/cm95pbf/builder.py index 576db77..a162b0a 100644 --- a/src/cm95pbf/builder.py +++ b/src/cm95pbf/builder.py @@ -25,6 +25,11 @@ from shapely.geometry import ( from shapely.ops import transform, unary_union from cm93db.config import DbConfig +from cm93db.database import ( + mysql_compat_mode_for_connection, + as_wkb_sql, + geom_from_text_sql, +) from . import BUILDER_CONTRACT_VERSION from .identity import ( @@ -335,6 +340,18 @@ class Cm95TileBuilder: raise RuntimeError(f"no scale ladder configured for z{z}") return [int(v) for v in raw] + def navigation_aid_scale_ladder_for_zoom(self, z: int) -> list[int]: + ladders = self.policy.get("navigation_aid_scale_ladders", {}) or {} + raw = ladders.get(str(z), ladders.get(z)) + if not raw: + return [] + values = [int(v) for v in raw] + if len(values) != len(set(values)): + raise RuntimeError( + f"duplicate navigation aid scale in z{z}: {values}" + ) + return values + def source_id(self, row: SourceFeature) -> str: return feature_source_key( self.dataset_id, @@ -582,6 +599,23 @@ class Cm95TileBuilder: trace_rows, ) + navigation_aid_scales = self.navigation_aid_scale_ladder_for_zoom( + job.z + ) + if navigation_aid_scales: + self._replace_navigation_marks_with_retention_ladder( + conn=conn, + manifest=manifest, + job=job, + scales=navigation_aid_scales, + query_bbox_wgs=query_bbox_wgs, + expanded_box_m=expanded_box_m, + tile_box_m=tile_box_m, + layer_features=layer_features, + trace_rows=trace_rows, + report=report, + ) + # Delivery label renderability gate: # label-only source layers must contain at least one property that the # frozen NavSea style can actually render. A label feature without @@ -1107,6 +1141,197 @@ class Cm95TileBuilder: regions.append((scale, visible)) return regions + def _replace_navigation_marks_with_retention_ladder( + self, + *, + conn, + manifest: FidManifest, + job: TileJob, + scales: list[int], + query_bbox_wgs: tuple[float, float, float, float], + expanded_box_m, + tile_box_m, + layer_features: dict[str, list[dict[str, Any]]], + trace_rows: list[dict[str, Any]], + report: dict[str, Any], + ) -> None: + default_features = list(layer_features.get("navigation_marks", [])) + default_ids = { + int(feature["id"]) + for feature in default_features + if feature.get("id") is not None + } + + trace_rows[:] = [ + item + for item in trace_rows + if str(item.get("output_layer") or "") != "navigation_marks" + ] + layer_features["navigation_marks"] = [] + + default_composites = dict(report.get("navigation_composites") or {}) + report["navigation_composites"] = { + "groups": 0, + "component_rows": 0, + "standalone_lights": 0, + "multi_light_groups": 0, + "light_presentation_records": 0, + "named_from_direct_source": 0, + "name_conflict_groups": 0, + } + + scale_regions = self._build_scale_regions( + conn, + scales, + query_bbox_wgs, + expanded_box_m, + tile_box_m, + ) + retention = { + "enabled": True, + "zoom": job.z, + "scale_ladder": list(scales), + "default_output_count": len(default_features), + "default_navigation_composites": default_composites, + "retained_output_count": 0, + "added_output_ids": 0, + "removed_output_ids": 0, + "input_component_rows": 0, + "scale_regions": [ + { + "scale": scale, + "coverage_area_fraction_of_tile": ( + float(region.area / tile_box_m.area) + if tile_box_m.area > 0 + else 0.0 + ), + } + for scale, region in scale_regions + if not region.is_empty + ], + } + report["navigation_aid_retention"] = retention + + input_component_rows = 0 + for scale, visible_region in scale_regions: + if visible_region.is_empty: + continue + + rows = self._fetch_scale_navigation_aid_features( + conn, + scale, + query_bbox_wgs, + ) + input_component_rows += len(rows) + + selected: list[SourceFeature] = [] + for row in rows: + geom_m = _to_mercator_geom(row.geom_wgs84) + clipped = _extract_dimensional( + geom_m.intersection(visible_region).intersection(expanded_box_m), + row.geometry_type, + ) + if clipped is None or clipped.is_empty: + continue + row.geom_wgs84 = clipped + selected.append(row) + + if not selected: + continue + + self._s52_context = self._build_s52_context(selected) + consumed: set[int] = set() + nav_outputs = self._compose_navigation_aids( + selected, + manifest, + scale, + consumed, + trace_rows, + ) + for output in nav_outputs: + layer_features[output["layer"]].append(output["feature"]) + + nav_report = report["navigation_composites"] + nav_report["groups"] += sum( + 1 for output in nav_outputs if output.get("composite") + ) + nav_report["component_rows"] += sum( + len(output.get("component_feature_ids", [])) + for output in nav_outputs + if output.get("composite") + ) + nav_report["standalone_lights"] += sum( + 1 for output in nav_outputs if output.get("standalone_light") + ) + nav_report["multi_light_groups"] += sum( + 1 + for output in nav_outputs + if output.get("composite") + and int(output.get("light_record_count") or 0) > 1 + ) + nav_report["light_presentation_records"] += sum( + int(output.get("light_presentation_count") or 0) + for output in nav_outputs + if output.get("composite") + ) + nav_report["named_from_direct_source"] += sum( + 1 + for output in nav_outputs + if output.get("name_source") == "direct_source" + ) + nav_report["name_conflict_groups"] += sum( + 1 + for output in nav_outputs + if output.get("name_conflict") + ) + + for row in selected: + if row.feature_id in consumed: + continue + rule = self.mapping.resolve(row.object_class, row.geometry_type) + if rule is None: + continue + if str(rule.get("target_source_layer") or "") != "navigation_marks": + continue + status = str(rule.get("mapping_status") or "") + if ( + self.profile == "delivery" + and status not in DELIVERY_MAIN_STATUSES + ): + continue + self._append_main_output( + layer_features, + manifest, + row, + rule, + trace_rows, + ) + + unique_features: dict[int, dict[str, Any]] = {} + for feature in layer_features.get("navigation_marks", []): + output_id = int(feature["id"]) + unique_features.setdefault(output_id, feature) + layer_features["navigation_marks"] = list(unique_features.values()) + + seen_trace_ids: set[int] = set() + deduped_traces: list[dict[str, Any]] = [] + for item in trace_rows: + if str(item.get("output_layer") or "") != "navigation_marks": + deduped_traces.append(item) + continue + output_id = int(item["output_id"]) + if output_id in seen_trace_ids: + continue + seen_trace_ids.add(output_id) + deduped_traces.append(item) + trace_rows[:] = deduped_traces + + retained_ids = set(unique_features) + retention["input_component_rows"] = input_component_rows + retention["retained_output_count"] = len(retained_ids) + retention["added_output_ids"] = len(retained_ids - default_ids) + retention["removed_output_ids"] = len(default_ids - retained_ids) + def _coverage_for_scale( self, conn, @@ -1114,6 +1339,8 @@ class Cm95TileBuilder: query_bbox_wgs: tuple[float, float, float, float], expanded_box_m, ): + db_mode = mysql_compat_mode_for_connection(conn) + coverage_wkb_sql = as_wkb_sql(db_mode, "geom") west, south, east, north = query_bbox_wgs with conn.cursor() as cur: cur.execute( @@ -1153,7 +1380,7 @@ class Cm95TileBuilder: f""" SELECT cell_id, - ST_AsWKB(geom) AS geom_wkb + {coverage_wkb_sql} AS geom_wkb FROM chart_coverage WHERE dataset_id=%s AND cell_id IN ({placeholders}) @@ -1198,12 +1425,102 @@ class Cm95TileBuilder: "))" ) + def _fetch_scale_navigation_aid_features( + self, + conn, + scale: int, + bbox_wgs: tuple[float, float, float, float], + ) -> list[SourceFeature]: + db_mode = mysql_compat_mode_for_connection(conn) + feature_wkb_sql = as_wkb_sql(db_mode, "g.geom") + envelope_geom_sql = geom_from_text_sql(db_mode) + west, south, east, north = bbox_wgs + envelope_wkt = self._tile_envelope_wkt(bbox_wgs) + classes = sorted(NAV_COMPONENT_CLASSES) + placeholders = ",".join(["%s"] * len(classes)) + + with conn.cursor() as cur: + cur.execute( + f""" + SELECT + f.id AS feature_id, + f.cell_id, + c.source_relative_path, + f.source_record_index, + f.object_class, + f.geometry_type, + f.native_scale, + f.coverage_id, + f.attributes_json, + {feature_wkb_sql} AS geom_wkb + FROM chart_point_geometry g FORCE INDEX (spx_point_geom) + STRAIGHT_JOIN chart_feature f ON f.id=g.feature_id + JOIN chart_cell c ON c.id=f.cell_id + WHERE MBRIntersects( + g.geom, + {envelope_geom_sql} + ) + AND f.dataset_id=%s + AND f.native_scale=%s + AND f.geometry_type='POINT' + AND f.object_class IN ({placeholders}) + AND f.bbox_max_lon >= %s + AND f.bbox_min_lon <= %s + AND f.bbox_max_lat >= %s + AND f.bbox_min_lat <= %s + """, + ( + envelope_wkt, + self.dataset_id, + scale, + *classes, + west, + east, + south, + north, + ), + ) + rows: list[SourceFeature] = [] + for raw in cur.fetchall(): + geom = _loads_wkb(raw["geom_wkb"]) + if geom is None or geom.is_empty: + continue + rows.append( + SourceFeature( + feature_id=int(raw["feature_id"]), + cell_id=int(raw["cell_id"]), + source_relative_path=str(raw["source_relative_path"]), + source_record_index=int(raw["source_record_index"]), + object_class=str(raw["object_class"]), + geometry_type=str(raw["geometry_type"]), + native_scale=int(raw["native_scale"]), + coverage_id=( + int(raw["coverage_id"]) + if raw["coverage_id"] is not None + else None + ), + attrs=_json_dict(raw["attributes_json"]), + geom_wgs84=geom, + ) + ) + + rows.sort( + key=lambda row: ( + row.source_relative_path, + row.source_record_index, + ) + ) + return rows + def _fetch_scale_features( self, conn, scale: int, bbox_wgs: tuple[float, float, float, float], ) -> list[SourceFeature]: + db_mode = mysql_compat_mode_for_connection(conn) + feature_wkb_sql = as_wkb_sql(db_mode, "g.geom") + envelope_geom_sql = geom_from_text_sql(db_mode) west, south, east, north = bbox_wgs envelope_wkt = self._tile_envelope_wkt(bbox_wgs) rows: list[SourceFeature] = [] @@ -1234,13 +1551,13 @@ class Cm95TileBuilder: f.native_scale, f.coverage_id, f.attributes_json, - ST_AsWKB(g.geom) AS geom_wkb + {feature_wkb_sql} AS geom_wkb FROM {table} g FORCE INDEX ({spatial_index}) STRAIGHT_JOIN chart_feature f ON f.id=g.feature_id JOIN chart_cell c ON c.id=f.cell_id WHERE MBRIntersects( g.geom, - ST_GeomFromText(%s, 4326) + {envelope_geom_sql} ) AND f.dataset_id=%s AND f.native_scale=%s @@ -1304,10 +1621,13 @@ class Cm95TileBuilder: scale: int, bbox_wgs: tuple[float, float, float, float], ) -> list[dict[str, Any]]: + db_mode = mysql_compat_mode_for_connection(conn) + sounding_wkb_sql = as_wkb_sql(db_mode, "s.geom") + envelope_geom_sql = geom_from_text_sql(db_mode) envelope_wkt = self._tile_envelope_wkt(bbox_wgs) with conn.cursor() as cur: cur.execute( - """ + f""" SELECT f.id AS feature_id, f.cell_id, @@ -1317,12 +1637,12 @@ class Cm95TileBuilder: f.attributes_json, s.point_index, s.depth, - ST_AsWKB(s.geom) AS geom_wkb + {sounding_wkb_sql} AS geom_wkb FROM chart_sounding_point s FORCE INDEX (spx_sounding_geom) STRAIGHT_JOIN chart_feature f ON f.id=s.feature_id JOIN chart_cell c ON c.id=f.cell_id WHERE MBRContains( - ST_GeomFromText(%s, 4326), + {envelope_geom_sql}, s.geom ) AND f.dataset_id=%s diff --git a/src/cm95pbf/production.py b/src/cm95pbf/production.py index cf093fb..f3376b1 100644 --- a/src/cm95pbf/production.py +++ b/src/cm95pbf/production.py @@ -21,6 +21,7 @@ from shapely.geometry import box from shapely.ops import transform from cm93db.config import DbConfig +from cm93db.database import mysql_compat_mode_for_connection, as_wkb_sql from . import BUILDER_CONTRACT_VERSION from .builder import Cm95TileBuilder, TileJob @@ -206,6 +207,8 @@ def _connect(db: DbConfig): def _fetch_plan_geometries(conn, dataset_id: int, scales: list[int]): + db_mode = mysql_compat_mode_for_connection(conn) + coverage_wkb_sql = as_wkb_sql(db_mode, "geom") placeholders = ",".join(["%s"] * len(scales)) with conn.cursor() as cur: cur.execute( @@ -234,7 +237,7 @@ def _fetch_plan_geometries(conn, dataset_id: int, scales: list[int]): with conn.cursor() as cur: cur.execute( f""" - SELECT cell_id, ST_AsWKB(geom) AS geom_wkb + SELECT cell_id, {coverage_wkb_sql} AS geom_wkb FROM chart_coverage WHERE dataset_id=%s AND cell_id IN ({ph}) diff --git a/tools/cm93_mysql_final_geometry_audit.py b/tools/cm93_mysql_final_geometry_audit.py index 61edee2..ecc0717 100644 --- a/tools/cm93_mysql_final_geometry_audit.py +++ b/tools/cm93_mysql_final_geometry_audit.py @@ -5,6 +5,8 @@ import argparse, csv, getpass, json, sys from datetime import datetime from pathlib import Path +from cm93db.database import as_wkb_sql, mysql_compat_mode_for_connection + try: import pymysql from pymysql.cursors import DictCursor @@ -33,16 +35,6 @@ def as_bytes(v): if isinstance(v, bytearray): return bytes(v) return v -def wkb_fn(cur, table): - for fn in ("ST_AsBinary","ST_AsWKB","AsBinary","AsWKB"): - try: - cur.execute(f"SELECT {fn}(geom) AS g FROM {table} LIMIT 1") - cur.fetchone() - return fn - except Exception: - pass - raise RuntimeError(f"{table}: 无可用 Geometry->WKB 函数") - def main(): ap = argparse.ArgumentParser(description="CM93 normalization 最终全库 Geometry 验收") ap.add_argument("--host", default="127.0.0.1") @@ -52,6 +44,7 @@ def main(): ap.add_argument("--database", default="cm93_chart") ap.add_argument("--dataset-id", type=int, default=1) ap.add_argument("--bbox-tolerance", type=float, default=5e-6) + ap.add_argument("--expected-cells", type=int, default=None, help="可选：要求 chart_cell 总数匹配，例如当前日本全量为 1492") a = ap.parse_args() if a.password is None: a.password = getpass.getpass("MySQL password: ") @@ -67,8 +60,10 @@ def main(): connect_timeout=20, read_timeout=900 ) + db_mode = mysql_compat_mode_for_connection(conn) + failures = [] - report = {"dataset_id": a.dataset_id, "generated_at": datetime.now().isoformat()} + report = {"dataset_id": a.dataset_id, "generated_at": datetime.now().isoformat(), "mysql_compat_mode": db_mode} try: with conn.cursor() as cur: @@ -90,20 +85,21 @@ def main(): report["counts"] = counts expected = { - "cells": 784, "features": 322465, - "feature_repaired": 55, "coverage_repaired": 3, - "feature_issues": 0, "coverage_issues": 0, - "geotype16_total": 4406, "geotype16_issues": 0, + "feature_issues": 0, + "coverage_issues": 0, + "geotype16_issues": 0, } + if a.expected_cells is not None: + expected["cells"] = a.expected_cells for k,v in expected.items(): if counts[k] != v: failures.append({"check":k,"expected":v,"actual":counts[k]}) - afn = wkb_fn(cur, "chart_area_geometry") + area_wkb = as_wkb_sql(db_mode, "g.geom") cur.execute(f""" SELECT f.id,f.object_class,f.geometry_repaired, f.bbox_min_lon,f.bbox_min_lat,f.bbox_max_lon,f.bbox_max_lat, - {afn}(g.geom) AS wkb + {area_wkb} AS wkb FROM chart_feature f JOIN chart_area_geometry g ON g.feature_id=f.id WHERE f.dataset_id=%s @@ -125,9 +121,9 @@ def main(): except Exception as e: failures.append({"check":"area_wkb","feature_id":r["id"],"error":repr(e)}) - cfn = wkb_fn(cur, "chart_coverage") + coverage_wkb = as_wkb_sql(db_mode, "geom") cur.execute(f""" - SELECT id,cell_id,source_record_index,geometry_repaired,{cfn}(geom) AS wkb + SELECT id,cell_id,source_record_index,geometry_repaired,{coverage_wkb} AS wkb FROM chart_coverage WHERE dataset_id=%s """,(a.dataset_id,)) cov_total=cov_valid=0 @@ -158,6 +154,7 @@ def main(): (outdir/"最终验收报告_v1.json").write_text(json.dumps(report,ensure_ascii=False,indent=2),encoding="utf-8") summary = f"""CM93 normalization 最终全库 Geometry 验收 v1 Dataset ID : {a.dataset_id} +MySQL compatibility : {report['mysql_compat_mode']} Cells : {report['counts']['cells']} Features : {report['counts']['features']} Feature repaired : {report['counts']['feature_repaired']} diff --git a/tools/cm95_bathymetry_build_national_v5.py b/tools/cm95_bathymetry_build_national_v5.py new file mode 100644 index 0000000..4986f4e --- /dev/null +++ b/tools/cm95_bathymetry_build_national_v5.py @@ -0,0 +1,528 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import json +import math +import os +import resource +import shutil +import sys +import time +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import mapbox_vector_tile +import numpy as np +from shapely.geometry import box, mapping +from shapely import union_all + +import cm95_bathymetry_build_v1 as v1 +import cm95_bathymetry_build_v3_dense as v3 +import cm95_bathymetry_build_v5_adaptive as v5 + +ROOT = Path(__file__).resolve().parents[1] +BASE_CONFIG = ROOT / "config" / "bathymetry" / "cm95_bathymetry_north_kyushu_v5_adaptive.yaml" +PLAN_PATH = ROOT / "build" / "production" / "cm95_s52_japan_production_plan_v1.json" +PROBE_PATH = ROOT / "build" / "bathymetry" / "national_z8_probe.json" +OUTPUT_ROOT = ROOT / "build" / "bathymetry" / "cm95-bathymetry-japan-v5-adaptive" +AOI = [122.0, 20.0, 154.5, 46.0] +BUILD_ID = "cm95-bathymetry-japan-v5-adaptive" +def parse_args() -> argparse.Namespace: + p = argparse.ArgumentParser(description="Build nationwide Japan v5 bathymetry PBF by z8 core chunks.") + p.add_argument("--only", default="", help="Comma separated z8 x/y pairs, e.g. 222/102,223/101") + p.add_argument("--limit", type=int, default=0) + p.add_argument("--overlap-m", type=float, default=6000.0) + p.add_argument("--resume", action="store_true") + p.add_argument("--workers", type=int, default=1) + p.add_argument("--shard-count", type=int, default=1) + p.add_argument("--shard-index", type=int, default=0) + p.add_argument("--no-finalize", action="store_true") + p.add_argument("--output-root", type=Path, default=OUTPUT_ROOT) + return p.parse_args() + + +def load_active_tiles() -> list[tuple[int, int, int, int]]: + plan = json.loads(PLAN_PATH.read_text(encoding="utf-8")) + plan_z8 = {(int(j["x"]), int(j["y"])) for j in plan["jobs"] if int(j["z"]) == 8} + probe = json.loads(PROBE_PATH.read_text(encoding="utf-8")) + active = [] + for row in probe: + x = int(float(row["x8"])) + y = int(float(row["y8"])) + shallow = int(row.get("shallow60") or 0) + n = int(row.get("n") or 0) + if shallow > 0 and (x, y) in plan_z8: + active.append((x, y, n, shallow)) + active.sort(key=lambda item: (item[1], item[0])) + return active + + +def mercator_aoi_box(): + west, south = v1.lonlat_to_mercator(AOI[0], AOI[1]) + east, north = v1.lonlat_to_mercator(AOI[2], AOI[3]) + return box(west, south, east, north) + + +def core_box_for_z8(x: int, y: int): + tile = v1.xyz_tile_bounds(8, x, y) + return box(tile.west_m, tile.south_m, tile.east_m, tile.north_m).intersection(mercator_aoi_box()) + + +def bbox_wgs84_for_mercator(geom) -> list[float]: + minx, miny, maxx, maxy = geom.bounds + west, south = v1.mercator_to_lonlat(minx, miny) + east, north = v1.mercator_to_lonlat(maxx, maxy) + return [float(west), float(south), float(east), float(north)] +def descendant_tiles(x8: int, y8: int, z: int): + factor = 1 << (z - 8) + for x in range(x8 * factor, (x8 + 1) * factor): + for y in range(y8 * factor, (y8 + 1) * factor): + yield x, y + + +def build_core_tiles( + contours: dict[float, Any], + cfg: dict[str, Any], + output_root: Path, + x8: int, + y8: int, +) -> dict[str, Any]: + pbf_root = output_root / "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 = BUILD_ID + major_contours = {float(v) for v in cfg.get("major_contours_m", cfg["contour_levels_m"])} + simplify_by_zoom = {int(k): float(v) for k, v in cfg.get("simplify_m_by_zoom", {}).items()} + + tile_count = feature_count = total_bytes = encoder_dropped = 0 + zoom_stats: dict[str, Any] = {} + + for z in range(int(cfg["minzoom"]), int(cfg["maxzoom"]) + 1): + z_tiles = z_features = z_bytes = z_drop = 0 + simplify_m = simplify_by_zoom.get(z, 0.0) + for x, y in descendant_tiles(x8, y8, z): + tile = v1.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 = [] + 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 v1.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 v1.line_survives_mvt_quantization(part, tile.mercator_tuple, extent): + continue + local_index += 1 + features.append({ + "id": v1.feature_id(build_id, z, x, y, depth, local_index, part), + "geometry": mapping(part), + "properties": { + "depth_m": float(depth), + "contour_class": "major" if float(depth) in major_contours else "minor", + "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}, + ) + actual = len(mapbox_vector_tile.decode(data).get(source_layer, {}).get("features", [])) + dropped = max(0, len(features)-actual) + encoder_dropped += dropped + z_drop += dropped + tile_path = pbf_root / str(z) / str(x) / f"{y}.pbf" + tile_path.parent.mkdir(parents=True, exist_ok=True) + if tile_path.exists(): + raise RuntimeError(f"tile ownership collision: {tile_path}") + 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_drop, "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 extract_contours_single_pass_fast( + xs: np.ndarray, + ys: np.ndarray, + grid: np.ndarray, + water_mask, + cfg: dict[str, Any], +) -> tuple[dict[float, Any], dict[str, Any]]: + levels = sorted(float(value) for value in cfg["contour_levels_m"]) + segments: dict[float, list[Any]] = {level: [] for level in levels} + fallback_counts: dict[float, int] = {level: 0 for level in levels} + + v00 = grid[:-1, :-1] + v10 = grid[:-1, 1:] + v11 = grid[1:, 1:] + v01 = grid[1:, :-1] + + finite_count = ( + np.isfinite(v00).astype(np.uint8) + + np.isfinite(v10).astype(np.uint8) + + np.isfinite(v11).astype(np.uint8) + + np.isfinite(v01).astype(np.uint8) + ) + lo = np.fmin(np.fmin(v00, v10), np.fmin(v11, v01)) + hi = np.fmax(np.fmax(v00, v10), np.fmax(v11, v01)) + candidate_mask = ( + (finite_count >= 3) + & np.isfinite(lo) + & np.isfinite(hi) + & (lo != hi) + & (hi >= levels[0]) + & (lo <= levels[-1]) + ) + candidate_cells = np.argwhere(candidate_mask) + + import bisect + import math + + for iy, ix in candidate_cells: + values = ( + float(grid[iy, ix]), + float(grid[iy, ix + 1]), + float(grid[iy + 1, ix + 1]), + float(grid[iy + 1, ix]), + ) + finite_values = [value for value in values if math.isfinite(value)] + cell_lo = min(finite_values) + cell_hi = max(finite_values) + start = bisect.bisect_left(levels, cell_lo) + stop = bisect.bisect_right(levels, cell_hi) + finite_n = len(finite_values) + y0 = float(ys[iy]) + y1 = float(ys[iy + 1]) + x0 = float(xs[ix]) + x1 = float(xs[ix + 1]) + for level in levels[start:stop]: + cell_segments = v3._cell_contours( + x0, + x1, + y0, + y1, + values, + level, + ) + segments[level].extend(cell_segments) + if finite_n == 3: + fallback_counts[level] += len(cell_segments) + + contours: dict[float, Any] = {} + stats: dict[str, Any] = {} + min_segment = float(cfg["min_segment_length_m"]) + simplify_m = float(cfg.get("regional_simplify_m", 0)) + smooth_iterations = int(cfg.get("chaikin_iterations", 0)) + + import cm95_bathymetry_build_v2_smooth as v2 + + for level in levels: + raw = segments[level] + if not raw: + stats[str(level)] = { + "raw_segments": 0, + "triangle_fallback_segments": 0, + "line_parts": 0, + "length_m": 0.0, + } + continue + + from shapely import line_merge + merged = line_merge(union_all(raw)) + parts: list[Any] = [] + for part in v1.line_parts(merged): + if part.length < min_segment: + continue + smoothed = v2.chaikin_smooth(part, smooth_iterations) + if simplify_m > 0: + smoothed = smoothed.simplify( + simplify_m, + preserve_topology=False, + ) + clipped = smoothed.intersection(water_mask) + parts.extend( + candidate + for candidate in v1.line_parts(clipped) + if candidate.length >= min_segment + ) + + if not parts: + stats[str(level)] = { + "raw_segments": len(raw), + "triangle_fallback_segments": fallback_counts[level], + "line_parts": 0, + "length_m": 0.0, + } + continue + + geom = union_all(parts) + contours[level] = geom + final_parts = v1.line_parts(geom) + stats[str(level)] = { + "raw_segments": len(raw), + "triangle_fallback_segments": fallback_counts[level], + "line_parts": len(final_parts), + "length_m": round( + sum(part.length for part in final_parts), + 3, + ), + } + + stats["_performance"] = { + "candidate_cells": int(len(candidate_cells)), + "all_cells": int((len(ys) - 1) * (len(xs) - 1)), + } + return contours, stats + + +def clip_contours_to_core(contours: dict[float, Any], core_box) -> dict[float, Any]: + out: dict[float, Any] = {} + for depth, geom in contours.items(): + clipped = geom.intersection(core_box) + parts = [part for part in v1.line_parts(clipped) if not part.is_empty and part.length > 0] + if parts: + out[depth] = union_all(parts) + return out + + +def build_chunk( + x8: int, + y8: int, + point_hint: int, + shallow_hint: int, + overlap_m: float, + output_root: Path, + resume: bool, +) -> dict[str, Any]: + chunk_dir = output_root / "chunks" + chunk_dir.mkdir(parents=True, exist_ok=True) + report_path = chunk_dir / f"z8_{x8}_{y8}.json" + if resume and report_path.exists(): + old = json.loads(report_path.read_text(encoding="utf-8")) + if old.get("status") == "PASS": + return old + + started = time.monotonic() + started_at = datetime.now(timezone.utc) + cfg = v1.load_config(BASE_CONFIG) + cfg["build_id"] = BUILD_ID + cfg["source_kind"] = "derived_sounding_adaptive_tin_v5" + core_box = core_box_for_z8(x8, y8) + if core_box.is_empty: + return {"status": "SKIP", "x8": x8, "y8": y8, "reason": "outside_aoi"} + + expanded = core_box.buffer(float(overlap_m), cap_style=3, join_style=2) + cfg["bbox_wgs84"] = bbox_wgs84_for_mercator(expanded) + + conn = v1.connect_db() + try: + points, sounding_stats = v1.fetch_soundings(conn, cfg) + if len(points) < 3: + result = {"status": "SKIP", "x8": x8, "y8": y8, "reason": "fewer_than_3_soundings", "duration_seconds": round(time.monotonic()-started,3)} + report_path.write_text(json.dumps(result, ensure_ascii=False, indent=2), encoding="utf-8") + return result + water_mask, water_rows = v1.fetch_water_mask(conn, cfg) + finally: + conn.close() + + xs, ys, raw_grid, water_grid, local_spacing_grid, surface_stats = v5.build_adaptive_tin_surface(points, water_mask, cfg) + repaired_grid, gap_stats = v3.repair_small_water_gaps(raw_grid, water_grid, cfg) + newly_filled = (~__import__("numpy").isfinite(raw_grid)) & __import__("numpy").isfinite(repaired_grid) + local_spacing_grid[newly_filled] = float(cfg["smooth_spacing_high_m"]) + smooth_grid, smooth_stats = v5.adaptive_smooth_surface(repaired_grid, water_grid, local_spacing_grid, cfg) + contours, contour_stats = extract_contours_single_pass_fast(xs, ys, smooth_grid, water_mask, cfg) + core_contours = clip_contours_to_core(contours, core_box) + tile_stats = build_core_tiles(core_contours, cfg, output_root, x8, y8) if core_contours else { + "pbf_tiles": 0, "tile_feature_instances": 0, "encoder_dropped_features": 0, "pbf_bytes": 0, "zooms": {} + } + + duration = round(time.monotonic()-started, 3) + peak_kib = int(resource.getrusage(resource.RUSAGE_SELF).ru_maxrss) + result = { + "status": "PASS", + "x8": x8, + "y8": y8, + "point_hint": point_hint, + "shallow60_hint": shallow_hint, + "input_bbox_wgs84": cfg["bbox_wgs84"], + "core_bbox_wgs84": bbox_wgs84_for_mercator(core_box), + "overlap_m": overlap_m, + "started_at": started_at.isoformat(), + "duration_seconds": duration, + "peak_rss_kib_process": peak_kib, + "sounding_stats": sounding_stats, + "water_area_rows": water_rows, + "surface_stats": surface_stats, + "gap_repair_stats": gap_stats, + "surface_smooth_stats": smooth_stats, + "contour_levels_generated": sorted(float(k) for k in core_contours.keys()), + "contour_stats": contour_stats, + "tile_stats": tile_stats, + } + report_path.write_text(json.dumps(result, ensure_ascii=False, indent=2), encoding="utf-8") + return result +def aggregate(output_root: Path, active_count: int, started_at: datetime, elapsed: float) -> dict[str, Any]: + reports = [] + for path in sorted((output_root / "chunks").glob("z8_*.json")): + reports.append(json.loads(path.read_text(encoding="utf-8"))) + passed = [r for r in reports if r.get("status") == "PASS"] + skipped = [r for r in reports if r.get("status") == "SKIP"] + tile_stats = { + "pbf_tiles": sum(int(r.get("tile_stats", {}).get("pbf_tiles", 0)) for r in passed), + "tile_feature_instances": sum(int(r.get("tile_stats", {}).get("tile_feature_instances", 0)) for r in passed), + "encoder_dropped_features": sum(int(r.get("tile_stats", {}).get("encoder_dropped_features", 0)) for r in passed), + "pbf_bytes": sum(int(r.get("tile_stats", {}).get("pbf_bytes", 0)) for r in passed), + } + actual_files = list((output_root / "pbf").rglob("*.pbf")) if (output_root / "pbf").exists() else [] + actual_bytes = sum(p.stat().st_size for p in actual_files) + return { + "schema": "navsea.cm95.bathymetry-national-result/v1", + "build_id": BUILD_ID, + "aoi_wgs84": AOI, + "algorithm_baseline": "v5-adaptive", + "active_z8_expected": active_count, + "chunk_reports": len(reports), + "chunks_pass": len(passed), + "chunks_skip": len(skipped), + "started_at": started_at.isoformat(), + "finished_at": datetime.now(timezone.utc).isoformat(), + "elapsed_seconds": round(elapsed, 3), + "sum_chunk_seconds": round(sum(float(r.get("duration_seconds", 0)) for r in reports), 3), + "max_chunk_seconds": round(max([float(r.get("duration_seconds", 0)) for r in reports] or [0]), 3), + "tile_stats_sum": tile_stats, + "actual_pbf_files": len(actual_files), + "actual_pbf_bytes": actual_bytes, + "disk_usage_bytes": sum(p.stat().st_size for p in output_root.rglob("*") if p.is_file()), + } + + +def main() -> int: + args = parse_args() + output_root = args.output_root.resolve() + output_root.mkdir(parents=True, exist_ok=True) + active = load_active_tiles() + + only = set() + if args.only.strip(): + for token in args.only.split(","): + x, y = token.strip().split("/") + only.add((int(x), int(y))) + active = [item for item in active if (item[0], item[1]) in only] + if args.limit > 0: + active = active[: args.limit] + + if args.shard_count < 1: + raise SystemExit("--shard-count must be >= 1") + if args.shard_index < 0 or args.shard_index >= args.shard_count: + raise SystemExit("--shard-index must satisfy 0 <= index < count") + active = [ + item + for index, item in enumerate(active) + if index % args.shard_count == args.shard_index + ] + + started_at = datetime.now(timezone.utc) + started = time.monotonic() + print(f"NATIONAL_BUILD {BUILD_ID}", flush=True) + print(f"ACTIVE_Z8 {len(active)}", flush=True) + print(f"OUTPUT {output_root}", flush=True) + + workers = max(1, int(args.workers)) + print(f"WORKERS {workers}", flush=True) + + if workers == 1: + for index, (x8, y8, n, shallow) in enumerate(active, start=1): + t0 = time.monotonic() + print(f"CHUNK_START {index}/{len(active)} z8/{x8}/{y8} points={n} shallow60={shallow}", flush=True) + result = build_chunk(x8, y8, n, shallow, args.overlap_m, output_root, args.resume) + print( + f"CHUNK_DONE {index}/{len(active)} z8/{x8}/{y8} status={result.get('status')} " + f"seconds={result.get('duration_seconds', round(time.monotonic()-t0,3))} " + f"tiles={result.get('tile_stats',{}).get('pbf_tiles',0)} " + f"bytes={result.get('tile_stats',{}).get('pbf_bytes',0)}", + flush=True, + ) + else: + from concurrent.futures import ProcessPoolExecutor, as_completed + + completed = 0 + with ProcessPoolExecutor(max_workers=workers) as executor: + futures = {} + for x8, y8, n, shallow in active: + print(f"CHUNK_QUEUE z8/{x8}/{y8} points={n} shallow60={shallow}", flush=True) + future = executor.submit( + build_chunk, + x8, + y8, + n, + shallow, + args.overlap_m, + output_root, + args.resume, + ) + futures[future] = (x8, y8, n, shallow) + + for future in as_completed(futures): + x8, y8, n, shallow = futures[future] + completed += 1 + try: + result = future.result() + except Exception as exc: + print( + f"CHUNK_FAIL {completed}/{len(active)} z8/{x8}/{y8} " + f"error={type(exc).__name__}:{exc}", + flush=True, + ) + raise + print( + f"CHUNK_DONE {completed}/{len(active)} z8/{x8}/{y8} " + f"status={result.get('status')} " + f"seconds={result.get('duration_seconds',0)} " + f"tiles={result.get('tile_stats',{}).get('pbf_tiles',0)} " + f"bytes={result.get('tile_stats',{}).get('pbf_bytes',0)}", + flush=True, + ) + + if not args.no_finalize: + manifest = aggregate( + output_root, + len(active), + started_at, + time.monotonic()-started, + ) + (output_root / "build_manifest.json").write_text( + json.dumps(manifest, ensure_ascii=False, indent=2), + encoding="utf-8", + ) + print( + "NATIONAL_DONE", + json.dumps(manifest, ensure_ascii=False), + flush=True, + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/cm95_bathymetry_build_v1.py b/tools/cm95_bathymetry_build_v1.py new file mode 100644 index 0000000..c466d4c --- /dev/null +++ b/tools/cm95_bathymetry_build_v1.py @@ -0,0 +1,603 @@ +#!/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 cm93db.database import as_wkb_sql, geom_from_text_sql, mysql_compat_mode_for_connection +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"]) + mode = mysql_compat_mode_for_connection(conn) + sql = f""" + SELECT {as_wkb_sql(mode, "s.geom")} AS geom_wkb, 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({geom_from_text_sql(mode)}, 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 + point = wkb.loads(bytes(row["geom_wkb"])) + if point.is_empty: + continue + mx, my = lonlat_to_mercator(float(point.x), float(point.y)) + 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)) + mode = mysql_compat_mode_for_connection(conn) + sql = f""" + SELECT {as_wkb_sql(mode, "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, {geom_from_text_sql(mode)}) + 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"]) + major_contours = { + float(v) + for v in cfg.get( + "major_contours_m", + cfg.get("contour_levels_m", []), + ) + } + 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), + "contour_class": ( + "major" + if float(depth) in major_contours + else "minor" + ), + "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 = [] + major_contours = { + float(v) + for v in cfg.get( + "major_contours_m", + cfg.get("contour_levels_m", []), + ) + } + 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), + "contour_class": ( + "major" + if float(depth) in major_contours + else "minor" + ), + "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", + "contour_class": "String", + "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/tools/cm95_bathymetry_build_v2_smooth.py b/tools/cm95_bathymetry_build_v2_smooth.py new file mode 100644 index 0000000..055ec65 --- /dev/null +++ b/tools/cm95_bathymetry_build_v2_smooth.py @@ -0,0 +1,502 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import json +import math +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import numpy as np +from shapely import contains_xy, line_merge, union_all +from shapely.geometry import LineString + +import cm95_bathymetry_build_v1 as v1 + +ROOT = Path(__file__).resolve().parents[1] +DEFAULT_CONFIG = ( + ROOT + / "config" + / "bathymetry" + / "cm95_bathymetry_north_kyushu_v2_smooth.yaml" +) + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Build smooth independent bathymetry PBF using local IDW grid interpolation." + ) + parser.add_argument("--config", type=Path, default=DEFAULT_CONFIG) + parser.add_argument("--output-dir", type=Path, default=None) + return parser.parse_args() + + +def _grid_axis(start: float, stop: float, spacing: float) -> np.ndarray: + count = int(math.ceil((stop - start) / spacing)) + 1 + return start + np.arange(count, dtype=np.float64) * spacing + + +def _bin_key( + x: float, + y: float, + minx: float, + miny: float, + bin_size: float, +) -> tuple[int, int]: + return ( + int(math.floor((x - minx) / bin_size)), + int(math.floor((y - miny) / bin_size)), + ) + +def build_idw_grid( + points: dict[tuple[float, float], float], + water_mask, + cfg: dict[str, Any], +) -> tuple[np.ndarray, np.ndarray, np.ndarray, dict[str, Any]]: + west, south, east, north = [float(v) for v in cfg["bbox_wgs84"]] + minx, miny = v1.lonlat_to_mercator(west, south) + maxx, maxy = v1.lonlat_to_mercator(east, north) + + spacing = float(cfg["grid_spacing_m"]) + radius = float(cfg["idw_search_radius_m"]) + max_nearest = float(cfg["idw_max_nearest_m"]) + min_neighbors = int(cfg["idw_min_neighbors"]) + max_neighbors = int(cfg["idw_max_neighbors"]) + power = float(cfg["idw_power"]) + + xs = _grid_axis(minx, maxx, spacing) + ys = _grid_axis(miny, maxy, spacing) + grid = np.full((len(ys), len(xs)), np.nan, dtype=np.float64) + + point_rows = np.array( + [(x, y, depth) for (x, y), depth in points.items()], + dtype=np.float64, + ) + bins: dict[tuple[int, int], list[int]] = {} + + bin_size = radius + for index, (x, y, _depth) in enumerate(point_rows): + bins.setdefault( + _bin_key(float(x), float(y), minx, miny, bin_size), + [], + ).append(index) + + radius2 = radius * radius + max_nearest2 = max_nearest * max_nearest + interpolated = 0 + rejected_no_neighbors = 0 + rejected_nearest = 0 + water_nodes = 0 + + for iy, y in enumerate(ys): + row_mask = contains_xy( + water_mask, + xs, + np.full(xs.shape, y, dtype=np.float64), + ) + water_indexes = np.flatnonzero(row_mask) + water_nodes += int(len(water_indexes)) + + for ix in water_indexes: + x = float(xs[ix]) + bx, by = _bin_key(x, float(y), minx, miny, bin_size) + candidate_indexes: list[int] = [] + + for ox in (-1, 0, 1): + for oy in (-1, 0, 1): + candidate_indexes.extend( + bins.get((bx + ox, by + oy), ()) + ) + + if not candidate_indexes: + rejected_no_neighbors += 1 + continue + + candidates = point_rows[candidate_indexes] + dx = candidates[:, 0] - x + dy = candidates[:, 1] - float(y) + dist2 = dx * dx + dy * dy + inside = np.flatnonzero(dist2 <= radius2) + + if len(inside) < min_neighbors: + rejected_no_neighbors += 1 + continue + + inside = inside[np.argsort(dist2[inside])] + if float(dist2[inside[0]]) > max_nearest2: + rejected_nearest += 1 + continue + + inside = inside[:max_neighbors] + chosen_d2 = dist2[inside] + chosen_depths = candidates[inside, 2] + + if float(chosen_d2[0]) < 1.0: + value = float(chosen_depths[0]) + + else: + weights = 1.0 / np.power( + np.maximum(chosen_d2, 1.0), + power / 2.0, + ) + value = float( + np.sum(weights * chosen_depths) + / np.sum(weights) + ) + + grid[iy, ix] = value + interpolated += 1 + + stats = { + "grid_width": int(len(xs)), + "grid_height": int(len(ys)), + "grid_nodes": int(grid.size), + "water_nodes": water_nodes, + "interpolated_nodes": interpolated, + "rejected_no_neighbors": rejected_no_neighbors, + "rejected_nearest": rejected_nearest, + "finite_ratio": round( + interpolated / float(grid.size), + 6, + ), + } + return xs, ys, grid, stats + + +def _edge_intersection( + p1: tuple[float, float], + p2: tuple[float, float], + + v1_value: float, + v2_value: float, + level: float, +) -> tuple[float, float] | None: + dv = v2_value - v1_value + if abs(dv) < 1e-12: + return None + t = (level - v1_value) / dv + if t < 0.0 or t > 1.0: + return None + return ( + p1[0] + t * (p2[0] - p1[0]), + p1[1] + t * (p2[1] - p1[1]), + ) + + +def _march_cell( + x0: float, + x1: float, + y0: float, + y1: float, + values: tuple[float, float, float, float], + level: float, +) -> list[LineString]: + v0, v1_value, v2, v3 = values + points = ( + (x0, y0), + (x1, y0), + (x1, y1), + (x0, y1), + ) + + edge_vertices = ((0, 1), (1, 2), (2, 3), (3, 0)) + hits: dict[int, tuple[float, float]] = {} + + for edge_index, (a, b) in enumerate(edge_vertices): + va, vb = values[a], values[b] + if (va < level and vb < level) or (va > level and vb > level): + continue + hit = _edge_intersection(points[a], points[b], va, vb, level) + if hit is not None: + hits[edge_index] = hit + + if len(hits) == 2: + ordered = list(hits.values()) + if ordered[0] != ordered[1]: + return [LineString(ordered)] + return [] + + if len(hits) != 4: + return [] + + case_index = 0 + for bit, value in enumerate(values): + if value >= level: + case_index |= 1 << bit + + center_high = sum(values) * 0.25 >= level + if case_index == 5: + pairs = ((0, 1), (2, 3)) if center_high else ((3, 0), (1, 2)) + elif case_index == 10: + + pairs = ((3, 0), (1, 2)) if center_high else ((0, 1), (2, 3)) + else: + pairs = ((0, 1), (2, 3)) + + return [ + LineString((hits[a], hits[b])) + for a, b in pairs + if hits[a] != hits[b] + ] + + +def chaikin_smooth(line: LineString, iterations: int) -> LineString: + coords = [(float(x), float(y)) for x, y, *rest in line.coords] + if len(coords) < 3 or iterations <= 0: + return line + + closed = coords[0] == coords[-1] + if closed: + coords = coords[:-1] + + for _ in range(iterations): + if len(coords) < 2: + break + output: list[tuple[float, float]] = [] + if not closed: + output.append(coords[0]) + pair_count = len(coords) if closed else len(coords) - 1 + for index in range(pair_count): + p0 = coords[index] + p1 = coords[(index + 1) % len(coords)] + + q = ( + 0.75 * p0[0] + 0.25 * p1[0], + 0.75 * p0[1] + 0.25 * p1[1], + ) + r = ( + 0.25 * p0[0] + 0.75 * p1[0], + 0.25 * p0[1] + 0.75 * p1[1], + ) + output.extend((q, r)) + if not closed: + output.append(coords[-1]) + coords = output + + if closed and coords: + coords.append(coords[0]) + return LineString(coords) + + +def extract_smooth_contours( + xs: np.ndarray, + ys: np.ndarray, + grid: np.ndarray, + water_mask, + cfg: dict[str, Any], +) -> tuple[dict[float, Any], dict[str, Any]]: + levels = sorted(float(v) for v in cfg["contour_levels_m"]) + min_segment = float(cfg["min_segment_length_m"]) + simplify_m = float(cfg.get("regional_simplify_m", 0)) + smooth_iterations = int(cfg.get("chaikin_iterations", 0)) + + contours: dict[float, Any] = {} + stats: dict[str, Any] = {} + + for level in levels: + raw_segments: list[LineString] = [] + for iy in range(len(ys) - 1): + y0 = float(ys[iy]) + y1 = float(ys[iy + 1]) + for ix in range(len(xs) - 1): + values = ( + float(grid[iy, ix]), + float(grid[iy, ix + 1]), + float(grid[iy + 1, ix + 1]), + float(grid[iy + 1, ix]), + ) + if not all(math.isfinite(value) for value in values): + continue + lo, hi = min(values), max(values) + if level < lo or level > hi or lo == hi: + continue + raw_segments.extend( + _march_cell( + float(xs[ix]), + float(xs[ix + 1]), + y0, + y1, + values, + level, + ) + ) + + if not raw_segments: + stats[str(level)] = { + "raw_segments": 0, + "line_parts": 0, + "length_m": 0.0, + } + continue + + merged = line_merge(union_all(raw_segments)) + parts = [] + for part in v1.line_parts(merged): + if part.length < min_segment: + continue + smooth = chaikin_smooth(part, smooth_iterations) + if simplify_m > 0: + smooth = smooth.simplify( + simplify_m, + preserve_topology=False, + ) + clipped = smooth.intersection(water_mask) + parts.extend( + candidate + for candidate in v1.line_parts(clipped) + if candidate.length >= min_segment + ) + + if not parts: + stats[str(level)] = { + "raw_segments": len(raw_segments), + "line_parts": 0, + "length_m": 0.0, + } + continue + + geom = union_all(parts) + contours[level] = geom + final_parts = v1.line_parts(geom) + stats[str(level)] = { + "raw_segments": len(raw_segments), + "line_parts": len(final_parts), + "length_m": round( + sum(part.length for part in final_parts), + 3, + ), + } + + return contours, stats + + +def main() -> int: + args = parse_args() + config_path = args.config.resolve() + cfg = v1.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']}", flush=True) + print(f"BBOX {cfg['bbox_wgs84']}", flush=True) + print(f"OUTPUT {output_dir}", flush=True) + + conn = v1.connect_db() + try: + points, sounding_stats = v1.fetch_soundings(conn, cfg) + print( + "SOUNDINGS", + json.dumps(sounding_stats, ensure_ascii=False), + flush=True, + ) + water_mask, water_rows = v1.fetch_water_mask(conn, cfg) + print(f"WATER_AREAS {water_rows}", flush=True) + finally: + conn.close() + + xs, ys, grid, surface_stats = build_idw_grid( + points, + water_mask, + cfg, + ) + print( + "SURFACE", + json.dumps(surface_stats, ensure_ascii=False), + flush=True, + ) + + contours, contour_stats = extract_smooth_contours( + xs, + ys, + grid, + water_mask, + cfg, + ) + + if not contours: + raise RuntimeError("no smooth contours generated") + print( + "CONTOURS", + json.dumps(contour_stats, ensure_ascii=False), + flush=True, + ) + + debug_feature_count = v1.write_debug_geojson( + contours, + cfg, + output_dir, + ) + v1.write_tilejson(cfg, output_dir) + tile_stats = v1.build_tiles(contours, cfg, output_dir) + print( + "TILES", + json.dumps(tile_stats, ensure_ascii=False), + flush=True, + ) + + finished = datetime.now(timezone.utc) + manifest = { + "schema": "navsea.cm95.bathymetry-build-result/v2", + "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"]), + "water_object_classes": [ + str(v) for v in cfg["water_object_classes"] + ], + "surface_method": "local_idw_regular_grid", + "surface_parameters": { + key: cfg[key] + for key in ( + "grid_spacing_m", + "idw_search_radius_m", + "idw_max_nearest_m", + "idw_min_neighbors", + "idw_max_neighbors", + "idw_power", + "chaikin_iterations", + ) + }, + + "sounding_stats": sounding_stats, + "water_area_rows": water_rows, + "surface_stats": surface_stats, + "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", + flush=True, + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/cm95_bathymetry_build_v3_dense.py b/tools/cm95_bathymetry_build_v3_dense.py new file mode 100644 index 0000000..0cbd4a3 --- /dev/null +++ b/tools/cm95_bathymetry_build_v3_dense.py @@ -0,0 +1,515 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import json +import math +from collections import deque +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import numpy as np +from shapely import contains_xy, line_merge, union_all +from shapely.geometry import LineString + +import cm95_bathymetry_build_v1 as v1 +import cm95_bathymetry_build_v2_smooth as v2 + +ROOT = Path(__file__).resolve().parents[1] +DEFAULT_CONFIG = ( + ROOT + / "config" + / "bathymetry" + / "cm95_bathymetry_north_kyushu_v3_dense.yaml" +) +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Build dense smooth bathymetry PBF with conservative water-gap repair." + ) + parser.add_argument("--config", type=Path, default=DEFAULT_CONFIG) + parser.add_argument("--output-dir", type=Path, default=None) + return parser.parse_args() + + +def make_water_grid_mask( + xs: np.ndarray, + ys: np.ndarray, + water_mask, +) -> np.ndarray: + mask = np.zeros((len(ys), len(xs)), dtype=bool) + for iy, y in enumerate(ys): + mask[iy, :] = contains_xy( + water_mask, + xs, + np.full(xs.shape, float(y), dtype=np.float64), + ) + return mask + + +_NEIGHBORS4 = ((-1, 0), (1, 0), (0, -1), (0, 1)) +_NEIGHBORS8 = ( + (-1, -1), (-1, 0), (-1, 1), + (0, -1), (0, 1), + (1, -1), (1, 0), (1, 1), +) +def _collect_component( + start_y: int, + start_x: int, + missing: np.ndarray, + visited: np.ndarray, +) -> list[tuple[int, int]]: + height, width = missing.shape + queue = deque([(start_y, start_x)]) + visited[start_y, start_x] = True + component: list[tuple[int, int]] = [] + + while queue: + iy, ix = queue.popleft() + component.append((iy, ix)) + for dy, dx in _NEIGHBORS4: + ny, nx = iy + dy, ix + dx + if ny < 0 or nx < 0 or ny >= height or nx >= width: + continue + if visited[ny, nx] or not missing[ny, nx]: + continue + visited[ny, nx] = True + queue.append((ny, nx)) + return component + + +def _component_touches_nonwater( + component: list[tuple[int, int]], + water: np.ndarray, +) -> bool: + height, width = water.shape + for iy, ix in component: + if iy == 0 or ix == 0 or iy == height - 1 or ix == width - 1: + return True + for dy, dx in _NEIGHBORS8: + ny, nx = iy + dy, ix + dx + if not water[ny, nx]: + return True + return False +def repair_small_water_gaps( + grid: np.ndarray, + water: np.ndarray, + cfg: dict[str, Any], +) -> tuple[np.ndarray, dict[str, Any]]: + repaired = grid.copy() + missing = water & ~np.isfinite(repaired) + visited = np.zeros(missing.shape, dtype=bool) + + max_nodes = int(cfg["gap_fill_max_component_nodes"]) + radius = int(cfg["gap_fill_radius_cells"]) + min_neighbors = int(cfg["gap_fill_min_neighbors"]) + + components_total = 0 + components_small = 0 + components_repaired = 0 + components_skipped_land = 0 + components_skipped_large = 0 + nodes_filled = 0 + + height, width = grid.shape + for iy in range(height): + for ix in range(width): + if visited[iy, ix] or not missing[iy, ix]: + continue + component = _collect_component(iy, ix, missing, visited) + components_total += 1 + + if len(component) > max_nodes: + components_skipped_large += 1 + continue + components_small += 1 + + if _component_touches_nonwater(component, water): + components_skipped_land += 1 + continue + pending: list[tuple[int, int, float]] = [] + component_set = set(component) + + for cy, cx in component: + samples: list[tuple[float, float]] = [] + y0 = max(0, cy - radius) + y1 = min(height - 1, cy + radius) + x0 = max(0, cx - radius) + x1 = min(width - 1, cx + radius) + + for ny in range(y0, y1 + 1): + for nx in range(x0, x1 + 1): + if (ny, nx) in component_set: + continue + value = repaired[ny, nx] + if not water[ny, nx] or not math.isfinite(float(value)): + continue + dy = ny - cy + dx = nx - cx + dist2 = float(dx * dx + dy * dy) + if dist2 <= 0.0 or dist2 > float(radius * radius): + continue + samples.append((dist2, float(value))) + + if len(samples) < min_neighbors: + pending = [] + break + + samples.sort(key=lambda item: item[0]) + samples = samples[:16] + weights = np.array( + [1.0 / max(item[0], 1.0) for item in samples], + dtype=np.float64, + ) + values = np.array( + [item[1] for item in samples], + dtype=np.float64, + ) + pending.append( + ( + cy, + cx, + float(np.sum(weights * values) / np.sum(weights)), + ) + ) + if len(pending) != len(component): + continue + + for cy, cx, value in pending: + repaired[cy, cx] = value + nodes_filled += 1 + components_repaired += 1 + + return repaired, { + "missing_water_nodes_before": int(np.count_nonzero(missing)), + "missing_water_nodes_after": int( + np.count_nonzero(water & ~np.isfinite(repaired)) + ), + "components_total": components_total, + "components_small": components_small, + "components_repaired": components_repaired, + "components_skipped_land": components_skipped_land, + "components_skipped_large": components_skipped_large, + "nodes_filled": nodes_filled, + } + + +def _triangle_contour( + points: list[tuple[float, float]], + values: list[float], + level: float, +) -> list[LineString]: + hits: list[tuple[float, float]] = [] + for a, b in ((0, 1), (1, 2), (2, 0)): + va, vb = values[a], values[b] + if (va < level and vb < level) or (va > level and vb > level): + continue + hit = v2._edge_intersection( + points[a], + points[b], + va, + vb, + level, + ) + if hit is not None: + hits.append(hit) + unique: list[tuple[float, float]] = [] + for hit in hits: + if not any( + math.hypot(hit[0] - old[0], hit[1] - old[1]) < 1e-6 + for old in unique + ): + unique.append(hit) + + if len(unique) != 2 or unique[0] == unique[1]: + return [] + return [LineString(unique)] + + +def _cell_contours( + x0: float, + x1: float, + y0: float, + y1: float, + values: tuple[float, float, float, float], + level: float, +) -> list[LineString]: + finite = [math.isfinite(value) for value in values] + count = sum(finite) + + if count == 4: + return v2._march_cell( + x0, + x1, + y0, + y1, + values, + level, + ) + if count != 3: + return [] + + all_points = [ + (x0, y0), + (x1, y0), + (x1, y1), + (x0, y1), + ] + indexes = [index for index, ok in enumerate(finite) if ok] + triangle_points = [all_points[index] for index in indexes] + triangle_values = [values[index] for index in indexes] + if level < min(triangle_values) or level > max(triangle_values): + return [] + return _triangle_contour( + triangle_points, + triangle_values, + level, + ) +def extract_dense_contours( + xs: np.ndarray, + ys: np.ndarray, + grid: np.ndarray, + water_mask, + cfg: dict[str, Any], +) -> tuple[dict[float, Any], dict[str, Any]]: + levels = sorted(float(value) for value in cfg["contour_levels_m"]) + min_segment = float(cfg["min_segment_length_m"]) + simplify_m = float(cfg.get("regional_simplify_m", 0)) + smooth_iterations = int(cfg.get("chaikin_iterations", 0)) + + contours: dict[float, Any] = {} + stats: dict[str, Any] = {} + + for level in levels: + raw_segments: list[LineString] = [] + triangle_fallback_segments = 0 + + for iy in range(len(ys) - 1): + y0 = float(ys[iy]) + y1 = float(ys[iy + 1]) + for ix in range(len(xs) - 1): + values = ( + float(grid[iy, ix]), + float(grid[iy, ix + 1]), + float(grid[iy + 1, ix + 1]), + float(grid[iy + 1, ix]), + ) + finite_count = sum(math.isfinite(value) for value in values) + if finite_count < 3: + continue + + finite_values = [value for value in values if math.isfinite(value)] + if level < min(finite_values) or level > max(finite_values): + continue + + segments = _cell_contours( + float(xs[ix]), + float(xs[ix + 1]), + y0, + y1, + values, + level, + ) + if finite_count == 3: + triangle_fallback_segments += len(segments) + raw_segments.extend(segments) + if not raw_segments: + stats[str(level)] = { + "raw_segments": 0, + "triangle_fallback_segments": 0, + "line_parts": 0, + "length_m": 0.0, + } + continue + + merged = line_merge(union_all(raw_segments)) + parts: list[LineString] = [] + for part in v1.line_parts(merged): + if part.length < min_segment: + continue + smoothed = v2.chaikin_smooth( + part, + smooth_iterations, + ) + if simplify_m > 0: + smoothed = smoothed.simplify( + simplify_m, + preserve_topology=False, + ) + clipped = smoothed.intersection(water_mask) + parts.extend( + candidate + for candidate in v1.line_parts(clipped) + if candidate.length >= min_segment + ) + + if not parts: + stats[str(level)] = { + "raw_segments": len(raw_segments), + "triangle_fallback_segments": triangle_fallback_segments, + "line_parts": 0, + "length_m": 0.0, + } + continue + geom = union_all(parts) + contours[level] = geom + final_parts = v1.line_parts(geom) + stats[str(level)] = { + "raw_segments": len(raw_segments), + "triangle_fallback_segments": triangle_fallback_segments, + "line_parts": len(final_parts), + "length_m": round( + sum(part.length for part in final_parts), + 3, + ), + } + + return contours, stats + + +def main() -> int: + args = parse_args() + config_path = args.config.resolve() + cfg = v1.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']}", flush=True) + print(f"BBOX {cfg['bbox_wgs84']}", flush=True) + print(f"OUTPUT {output_dir}", flush=True) + + conn = v1.connect_db() + try: + points, sounding_stats = v1.fetch_soundings(conn, cfg) + print( + "SOUNDINGS", + json.dumps(sounding_stats, ensure_ascii=False), + flush=True, + ) + water_mask, water_rows = v1.fetch_water_mask(conn, cfg) + print(f"WATER_AREAS {water_rows}", flush=True) + finally: + conn.close() + xs, ys, grid, surface_stats = v2.build_idw_grid( + points, + water_mask, + cfg, + ) + water_grid = make_water_grid_mask(xs, ys, water_mask) + repaired_grid, gap_stats = repair_small_water_gaps( + grid, + water_grid, + cfg, + ) + + print( + "SURFACE", + json.dumps(surface_stats, ensure_ascii=False), + flush=True, + ) + print( + "GAP_REPAIR", + json.dumps(gap_stats, ensure_ascii=False), + flush=True, + ) + + contours, contour_stats = extract_dense_contours( + xs, + ys, + repaired_grid, + water_mask, + cfg, + ) + if not contours: + raise RuntimeError("no dense contours generated") + + print( + "CONTOURS", + json.dumps(contour_stats, ensure_ascii=False), + flush=True, + ) + + debug_feature_count = v1.write_debug_geojson( + contours, + cfg, + output_dir, + ) + v1.write_tilejson(cfg, output_dir) + tile_stats = v1.build_tiles(contours, cfg, output_dir) + print( + "TILES", + json.dumps(tile_stats, ensure_ascii=False), + flush=True, + ) + + finished = datetime.now(timezone.utc) + manifest = { + "schema": "navsea.cm95.bathymetry-build-result/v3", + "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(value) for value in cfg["contour_levels_m"] + ], + "major_contours_m": [ + float(value) for value in cfg["major_contours_m"] + ], + "native_scale_max": int(cfg["native_scale_max"]), + "water_object_classes": [ + str(value) for value in cfg["water_object_classes"] + ], + "surface_method": "local_idw_regular_grid_with_small_gap_repair", + "surface_parameters": { + key: cfg[key] + for key in ( + "grid_spacing_m", + "idw_search_radius_m", + "idw_max_nearest_m", + "idw_min_neighbors", + "idw_max_neighbors", + "idw_power", + "gap_fill_max_component_nodes", + "gap_fill_radius_cells", + "gap_fill_min_neighbors", + "chaikin_iterations", + ) + }, + "sounding_stats": sounding_stats, + "water_area_rows": water_rows, + "surface_stats": surface_stats, + "gap_repair_stats": gap_stats, + "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", + flush=True, + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/cm95_bathymetry_build_v4_tin_dense.py b/tools/cm95_bathymetry_build_v4_tin_dense.py new file mode 100644 index 0000000..a5764a0 --- /dev/null +++ b/tools/cm95_bathymetry_build_v4_tin_dense.py @@ -0,0 +1,399 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import json +import math +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import numpy as np +from shapely.geometry import MultiPoint + +import cm95_bathymetry_build_v1 as v1 +import cm95_bathymetry_build_v2_smooth as v2 +import cm95_bathymetry_build_v3_dense as v3 + +ROOT = Path(__file__).resolve().parents[1] +DEFAULT_CONFIG = ( + ROOT + / "config" + / "bathymetry" + / "cm95_bathymetry_north_kyushu_v4_tin_dense.yaml" +) +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Build dense bathymetry contours from a TIN-constrained continuous depth surface." + ) + parser.add_argument("--config", type=Path, default=DEFAULT_CONFIG) + parser.add_argument("--output-dir", type=Path, default=None) + return parser.parse_args() + + +def _edge_length( + a: tuple[float, float], + b: tuple[float, float], +) -> float: + return math.hypot(a[0] - b[0], a[1] - b[1]) + + +def _triangle_grid_bounds( + coords: list[tuple[float, float]], + xs: np.ndarray, + ys: np.ndarray, +) -> tuple[int, int, int, int] | None: + minx = min(point[0] for point in coords) + maxx = max(point[0] for point in coords) + miny = min(point[1] for point in coords) + maxy = max(point[1] for point in coords) + + ix0 = max(0, int(np.searchsorted(xs, minx, side="left"))) + ix1 = min(len(xs) - 1, int(np.searchsorted(xs, maxx, side="right") - 1)) + iy0 = max(0, int(np.searchsorted(ys, miny, side="left"))) + iy1 = min(len(ys) - 1, int(np.searchsorted(ys, maxy, side="right") - 1)) + + if ix0 > ix1 or iy0 > iy1: + return None + return ix0, ix1, iy0, iy1 +def _rasterize_triangle( + grid: np.ndarray, + water_grid: np.ndarray, + xs: np.ndarray, + ys: np.ndarray, + coords: list[tuple[float, float]], + depths: list[float], +) -> int: + bounds = _triangle_grid_bounds(coords, xs, ys) + if bounds is None: + return 0 + ix0, ix1, iy0, iy1 = bounds + + x0, y0 = coords[0] + x1, y1 = coords[1] + x2, y2 = coords[2] + d0, d1, d2 = depths + + denom = ( + (y1 - y2) * (x0 - x2) + + (x2 - x1) * (y0 - y2) + ) + if abs(denom) < 1e-12: + return 0 + + sub_xs = xs[ix0 : ix1 + 1] + sub_ys = ys[iy0 : iy1 + 1] + xx, yy = np.meshgrid(sub_xs, sub_ys) + + w0 = ( + (y1 - y2) * (xx - x2) + + (x2 - x1) * (yy - y2) + ) / denom + w1 = ( + (y2 - y0) * (xx - x2) + + (x0 - x2) * (yy - y2) + ) / denom + w2 = 1.0 - w0 - w1 + + eps = 1e-9 + inside = ( + (w0 >= -eps) + & (w1 >= -eps) + & (w2 >= -eps) + ) + water = water_grid[iy0 : iy1 + 1, ix0 : ix1 + 1] + target = grid[iy0 : iy1 + 1, ix0 : ix1 + 1] + writable = inside & water & ~np.isfinite(target) + + if not np.any(writable): + return 0 + + values = w0 * d0 + w1 * d1 + w2 * d2 + target[writable] = values[writable] + return int(np.count_nonzero(writable)) + + +def build_tin_surface( + points: dict[tuple[float, float], float], + water_mask, + cfg: dict[str, Any], +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, dict[str, Any]]: + west, south, east, north = [float(value) for value in cfg["bbox_wgs84"]] + minx, miny = v1.lonlat_to_mercator(west, south) + maxx, maxy = v1.lonlat_to_mercator(east, north) + + spacing = float(cfg["grid_spacing_m"]) + max_edge = float(cfg["max_triangle_edge_m"]) + + xs = v2._grid_axis(minx, maxx, spacing) + ys = v2._grid_axis(miny, maxy, spacing) + water_grid = v3.make_water_grid_mask(xs, ys, water_mask) + grid = np.full((len(ys), len(xs)), np.nan, dtype=np.float64) + + depth_by_coord = dict(points) + triangles = v1.delaunay_triangles(MultiPoint(list(points.keys()))) + + triangle_total = 0 + triangle_accepted = 0 + triangle_rejected_long_edge = 0 + triangle_rejected_lookup = 0 + raster_writes = 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 = [ + float(depth_by_coord[v1.coord_key(*coord)]) + for coord in coords + ] + except KeyError: + triangle_rejected_lookup += 1 + 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: + triangle_rejected_long_edge += 1 + continue + + triangle_accepted += 1 + raster_writes += _rasterize_triangle( + grid, + water_grid, + xs, + ys, + coords, + depths, + ) + + finite = np.isfinite(grid) + water_nodes = int(np.count_nonzero(water_grid)) + covered_water_nodes = int(np.count_nonzero(finite & water_grid)) + uncovered_water_nodes = int(np.count_nonzero(water_grid & ~finite)) + + stats = { + "grid_width": int(len(xs)), + "grid_height": int(len(ys)), + "grid_nodes": int(grid.size), + "water_nodes": water_nodes, + "covered_water_nodes": covered_water_nodes, + "uncovered_water_nodes": uncovered_water_nodes, + "coverage_ratio_of_water": round( + covered_water_nodes / float(water_nodes), + 6, + ) if water_nodes else 0.0, + "triangle_total": triangle_total, + "triangle_accepted": triangle_accepted, + "triangle_rejected_long_edge": triangle_rejected_long_edge, + "triangle_rejected_lookup": triangle_rejected_lookup, + "raster_writes": raster_writes, + } + return xs, ys, grid, water_grid, stats + + +def smooth_surface( + grid: np.ndarray, + water_grid: np.ndarray, + cfg: dict[str, Any], +) -> tuple[np.ndarray, dict[str, Any]]: + iterations = int(cfg.get("surface_smooth_iterations", 0)) + strength = float(cfg.get("surface_smooth_strength", 0.0)) + if iterations <= 0 or strength <= 0.0: + return grid.copy(), { + "iterations": 0, + "strength": 0.0, + "smoothed_nodes": 0, + } + + current = grid.copy() + valid_base = water_grid & np.isfinite(current) + smoothed_nodes = 0 + + for _ in range(iterations): + values_sum = np.zeros(current.shape, dtype=np.float64) + values_count = np.zeros(current.shape, dtype=np.int16) + + for dy, dx in v3._NEIGHBORS8: + src_y0 = max(0, -dy) + src_y1 = current.shape[0] - max(0, dy) + src_x0 = max(0, -dx) + src_x1 = current.shape[1] - max(0, dx) + + dst_y0 = max(0, dy) + dst_y1 = current.shape[0] - max(0, -dy) + dst_x0 = max(0, dx) + dst_x1 = current.shape[1] - max(0, -dx) + + source = current[src_y0:src_y1, src_x0:src_x1] + source_valid = np.isfinite(source) + target_sum = values_sum[dst_y0:dst_y1, dst_x0:dst_x1] + target_count = values_count[dst_y0:dst_y1, dst_x0:dst_x1] + target_sum[source_valid] += source[source_valid] + target_count[source_valid] += 1 + + can_smooth = valid_base & (values_count >= 3) + neighbor_mean = np.zeros(current.shape, dtype=np.float64) + neighbor_mean[can_smooth] = ( + values_sum[can_smooth] + / values_count[can_smooth] + ) + + updated = current.copy() + updated[can_smooth] = ( + current[can_smooth] * (1.0 - strength) + + neighbor_mean[can_smooth] * strength + ) + current = updated + smoothed_nodes += int(np.count_nonzero(can_smooth)) + + return current, { + "iterations": iterations, + "strength": strength, + "smoothed_nodes": smoothed_nodes, + } + + +def main() -> int: + args = parse_args() + config_path = args.config.resolve() + cfg = v1.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']}", flush=True) + print(f"BBOX {cfg['bbox_wgs84']}", flush=True) + print(f"OUTPUT {output_dir}", flush=True) + conn = v1.connect_db() + try: + points, sounding_stats = v1.fetch_soundings(conn, cfg) + print( + "SOUNDINGS", + json.dumps(sounding_stats, ensure_ascii=False), + flush=True, + ) + water_mask, water_rows = v1.fetch_water_mask(conn, cfg) + print(f"WATER_AREAS {water_rows}", flush=True) + finally: + conn.close() + + xs, ys, raw_grid, water_grid, surface_stats = build_tin_surface( + points, + water_mask, + cfg, + ) + print( + "TIN_SURFACE", + json.dumps(surface_stats, ensure_ascii=False), + flush=True, + ) + + smooth_grid, smooth_stats = smooth_surface( + raw_grid, + water_grid, + cfg, + ) + print( + "SURFACE_SMOOTH", + json.dumps(smooth_stats, ensure_ascii=False), + flush=True, + ) + + contours, contour_stats = v3.extract_dense_contours( + xs, + ys, + smooth_grid, + water_mask, + cfg, + ) + if not contours: + raise RuntimeError("no v4 contours generated") + print( + "CONTOURS", + json.dumps(contour_stats, ensure_ascii=False), + flush=True, + ) + debug_feature_count = v1.write_debug_geojson( + contours, + cfg, + output_dir, + ) + v1.write_tilejson(cfg, output_dir) + tile_stats = v1.build_tiles(contours, cfg, output_dir) + print( + "TILES", + json.dumps(tile_stats, ensure_ascii=False), + flush=True, + ) + + finished = datetime.now(timezone.utc) + manifest = { + "schema": "navsea.cm95.bathymetry-build-result/v4", + "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(value) for value in cfg["bbox_wgs84"]], + "source_layer": str(cfg["source_layer"]), + "source_kind": str(cfg["source_kind"]), + "reliability": str(cfg["reliability"]), + "contour_levels_m": [ + float(value) for value in cfg["contour_levels_m"] + ], + "major_contours_m": [ + float(value) for value in cfg["major_contours_m"] + ], + "native_scale_max": int(cfg["native_scale_max"]), + "water_object_classes": [ + str(value) for value in cfg["water_object_classes"] + ], + "surface_method": "tin_barycentric_grid_with_light_smoothing", + "surface_parameters": { + key: cfg[key] + for key in ( + "grid_spacing_m", + "max_triangle_edge_m", + "surface_smooth_iterations", + "surface_smooth_strength", + "chaikin_iterations", + ) + }, + "sounding_stats": sounding_stats, + "water_area_rows": water_rows, + "surface_stats": surface_stats, + "surface_smooth_stats": smooth_stats, + "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", + flush=True, + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/cm95_bathymetry_build_v5_adaptive.py b/tools/cm95_bathymetry_build_v5_adaptive.py new file mode 100644 index 0000000..9abd862 --- /dev/null +++ b/tools/cm95_bathymetry_build_v5_adaptive.py @@ -0,0 +1,588 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import bisect +import json +import math +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import numpy as np +from shapely import line_merge, union_all +from shapely.geometry import LineString, MultiPoint + +import cm95_bathymetry_build_v1 as v1 +import cm95_bathymetry_build_v2_smooth as v2 +import cm95_bathymetry_build_v3_dense as v3 +import cm95_bathymetry_build_v4_tin_dense as v4 + +ROOT = Path(__file__).resolve().parents[1] +DEFAULT_CONFIG = ( + ROOT + / "config" + / "bathymetry" + / "cm95_bathymetry_north_kyushu_v5_adaptive.yaml" +) +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="Build adaptive-TIN dense bathymetry contours." + ) + parser.add_argument("--config", type=Path, default=DEFAULT_CONFIG) + parser.add_argument("--output-dir", type=Path, default=None) + return parser.parse_args() + + +def _percentiles(values: list[float]) -> dict[str, float]: + if not values: + return {} + ordered = sorted(float(value) for value in values) + result: dict[str, float] = {} + for label, fraction in ( + ("p10", 0.10), + ("p25", 0.25), + ("p50", 0.50), + ("p75", 0.75), + ("p90", 0.90), + ("p95", 0.95), + ): + index = min( + len(ordered) - 1, + int(round((len(ordered) - 1) * fraction)), + ) + result[label] = round(ordered[index], 3) + return result +def prepare_adaptive_triangles( + points: dict[tuple[float, float], float], + cfg: dict[str, Any], +) -> tuple[list[dict[str, Any]], dict[tuple[float, float], float], dict[str, Any]]: + triangle_geoms = v1.delaunay_triangles(MultiPoint(list(points.keys()))) + records: list[dict[str, Any]] = [] + local_spacing: dict[tuple[float, float], float] = {} + + lookup_failures = 0 + for triangle in triangle_geoms.geoms: + coords = [ + (float(x), float(y)) + for x, y, *rest in list(triangle.exterior.coords)[:3] + ] + try: + depths = [ + float(points[v1.coord_key(*coord)]) + for coord in coords + ] + except KeyError: + lookup_failures += 1 + continue + + edges = ( + v4._edge_length(coords[0], coords[1]), + v4._edge_length(coords[1], coords[2]), + v4._edge_length(coords[2], coords[0]), + ) + records.append( + { + "coords": coords, + "depths": depths, + "edges": edges, + "longest": max(edges), + } + ) + + for index, coord in enumerate(coords): + key = v1.coord_key(*coord) + if index == 0: + incident = min(edges[0], edges[2]) + elif index == 1: + incident = min(edges[0], edges[1]) + else: + incident = min(edges[1], edges[2]) + previous = local_spacing.get(key) + if previous is None or incident < previous: + local_spacing[key] = incident + + factor = float(cfg["adaptive_edge_factor"]) + min_limit = float(cfg["adaptive_edge_min_m"]) + max_limit = float(cfg["adaptive_edge_max_m"]) + limits: list[float] = [] + primary = 0 + fallback = 0 + rejected_global = 0 + + for record in records: + spacings = sorted( + local_spacing[v1.coord_key(*coord)] + for coord in record["coords"] + ) + median_spacing = float(spacings[1]) + adaptive_limit = max( + min_limit, + min(max_limit, median_spacing * factor), + ) + record["local_spacing"] = median_spacing + record["adaptive_limit"] = adaptive_limit + + if record["longest"] <= adaptive_limit: + triangle_class = "primary" + primary += 1 + elif record["longest"] <= max_limit: + triangle_class = "fallback" + fallback += 1 + else: + triangle_class = "rejected_global" + rejected_global += 1 + + record["triangle_class"] = triangle_class + limits.append(adaptive_limit) + + spacing_values = list(local_spacing.values()) + stats = { + "triangle_total": int(len(triangle_geoms.geoms)), + "triangle_records": len(records), + "triangle_lookup_failures": lookup_failures, + "triangle_primary": primary, + "triangle_fallback": fallback, + "triangle_rejected_global": rejected_global, + "point_spacing_percentiles_m": _percentiles(spacing_values), + "adaptive_limit_percentiles_m": _percentiles(limits), + "adaptive_edge_factor": factor, + "adaptive_edge_min_m": min_limit, + "adaptive_edge_max_m": max_limit, + } + return records, local_spacing, stats +def build_adaptive_tin_surface( + points: dict[tuple[float, float], float], + water_mask, + cfg: dict[str, Any], +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, dict[str, Any]]: + west, south, east, north = [float(value) for value in cfg["bbox_wgs84"]] + minx, miny = v1.lonlat_to_mercator(west, south) + maxx, maxy = v1.lonlat_to_mercator(east, north) + + spacing = float(cfg["grid_spacing_m"]) + xs = v2._grid_axis(minx, maxx, spacing) + ys = v2._grid_axis(miny, maxy, spacing) + water_grid = v3.make_water_grid_mask(xs, ys, water_mask) + + grid = np.full((len(ys), len(xs)), np.nan, dtype=np.float64) + local_spacing_grid = np.full( + (len(ys), len(xs)), + np.nan, + dtype=np.float64, + ) + + records, _point_spacing, triangle_stats = prepare_adaptive_triangles( + points, + cfg, + ) + + primary_writes = 0 + fallback_writes = 0 + fallback_spacing = float(cfg["smooth_spacing_high_m"]) + + for triangle_class in ("primary", "fallback"): + for record in records: + if record["triangle_class"] != triangle_class: + continue + + before = np.isfinite(grid) + written_count = v4._rasterize_triangle( + grid, + water_grid, + xs, + ys, + record["coords"], + record["depths"], + ) + written = np.isfinite(grid) & ~before + + if triangle_class == "primary": + primary_writes += written_count + if np.any(written): + local_spacing_grid[written] = float( + record["local_spacing"] + ) + else: + fallback_writes += written_count + if np.any(written): + local_spacing_grid[written] = fallback_spacing + + finite = np.isfinite(grid) + water_nodes = int(np.count_nonzero(water_grid)) + covered = int(np.count_nonzero(finite & water_grid)) + uncovered = int(np.count_nonzero(water_grid & ~finite)) + + stats = { + **triangle_stats, + "grid_width": int(len(xs)), + "grid_height": int(len(ys)), + "grid_nodes": int(grid.size), + "water_nodes": water_nodes, + "covered_water_nodes": covered, + "uncovered_water_nodes": uncovered, + "coverage_ratio_of_water": round( + covered / float(water_nodes), + 6, + ) if water_nodes else 0.0, + "primary_raster_writes": primary_writes, + "fallback_raster_writes": fallback_writes, + "raster_writes": primary_writes + fallback_writes, + } + + return xs, ys, grid, water_grid, local_spacing_grid, stats + + +def adaptive_smooth_surface( + grid: np.ndarray, + water_grid: np.ndarray, + local_spacing_grid: np.ndarray, + cfg: dict[str, Any], +) -> tuple[np.ndarray, dict[str, Any]]: + iterations = int(cfg.get("surface_smooth_iterations", 1)) + spacing_low = float(cfg["smooth_spacing_low_m"]) + spacing_high = float(cfg["smooth_spacing_high_m"]) + strength_low = float(cfg["smooth_strength_low"]) + strength_high = float(cfg["smooth_strength_high"]) + + current = grid.copy() + valid_base = water_grid & np.isfinite(current) + total_smoothed = 0 + + spacing_for_strength = local_spacing_grid.copy() + spacing_for_strength[ + valid_base & ~np.isfinite(spacing_for_strength) + ] = spacing_high + denominator = max(spacing_high - spacing_low, 1.0) + t = np.clip( + (spacing_for_strength - spacing_low) / denominator, + 0.0, + 1.0, + ) + strength_grid = strength_low + t * ( + strength_high - strength_low + ) + + for _ in range(iterations): + values_sum = np.zeros(current.shape, dtype=np.float64) + values_count = np.zeros(current.shape, dtype=np.int16) + + for dy, dx in v3._NEIGHBORS8: + src_y0 = max(0, -dy) + src_y1 = current.shape[0] - max(0, dy) + src_x0 = max(0, -dx) + src_x1 = current.shape[1] - max(0, dx) + + dst_y0 = max(0, dy) + dst_y1 = current.shape[0] - max(0, -dy) + dst_x0 = max(0, dx) + dst_x1 = current.shape[1] - max(0, -dx) + + source = current[src_y0:src_y1, src_x0:src_x1] + source_valid = np.isfinite(source) + target_sum = values_sum[ + dst_y0:dst_y1, + dst_x0:dst_x1, + ] + target_count = values_count[ + dst_y0:dst_y1, + dst_x0:dst_x1, + ] + target_sum[source_valid] += source[source_valid] + target_count[source_valid] += 1 + can_smooth = valid_base & (values_count >= 3) + neighbor_mean = np.zeros(current.shape, dtype=np.float64) + neighbor_mean[can_smooth] = ( + values_sum[can_smooth] + / values_count[can_smooth] + ) + + updated = current.copy() + local_strength = strength_grid[can_smooth] + updated[can_smooth] = ( + current[can_smooth] * (1.0 - local_strength) + + neighbor_mean[can_smooth] * local_strength + ) + current = updated + total_smoothed += int(np.count_nonzero(can_smooth)) + + valid_strengths = strength_grid[ + valid_base & np.isfinite(strength_grid) + ] + return current, { + "iterations": iterations, + "spacing_low_m": spacing_low, + "spacing_high_m": spacing_high, + "strength_low": strength_low, + "strength_high": strength_high, + "smoothed_nodes": total_smoothed, + "strength_percentiles": _percentiles( + [float(value) for value in valid_strengths] + ), + } +def extract_contours_single_pass( + xs: np.ndarray, + ys: np.ndarray, + grid: np.ndarray, + water_mask, + cfg: dict[str, Any], +) -> tuple[dict[float, Any], dict[str, Any]]: + levels = sorted(float(value) for value in cfg["contour_levels_m"]) + segments: dict[float, list[LineString]] = { + level: [] for level in levels + } + fallback_counts: dict[float, int] = { + level: 0 for level in levels + } + + for iy in range(len(ys) - 1): + y0 = float(ys[iy]) + y1 = float(ys[iy + 1]) + for ix in range(len(xs) - 1): + values = ( + float(grid[iy, ix]), + float(grid[iy, ix + 1]), + float(grid[iy + 1, ix + 1]), + float(grid[iy + 1, ix]), + ) + finite_values = [ + value for value in values if math.isfinite(value) + ] + if len(finite_values) < 3: + continue + lo = min(finite_values) + hi = max(finite_values) + if lo == hi: + continue + + start = bisect.bisect_left(levels, lo) + stop = bisect.bisect_right(levels, hi) + if start >= stop: + continue + + finite_count = len(finite_values) + for level in levels[start:stop]: + cell_segments = v3._cell_contours( + float(xs[ix]), + float(xs[ix + 1]), + y0, + y1, + values, + level, + ) + segments[level].extend(cell_segments) + if finite_count == 3: + fallback_counts[level] += len(cell_segments) + contours: dict[float, Any] = {} + stats: dict[str, Any] = {} + min_segment = float(cfg["min_segment_length_m"]) + simplify_m = float(cfg.get("regional_simplify_m", 0)) + smooth_iterations = int(cfg.get("chaikin_iterations", 0)) + + for level in levels: + raw = segments[level] + if not raw: + stats[str(level)] = { + "raw_segments": 0, + "triangle_fallback_segments": 0, + "line_parts": 0, + "length_m": 0.0, + } + continue + + merged = line_merge(union_all(raw)) + parts: list[LineString] = [] + for part in v1.line_parts(merged): + if part.length < min_segment: + continue + smoothed = v2.chaikin_smooth( + part, + smooth_iterations, + ) + if simplify_m > 0: + smoothed = smoothed.simplify( + simplify_m, + preserve_topology=False, + ) + clipped = smoothed.intersection(water_mask) + parts.extend( + candidate + for candidate in v1.line_parts(clipped) + if candidate.length >= min_segment + ) + + if not parts: + stats[str(level)] = { + "raw_segments": len(raw), + "triangle_fallback_segments": fallback_counts[level], + "line_parts": 0, + "length_m": 0.0, + } + continue + geom = union_all(parts) + contours[level] = geom + final_parts = v1.line_parts(geom) + stats[str(level)] = { + "raw_segments": len(raw), + "triangle_fallback_segments": fallback_counts[level], + "line_parts": len(final_parts), + "length_m": round( + sum(part.length for part in final_parts), + 3, + ), + } + + return contours, stats + + +def main() -> int: + args = parse_args() + cfg = v1.load_config(args.config.resolve()) + 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']}", flush=True) + print(f"BBOX {cfg['bbox_wgs84']}", flush=True) + print(f"OUTPUT {output_dir}", flush=True) + + conn = v1.connect_db() + try: + points, sounding_stats = v1.fetch_soundings(conn, cfg) + print( + "SOUNDINGS", + json.dumps(sounding_stats, ensure_ascii=False), + flush=True, + ) + water_mask, water_rows = v1.fetch_water_mask(conn, cfg) + print(f"WATER_AREAS {water_rows}", flush=True) + finally: + conn.close() + ( + xs, + ys, + raw_grid, + water_grid, + local_spacing_grid, + surface_stats, + ) = build_adaptive_tin_surface( + points, + water_mask, + cfg, + ) + print( + "ADAPTIVE_TIN_SURFACE", + json.dumps(surface_stats, ensure_ascii=False), + flush=True, + ) + + repaired_grid, gap_stats = v3.repair_small_water_gaps( + raw_grid, + water_grid, + cfg, + ) + newly_filled = ( + np.isfinite(repaired_grid) + & ~np.isfinite(raw_grid) + ) + local_spacing_grid[newly_filled] = float( + cfg["smooth_spacing_high_m"] + ) + print( + "GAP_REPAIR", + json.dumps(gap_stats, ensure_ascii=False), + flush=True, + ) + + smooth_grid, smooth_stats = adaptive_smooth_surface( + repaired_grid, + water_grid, + local_spacing_grid, + cfg, + ) + print( + "ADAPTIVE_SMOOTH", + json.dumps(smooth_stats, ensure_ascii=False), + flush=True, + ) + + contours, contour_stats = extract_contours_single_pass( + xs, + ys, + smooth_grid, + water_mask, + cfg, + ) + if not contours: + raise RuntimeError("no v5 contours generated") + print( + "CONTOURS", + json.dumps(contour_stats, ensure_ascii=False), + flush=True, + ) + + debug_feature_count = v1.write_debug_geojson( + contours, + cfg, + output_dir, + ) + v1.write_tilejson(cfg, output_dir) + tile_stats = v1.build_tiles(contours, cfg, output_dir) + print( + "TILES", + json.dumps(tile_stats, ensure_ascii=False), + flush=True, + ) + + finished = datetime.now(timezone.utc) + manifest = { + "schema": "navsea.cm95.bathymetry-build-result/v5", + "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(value) for value in cfg["bbox_wgs84"]], + "source_layer": str(cfg["source_layer"]), + "source_kind": str(cfg["source_kind"]), + "reliability": str(cfg["reliability"]), + "contour_levels_m": [ + float(value) for value in cfg["contour_levels_m"] + ], + "major_contours_m": [ + float(value) for value in cfg["major_contours_m"] + ], + "surface_method": "adaptive_tin_barycentric_with_density_smoothing", + "sounding_stats": sounding_stats, + "water_area_rows": water_rows, + "surface_stats": surface_stats, + "gap_repair_stats": gap_stats, + "surface_smooth_stats": smooth_stats, + "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", + flush=True, + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tools/cm95_linux_full_rebuild_preflight.py b/tools/cm95_linux_full_rebuild_preflight.py new file mode 100755 index 0000000..35b0b90 --- /dev/null +++ b/tools/cm95_linux_full_rebuild_preflight.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +from __future__ import annotations + +from pathlib import Path + +EXPECTED_PROJECT_VERSION = "0.4.7" +EXPECTED_RUNTIME_VERSION = "0.4.7-production-audited-nav-name-conflict-omission" +EXPECTED_BUILDER = "cm95pbf-builder/0.4.6-spatial-candidate-query" + + +def fail(message: str) -> None: + print(f"FAIL: {message}") + raise SystemExit(2) + + +def require_fragment(path: Path, fragment: str) -> None: + text = path.read_text(encoding="utf-8") + if fragment not in text: + fail(f"{path} 缺少: {fragment}") + + +def main() -> None: + root = Path.cwd() + pyproject = root / "pyproject.toml" + init_py = root / "src/cm95pbf/__init__.py" + builder_py = root / "src/cm95pbf/builder.py" + release_py = root / "src/cm95pbf/release.py" + database_py = root / "src/cm93db/database.py" + + required_paths = [ + pyproject, + init_py, + builder_py, + release_py, + database_py, + root / "config/mappings/cm95_to_navsea_v1.yaml", + root / "config/mappings/cm95_gap_disposition_japan_dataset_1_v1.yaml", + root / "config/pbf/cm95_tile_policy_japan_full_detail_v1.yaml", + root / "config/style/style.navsea-delivery-full-semantic-unified-20260808-cm95-v1.json", + root / "config/sprite/semantic-sprite-20260808.json", + ] + for path in required_paths: + if not path.exists(): + fail(f"缺少文件: {path}") + + require_fragment(pyproject, f'version = "{EXPECTED_PROJECT_VERSION}"') + require_fragment(init_py, EXPECTED_RUNTIME_VERSION) + require_fragment(init_py, EXPECTED_BUILDER) + + for fragment in [ + "mysql_compat_mode_for_connection", + "axis-order=long-lat", + "geom_from_wkb_sql", + "geom_from_text_sql", + "as_wkb_sql", + ]: + require_fragment(database_py, fragment) + + for fragment in [ + '"POINT": ("chart_point_geometry", "spx_point_geom")', + '"LINE": ("chart_line_geometry", "spx_line_geom")', + '"AREA": ("chart_area_geometry", "spx_area_geom")', + "FORCE INDEX ({spatial_index})", + "STRAIGHT_JOIN chart_feature f ON f.id=g.feature_id", + ]: + require_fragment(builder_py, fragment) + + for fragment in [ + "AMBIGUOUS_DIRECT_NAMES_OMIT_LABEL_FAIL_CLOSED", + "_collect_nav_name_conflict_expectations", + "_observe_nav_name_conflict_feature", + "_finalize_nav_name_conflict_audit", + "ambiguous_conflict_was_labeled", + 'failures.append("nav_name_conflict_audit")', + '"nav_name_conflict_omissions"', + '"nav_name_conflict_unique_outputs"', + ]: + require_fragment(release_py, fragment) + + release_text = release_py.read_text(encoding="utf-8") + if 'failures.append("nav_name_conflicts")' in release_text: + fail("仍存在旧 nav_name_conflicts unconditional hard-fail") + + try: + import cm95pbf + except Exception as exc: + fail(f"无法 import cm95pbf: {exc}") + + if getattr(cm95pbf, "__version__", None) != EXPECTED_RUNTIME_VERSION: + fail(f"installed runtime={getattr(cm95pbf, '__version__', None)!r}") + if getattr(cm95pbf, "BUILDER_CONTRACT_VERSION", None) != EXPECTED_BUILDER: + fail(f"installed builder={getattr(cm95pbf, 'BUILDER_CONTRACT_VERSION', None)!r}") + + print("CM95 Linux full rebuild preflight") + print("=" * 72) + print(f"Project version : {EXPECTED_PROJECT_VERSION} PASS") + print(f"Runtime version : {EXPECTED_RUNTIME_VERSION} PASS") + print(f"Builder contract : {EXPECTED_BUILDER} PASS") + print("Japan mapping/disposition : PASS") + print("Japan tile policy : PASS") + print("Style / sprite contract inputs : PASS") + print("Ambiguous-name finalizer gate : PASS") + print("Spatial-first builder invariant : PASS") + print("MySQL 5.7/8.0 spatial compat : PASS") + print("OVERALL : PASS") + + +if __name__ == "__main__": + main() diff --git a/tools/cm95_v1_4_6_query_equivalence.py b/tools/cm95_v1_4_6_query_equivalence.py index 5f500f4..5f05368 100644 --- a/tools/cm95_v1_4_6_query_equivalence.py +++ b/tools/cm95_v1_4_6_query_equivalence.py @@ -11,6 +11,7 @@ import pymysql from dotenv import load_dotenv from cm93db.config import load_config +from cm93db.database import geom_from_text_sql, mysql_compat_mode_for_connection TABLES = { @@ -37,6 +38,11 @@ def connect() -> pymysql.connections.Connection: ) +def envelope_geom_sql(conn) -> str: + mode = mysql_compat_mode_for_connection(conn) + return geom_from_text_sql(mode, "%s", 4326) + + def envelope_wkt(bbox: tuple[float, float, float, float]) -> str: west, south, east, north = bbox return ( @@ -76,13 +82,14 @@ def fetch_new(conn, dataset_id: int, scale: int, geometry_type: str, bbox): west, south, east, north = bbox table, idx = TABLES[geometry_type] wkt = envelope_wkt(bbox) + geom_sql = envelope_geom_sql(conn) with conn.cursor() as cur: cur.execute( f""" SELECT f.id AS feature_id FROM {table} g FORCE INDEX ({idx}) STRAIGHT_JOIN chart_feature f ON f.id=g.feature_id - WHERE MBRIntersects(g.geom, ST_GeomFromText(%s, 4326)) + WHERE MBRIntersects(g.geom, {geom_sql}) AND f.dataset_id=%s AND f.native_scale=%s AND f.geometry_type=%s @@ -100,6 +107,7 @@ def explain_new(conn, dataset_id: int, scale: int, geometry_type: str, bbox): west, south, east, north = bbox table, idx = TABLES[geometry_type] wkt = envelope_wkt(bbox) + geom_sql = envelope_geom_sql(conn) with conn.cursor() as cur: cur.execute( f""" @@ -107,7 +115,7 @@ def explain_new(conn, dataset_id: int, scale: int, geometry_type: str, bbox): SELECT f.id AS feature_id FROM {table} g FORCE INDEX ({idx}) STRAIGHT_JOIN chart_feature f ON f.id=g.feature_id - WHERE MBRIntersects(g.geom, ST_GeomFromText(%s, 4326)) + WHERE MBRIntersects(g.geom, {geom_sql}) AND f.dataset_id=%s AND f.native_scale=%s AND f.geometry_type=%s @@ -123,15 +131,16 @@ def explain_new(conn, dataset_id: int, scale: int, geometry_type: str, bbox): def soundings_old(conn, dataset_id: int, scale: int, bbox): wkt = envelope_wkt(bbox) + geom_sql = envelope_geom_sql(conn) with conn.cursor() as cur: cur.execute( - """ + f""" SELECT s.id FROM chart_sounding_point s JOIN chart_feature f ON f.id=s.feature_id WHERE f.dataset_id=%s AND f.native_scale=%s - AND MBRContains(ST_GeomFromText(%s, 4326), s.geom) + AND MBRContains({geom_sql}, s.geom) """, (dataset_id, scale, wkt), ) @@ -140,13 +149,14 @@ def soundings_old(conn, dataset_id: int, scale: int, bbox): def soundings_new(conn, dataset_id: int, scale: int, bbox): wkt = envelope_wkt(bbox) + geom_sql = envelope_geom_sql(conn) with conn.cursor() as cur: cur.execute( - """ + f""" SELECT s.id 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) + WHERE MBRContains({geom_sql}, s.geom) AND f.dataset_id=%s AND f.native_scale=%s """, @@ -158,14 +168,15 @@ def soundings_new(conn, dataset_id: int, scale: int, bbox): def soundings_explain_new(conn, dataset_id: int, scale: int, bbox): wkt = envelope_wkt(bbox) + geom_sql = envelope_geom_sql(conn) with conn.cursor() as cur: cur.execute( - """ + f""" EXPLAIN SELECT s.id 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) + WHERE MBRContains({geom_sql}, s.geom) AND f.dataset_id=%s AND f.native_scale=%s """, diff --git a/web/bathymetry_test.html b/web/bathymetry_test.html new file mode 100644 index 0000000..eb31ed1 --- /dev/null +++ b/web/bathymetry_test.html @@ -0,0 +1,191 @@ + + + + + +CM95 派生等深线独立 PBF 测试 + + + + +
+ + + + + diff --git a/web/bathymetry_test_server.py b/web/bathymetry_test_server.py new file mode 100644 index 0000000..e86eae9 --- /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-v5-adaptive" + + +def parse_args(): + p = argparse.ArgumentParser(description="CM95 independent bathymetry PBF test server") + p.add_argument("--host", default="192.168.200.184") + 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()