From 2513dd56f8ae87c66b616533ca61a13329c93481 Mon Sep 17 00:00:00 2001 From: OpenAI Codex Date: Sat, 2 May 2026 14:32:06 +0800 Subject: [PATCH] =?UTF-8?q?=E6=B7=BB=E5=8A=A0=E5=85=A8=E5=9B=BD=E4=B8=89?= =?UTF-8?q?=E5=B1=82=E7=94=9F=E6=88=90=E4=B8=8E=E6=B8=AF=E5=90=8D=E6=9F=A5?= =?UTF-8?q?FPC=E5=B7=A5=E5=85=B7?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- NavSea_全国三层数据生成与查看说明.md | 290 +++ STEP_RECORD.md | 1501 +++++++++++++- .../build_fish_port_20m_full_mysql_resume.py | 1823 +++++++++++++++++ coastline/build_japan_coast_grid_mysql.py | 565 +++++ coastline/build_japan_hazard_50m_mysql.py | 516 +++++ coastline/export_navgrid_mysql_assets.py | 485 +++++ find_fpc_by_port_name.py | 282 +++ .../navsea-coastline-fukuoka-saga-200m.html | 459 +++++ 8 files changed, 5920 insertions(+), 1 deletion(-) create mode 100644 NavSea_全国三层数据生成与查看说明.md create mode 100644 coastline/build_fish_port_20m_full_mysql_resume.py create mode 100644 coastline/build_japan_coast_grid_mysql.py create mode 100644 coastline/build_japan_hazard_50m_mysql.py create mode 100644 coastline/export_navgrid_mysql_assets.py create mode 100644 find_fpc_by_port_name.py create mode 100644 src/pbf/navsea-coastline-fukuoka-saga-200m.html diff --git a/NavSea_全国三层数据生成与查看说明.md b/NavSea_全国三层数据生成与查看说明.md new file mode 100644 index 0000000..e595fbe --- /dev/null +++ b/NavSea_全国三层数据生成与查看说明.md @@ -0,0 +1,290 @@ +# NavSea 全国四脚本一页面数据生成与查看说明 + +本文说明全国版三层数据的重算、导出、整体密度图生成和页面查看方式。 + +当前统一数据库： + +- `navsea_japan_coast_grid` + +当前全国三层： + +- `coast_200m`：全国海岸 200x200 +- `fish_port_20m`：全国渔港 20x20 +- `hazard_50m`：全国海上障碍 50x50 + +当前整体密度图： + +- `density_overview`：三档密度整体格子图 + +## 1. 生成顺序 + +建议按下面顺序重算： + +1. 先生成 `200x200` +2. 再生成 `20x20` +3. 最后生成 `50x50` +4. 再导出全国静态 JSON 资产 +5. 最后查看整体密度图 + +这样可以先把底盘准备好，再做上层数据。 + +## 2. 200x200 全国海岸重算 + +脚本： + +- [`coastline/build_japan_coast_grid_mysql.py`](/root/sourceserver/pbf/coastline/build_japan_coast_grid_mysql.py) + +默认会写入： + +- 数据库：`navsea_japan_coast_grid` +- 图层：`coast_200m` +- 导出目录：`out/coastline/japan_coast_grid_mysql` + +执行命令： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python coastline/build_japan_coast_grid_mysql.py +``` + +## 3. 20x20 全国渔港重算 + +脚本： + +- [`coastline/build_fish_port_20m_full_mysql_resume.py`](/root/sourceserver/pbf/coastline/build_fish_port_20m_full_mysql_resume.py) + +默认会写入： + +- 数据库：`navsea_japan_coast_grid` +- 图层：`fish_port_20m` + +推荐执行命令： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py \ + --resume \ + --max-scan-cells 200000 +``` + +如果只想重算某个 `PRC`，可以显式指定，例如： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py \ + --prc 02 \ + --max-scan-cells 200000 +``` + +如果只想重算某个具体渔港，可以直接用 `FPC`，例如： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py \ + --fpc 1210100 \ + --max-scan-cells 200000 +``` + +说明： + +- `PRC` 需要按两位字符串理解，`2` 和 `02` 等价 +- `FPC` 是单个渔港的港码，脚本会先反查它所属的 `PRC`，再只重算这个港口 +- 单港口模式会在日志里同时输出： + - 命中的 `200x200` 粗格数量 + - 最终写入的 `20x20` 格子数量 +- 默认先按渔港 bbox 找同区域 `coast_200m` 粗格，再把粗格切成 `20x20` +- 默认不会再允许无边界整片扫描 +- 如果确实要回退旧行为，需要手动加脚本里提供的相应开关 + +## 4. 根据港名查 FPC + +脚本： + +- [`find_fpc_by_port_name.py`](/root/sourceserver/pbf/find_fpc_by_port_name.py) + +用途： + +- 已知渔港名字时，先查出可能的 `FPC` +- 这个脚本会扫描 `coastline/C09-06.zip` 里的渔港要素 +- 支持模糊匹配和精确匹配 + +推荐用法： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python find_fpc_by_port_name.py 無垢島 +``` + +如果你想只看精确命中： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python find_fpc_by_port_name.py 無垢島 --exact +``` + +如果你想要 JSON 结果： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python find_fpc_by_port_name.py 無垢島 --json +``` + +如果原始包里没有港名字段，建议改用港名对照表： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python find_fpc_by_port_name.py 無垢島 --catalog out/port_name_catalog.csv +``` + +对照表最低需要这些列里的任意几列： + +- `name` / `港名` / `名称` +- `FPC` +- 可选 `PRC` + +说明： + +- 这份脚本主要查 `C09-06.zip` 里的渔港名称字段 +- 默认会把 `NA2`、`NA4`、`FCF` 等候选字段一起拿来比对 +- 输出里会给出： + - `FPC` + - `PRC` + - 主名称 + - 辅助名称 + - 源记录 id +- 如果原始包里没有这条名字，脚本会提示你改用 `--catalog` +## 5. 50x50 全国障碍重算 + +脚本： + +- [`coastline/build_japan_hazard_50m_mysql.py`](/root/sourceserver/pbf/coastline/build_japan_hazard_50m_mysql.py) + +默认会写入： + +- 数据库：`navsea_japan_coast_grid` +- 图层：`hazard_50m` +- PBF 根目录：`/home/wwwroot/pbf-delivery-full-20260418-rebuild` + +执行命令： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python coastline/build_japan_hazard_50m_mysql.py +``` + +这个脚本当前会纳入： + +- `navigation_hazard_area` +- `fixed_fishing_gear_area` +- `anchor_caution_hazard_area` +- `navigation_hazard_point` +- `anchor_caution_hazard_point` +- `navigation_marks` +- `baseline_area` 中 `canonical_object_type=breakwater` + +## 6. 全国三层静态 JSON 导出 + +脚本： + +- [`coastline/export_navgrid_mysql_assets.py`](/root/sourceserver/pbf/coastline/export_navgrid_mysql_assets.py) + +默认导出到： + +- `src/pbf/coastline-mysql/japan_national/` + +执行命令： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python coastline/export_navgrid_mysql_assets.py \ + --db-name navsea_japan_coast_grid \ + --out-dir src/pbf/coastline-mysql/japan_national +``` + +导出结果包括： + +- `coast_200m_grid.geojson` +- `fish_port_20m_grid.geojson` +- `hazard_50m_grid.geojson` +- `density_overview_grid.geojson` +- `manifest.json` + +说明： + +- `density_overview_grid.geojson` 是按密度优先级合并后的整体格子图 +- 同一区域如果存在更高密度格子，就优先保留更高密度 +- 当前密度优先级： + - `20x20` = 高密度，淡绿色 + - `50x50` = 中密度，黄色 + - `200x200` = 低密度，淡红色 + +## 7. HTML 预览页面 + +当前全国预览页面： + +- `http://192.168.200.184/newpec/navsea-coastline-fukuoka-saga-200m.html` + +说明： + +- 页面标题已经改成全国口径 +- 页面内有四个按钮： + - `200x200` + - `20x20` + - `50x50` + - `整体密度` + +页面对应的本地文件： + +- [`src/pbf/navsea-coastline-fukuoka-saga-200m.html`](/root/sourceserver/pbf/src/pbf/navsea-coastline-fukuoka-saga-200m.html) + +页面加载的数据源： + +- `./coastline-mysql/japan_coast_200m/manifest.json` +- `./coastline-mysql/japan_coast_200m/coast_200m_grid.geojson` +- `./coastline-mysql/japan_national/manifest.json` +- `./coastline-mysql/japan_national/fish_port_20m_grid.geojson` +- `./coastline-mysql/japan_national/hazard_50m_grid.geojson` +- `./coastline-mysql/japan_national/density_overview_grid.geojson` + +## 8. 页面查看方式 + +如果只是看全国 200x200： + +1. 打开页面 +2. 默认就是 `200x200` +3. 页面会自动定位到全国底盘范围 + +如果想切换看渔港或障碍： + +1. 点击 `20x20` +2. 点击 `50x50` +3. 点击 `整体密度` + +如果想强制刷新： + +1. 点击 `重新加载` + +## 9. 线下同步到网页目录 + +如果你重新生成了静态资产，想让线上页面立即看到新内容，通常需要把生成结果同步到网页目录： + +```bash +cp -r src/pbf/coastline-mysql/japan_national /mnt/sda1/www/newpec/ +cp src/pbf/navsea-coastline-fukuoka-saga-200m.html /mnt/sda1/www/newpec/navsea-coastline-fukuoka-saga-200m.html +``` + +## 10. 一句话版 + +最常用的整套命令是： + +```bash +cd /root/sourceserver/pbf +./.venv/bin/python coastline/build_japan_coast_grid_mysql.py +./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --resume --max-scan-cells 200000 +./.venv/bin/python coastline/build_japan_hazard_50m_mysql.py +./.venv/bin/python coastline/export_navgrid_mysql_assets.py --db-name navsea_japan_coast_grid --out-dir src/pbf/coastline-mysql/japan_national +``` + +然后打开： + +- `http://192.168.200.184/newpec/navsea-coastline-fukuoka-saga-200m.html` diff --git a/STEP_RECORD.md b/STEP_RECORD.md index 9c64c76..096e8c0 100644 --- a/STEP_RECORD.md +++ b/STEP_RECORD.md @@ -1,6 +1,6 @@ # 项目步骤记录 -最后更新：`2026-04-18` +最后更新：`2026-05-02` 仓库：`/root/sourceserver/pbf` 远端：`ssh://git@nas:2222/tei/pbf.git` @@ -24,6 +24,1223 @@ ## 最近 3 天摘要 +### 2026-05-02 全国 50x50 重算脚本 + +- 新增全国 50x50 海上障碍重算脚本： + - `coastline/build_japan_hazard_50m_mysql.py` +- 默认行为： + - 从全国 PBF 根目录 `'/home/wwwroot/pbf-delivery-full-20260418-rebuild'` 扫描 `z12` 瓦片 + - 识别 `navigation_hazard_area`、`fixed_fishing_gear_area`、`anchor_caution_hazard_area`、`navigation_hazard_point`、`anchor_caution_hazard_point`、`navigation_marks` + - 追加把 `baseline_area` 里的 `canonical_object_type=breakwater` 视作 50m 障碍底盘 + - 过滤 `fish_reef` 类对象后，重算 `hazard_50m` + - 写回统一库 `navsea_japan_coast_grid` + - 导出全国静态资产到 `src/pbf/coastline-mysql/japan_national/` +- 已验证： + - `./.venv/bin/python -m py_compile coastline/build_japan_hazard_50m_mysql.py` + - `./.venv/bin/python coastline/build_japan_hazard_50m_mysql.py --help` + - 抽样全国瓦片可以解析到 hazard layer + +### 2026-05-02 全国三层生成与查看说明 + +- 新增中文说明文档： + - `NavSea_全国四脚本一页面数据生成与查看说明.md` +- 说明内容包括： + - 全国 `200x200`、`20x20`、`50x50` 的重算脚本 + - 整体密度图 `density_overview` 的导出与查看 + - 各脚本的完整执行命令 + - 全国静态 JSON 导出命令 + - 全国预览 HTML 的完整 URL + - 线下同步到网页目录的命令 +- 已同步最新全国预览页与整体密度资产到： + - `/mnt/sda1/www/newpec/navsea-coastline-fukuoka-saga-200m.html` + - `/mnt/sda1/www/newpec/coastline-mysql/japan_national/density_overview_grid.geojson` +- 已修正全国预览页对 bbox 的容错： + - 现在遇到旧版 mercator bbox 时会自动回退到预设中心点，不再抛 `Invalid LngLat latitude value` +- 同步修正了导出脚本与 200m 构建脚本的 bbox 口径： + - `coastline/export_navgrid_mysql_assets.py` + - `coastline/build_japan_coast_grid_mysql.py` +- 已把 20m 重算流程收紧为： + - 先按渔港 bbox 找同区域 `coast_200m` 粗格 + - 再把粗格按 20m 直切并做海岸线重叠判断 + - 默认粗筛边距收紧为 0m，避免再做不必要的大范围外扩 + +### 2026-05-02 全国海岸/渔港/障碍三按钮预览页 + +- 已将 `src/pbf/navsea-coastline-fukuoka-saga-200m.html` 改造成全国入口页： + - 页面标题改为全国 + - 默认加载全国 `200x200` + - 增加四个切换按钮： + - `200x200` + - `20x20` + - `50x50` + - `整体密度` +- 已重新导出全国静态资产到： + - `src/pbf/coastline-mysql/japan_national/` +- 三个按钮对应的数据源现在是： + - `coast_200m`：全国海岸 200x200 + - `fish_port_20m`：全国渔港 20x20 + - `hazard_50m`：全国海上障碍 50x50 + - `density_overview`：三档密度整体格子图 +- 已同步到线上页面： + - `http://192.168.200.184/newpec/navsea-coastline-fukuoka-saga-200m.html` +- 已纳入固定基线的一部分，后续全国相关视觉确认优先看这页 + +### 2026-05-02 PRC 编号补零修正 + +- 已修正 `coastline/build_fish_port_20m_full_mysql_resume.py` 的 `PRC` 选择逻辑： + - 现在会把 `2` 自动归一成 `02` + - `PRC` 读取、选择、断点恢复和源名删除都统一按两位字符串处理 +- 这次修正的直接原因： + - `C09-06.zip` 里的原始 `PRC` 是两位编码 + - 传 `--prc 2` 时之前会匹配不到任何渔港线要素 +- 已验证： + - `normalize_prc_token('2') -> '02'` + - `discover_prcs(..., {'2'})` 不再漏掉 `02` + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + +### 2026-05-02 全国 200m / 20m / 障碍层统一汇聚脚本 + +- 新增可直接运行的汇聚脚本： + - `navsea_national_navigation_fusion.py` +- 当前默认汇聚的三层是： + - `navsea_japan_coast_grid.navsea_grid_cell` 里的 `coast_200m` + - `navsea_fukuoka_saga_grid.navsea_grid_cell` 里的 `fish_port_20m` + - `navsea_fukuoka_saga_grid.navsea_grid_cell` 里的 `hazard_50m` +- 输出格式： + - 单个 SQLite：`out/navsea_national_navigation_fusion.sqlite` + - 清单：`out/navsea_national_navigation_fusion.manifest.json` +- 脚本特性： + - 记录来源注册、图层注册、每层汇总与合并后的统一格网表 + - 支持 `--recreate` + - 支持 `--batch-size` + - 支持 `--emit-template` +- 已验证： + - `./.venv/bin/python -m py_compile navsea_national_navigation_fusion.py` + - `./.venv/bin/python navsea_national_navigation_fusion.py --help` +- 说明： + - 这只是 SQLite 版的全国融合骨架 + - 当前 MySQL 侧已经合并为单库 `navsea_japan_coast_grid` + +### 2026-05-02 渔港 20m 支持按单港 FPC 重算 + +- 已为 `coastline/build_fish_port_20m_full_mysql_resume.py` 增加单港口入口： + - 新增 `--fpc` + - 脚本会先按 `FPC` 反查所属 `PRC` + - 只处理该港口对应的线要素 +- 单港口模式会在日志里同时给出： + - 命中的 `200x200` 粗格数量 + - 最终写入的 `20x20` 格子数量 +- 单港口模式下，每个批次会先按当前 cell_id 清掉旧记录，再写回，便于重算同一港口 +- 已验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `discover_target_prcs(..., '1112040')` 可以反查到对应 `PRC` + - `collect_fish_lines_for_prc(..., '01', '1112040')` 可以只筛出这个港口的线要素 + +### 2026-05-02 按港名查 FPC 小工具 + +- 新增脚本： + - `find_fpc_by_port_name.py` +- 用途： + - 输入港名，扫描 `coastline/C09-06.zip` 里的渔港要素 + - 通过 `NA2`、`NA4`、`FCF` 等字段找候选 `FPC` + - 如果原始包里没有港名字段，也可以改用外部港名对照表 `--catalog` + - 支持精确匹配和模糊匹配 +- 文档已补到： + - `NavSea_全国三层数据生成与查看说明.md` +- 已验证： + - `./.venv/bin/python -m py_compile find_fpc_by_port_name.py` + +### 2026-05-02 无垢岛 PRC=50 改为单库运行 + +- 已把 `coastline/build_fish_port_20m_full_mysql_resume.py` 从双库连接收回到单库模式： + - 现在只打开一个 MySQL 连接 + - 同一个库里既读 `coast_200m`，也写 `fish_port_20m` + - 当前默认库改为 `navsea_japan_coast_grid` +- 这次收口的原因： + - 用户明确要求 `PRC=50` 的 20m 生成不要再同时开两个数据库 + - 以单库方式更容易保持 `coast_200m` / `fish_port_20m` 的同库可追溯性 +- 已验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` + +### 2026-05-02 两个 MySQL 库合并为一个 + +- 已执行合并脚本： + - `navsea_merge_navgrid_mysql.py` +- 合并方向： + - 保留 `navsea_japan_coast_grid` + - 把 `navsea_fukuoka_saga_grid` 里的 `fish_port_20m`、`hazard_50m`、`navsea_grid_import_state`、`navsea_grid_import_progress` 搬入全国库 + - 合并完成后删除 `navsea_fukuoka_saga_grid` +- 合并后目标库现状： + - `coast_200m`: `336215` + - `fish_port_20m`: `16300` + - `hazard_50m`: `176454` + - `navsea_grid_import_state`: `1` + - `navsea_grid_import_progress`: `39` +- 结果： + - 当前只保留一个 MySQL 库：`navsea_japan_coast_grid` + - 后续 `PRC=50` 的 20m 以及相关导出都以这个库为准 +- 已验证： + - `./.venv/bin/python navsea_merge_navgrid_mysql.py --dry-run` + - `./.venv/bin/python navsea_merge_navgrid_mysql.py --drop-source-db` + +### 2026-05-02 无垢岛漁港海岸线回退接入修正 + +- 已修正 `coastline/build_fish_port_20m_full_mysql_resume.py` 的海岸线加载逻辑： + - 先尝试 `C23-06_{PRC}_GML.zip` + - 如果同号包不存在，则按当前港区 bbox 退回扫描 `coastline/C23-06_*_GML.zip` + - 这样 `PRC=50` 不会再卡死在缺失的同号文件上 +- 追加修正： + - 回退扫描时使用的是海岸线函数内部的墨卡托米制筛选 + - 调用点必须先把港区 bbox 从经纬度转成墨卡托，否则会出现 `(-68, -166, ...)` 这种错误 bbox +- 进一步修正： + - 粗筛底盘现在改成从全国库 `navsea_japan_coast_grid` 读取 + - `navsea_japan_coast_grid.navsea_grid_cell` 实际存的是墨卡托米制坐标，查询时必须按米制直查，不能再转回经纬度 +- 已实测 `PRC=50`： + - 退回后从 `C23-06_44_GML.zip` 命中无垢岛附近海岸线 + - 命中线数 `6` + - 这说明无垢岛漁港可以借助现有海岸线底盘继续做 20m +- 最新复核： + - 无垢岛 bbox 在全国粗格库里可命中 `13` 个 `coast_200m` 窗口 +- 已验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + +### 2026-05-02 全国 200x200 底盘已存在 + +- 仓库里已经有全国海岸 200x200 相关资产与预览页： + - `src/pbf/coastline-mysql/japan_coast_200m/` + - `src/pbf/navsea-coastline-fukuoka-saga-200m.html` +- 数据库里也确实存在全国级 `coast_200m`： + - `navsea_japan_coast_grid` + - 当前 `coast_200m` 记录数约 `336215` +- 说明： + - 全国 200x200 地图是有的，不只是九州样区 + - 但它仍取决于当前已加载的海岸线来源包集合，缺包的区域不会凭空补出来 + +### 2026-05-02 无垢岛漁港 200m 覆盖复核 + +- 已复核 `PRC=50` 对应的无垢岛漁港坐标范围： + - 原始边界线 bbox 约为 `131.975066922648..131.977924005868, 33.1589557725031..33.160819912381` +- 已在两个当前库里逐一查询该点是否落入 `coast_200m`： + - `navsea_fukuoka_saga_grid` 命中数 `0` + - `navsea_japan_coast_grid` 命中数 `0` +- 结论： + - 当前仓库里能直接查到的两套 `coast_200m` 都没有覆盖无垢岛漁港 + - `navsea_japan_coast_grid` 的 `coast_200m` 来源分布里也没有 `C23-06_50` + - 所以无垢岛漁港并不在当前已加载的 200x200 底盘里 + +### 2026-05-01 渔港 20m 超大范围扫描僵死止损 + +- 用户现场日志显示任务卡在： + - `build_cells_for_geometry mask_hit row_idx=13315868 candidate_mask=3` + - `build_cells_for_geometry tile_hit idx=162390 bbox=(16232000.000, 5368000.000, 16232360.000, 5370000.000)` +- 判断： + - 这不是普通慢查询，而是 `build_cells_for_geometry` 已进入超大范围网格扫描 + - 直接原因是某些 `part` 没有有效 `coast_200m` 小窗时，旧逻辑会回退扫描完整 `fish_part` + - 该回退路径会产生十万级 tile / 千万级 row 检查，足以把 Linux 拖到近似僵死 +- 已修正 `coastline/build_fish_port_20m_full_mysql_resume.py`： + - 默认关闭无边界回退：没有有效 `coast_200m` 交集的 `part` 作为异常中止，不再静默跳过 + - 如确需旧行为，必须显式加 `--allow-unbounded-fallback` + - 如人工确认允许漏算某个无底盘 `part`，必须显式加 `--skip-missing-coarse` + - `coast_200m` 命中后不再扫描整个粗窗 box，而是扫描 `fish_part ∩ coast_200m窗口` + - 新增 `--max-scan-cells`，默认单次 `build_cells_for_geometry` 估算扫描格数超过 `200000` 直接中止 + - 日志新增 `grid=列x行` 和 `estimated_cells`，方便在真正扫描前判断风险 +- 已验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` +- 当前注意： + - `coastline/` 目录在当前工作树里仍是未跟踪目录，这次脚本修改不会出现在普通 `git diff` 中 + - 后续如果要正式提交渔港 20m 代码，需要单独决定是否把 `coastline/` 相关源码纳入版本控制，避免把导出大文件一起扫入提交 + +### 2026-05-01 渔港 20m 当前完成度复核 + +- 复核结果： + - `coastline/C09-06.zip` 里一共能识别出 `40` 个有效 PRC + - 当前数据库 `navsea_fukuoka_saga_grid` 里只有 `coast_200m` 和 `hazard_50m` 两层 + - 当前库里**没有** `fish_port_20m` 层数据 + - 当前磁盘上的 `out/fish_port_20m_full_resume.current_prc.json` 仍停在 `current_prc=01` +- 结论： + - 你贴出来的 `C23-06_01 ... C23-06_47` 统计是 **coast_200m** 的数据，不是渔港 20m + - 这份 `fish_port_20m` 任务在当前库里**还没有完成落库** + - 如果要判断 20x20 是否完成，应该查 `layer_name='fish_port_20m'`，当前结果为空 +- 额外更正： + - 用户截图里当前选中的库是 `navsea_japan_coast_grid` + - 这份库里只有 `coast_200m`，总数 `336215`，按 `source_name` 分布在 39 个县包 + - 这不是渔港 20m 任务库，不能拿来判断 `fish_port_20m` 是否完成 + +### 2026-05-01 渔港 20m 改为 part 内流式写库 + +- 已把 `coastline/build_fish_port_20m_full_mysql_resume.py` 的写库方式再收紧一层： + - 不再把一个 `part` 的所有 `rows` 整块攒起来 + - 改为在每个 `window` 生成后就直接进入批量写库 + - `commit_every` 仍然保留，负责控制每批提交多少条 +- 这样可以明显降低大 `part` 的峰值内存占用 +- 当前仍保留： + - `window_rows` 这一小段局部结果 + - `build_cells_for_geometry` 的返回值 + - 但已经不再把整个 `part` 的结果常驻在一个大列表里 + +### 2026-05-01 全国 full-audit 页 z12 右侧空白回归定位与临时修复 + +- 用户在 `http://192.168.200.184/newpec/navsea-compare-full-audit.html` 发现： + - 切到语义版 sprite / style 后，右侧 delivery 在 `zoom=12` 出现大块空白，表面像 PBF 损坏 +- 已按 2026-04-18 历史记录复核，确认这是同一类旧问题： + - 不是 sprite atlas 本身导致 + - 当前默认页仍指向旧的 `/home/wwwroot/pbf-delivery-full-semantic-20260416` + - 该目录部分 z12 高 extent tile 存在几何错移 +- 代表 tile 复核： + - `12/3527/1641.pbf` + - 原始 newpec / 20260418 rebuild：`y` 范围约 `-20480..1069056` + - 旧 semantic 20260416：`y` 范围约 `1024000..2113536` + - 这正是 2026-04-18 记录里的高 extent 重编码错移签名 +- 已把全国 full-audit 页默认 PBF 切到 2026-04-18 修复版： + - `/home/wwwroot/pbf-delivery-full-20260418-rebuild` + - URL：`http://192.168.200.184/pbf-delivery-full-20260418-rebuild/{z}/{x}/{y}.pbf` +- 已补充可供客户端直接远程加载的 extentfix style： + - 语义版：`http://192.168.200.184/newpec/domain/style.navsea-delivery-full-semantic-extentfix.json` + - 普通版：`http://192.168.200.184/newpec/domain/style.navsea-delivery-full-extentfix.json` + - 两者的 `sources.navsea_delivery.tiles` 都已指向 `/pbf-delivery-full-20260418-rebuild` +- 已按客户端一致性要求再次调整 full-audit 审计页： + - 默认右侧不再由 HTML 覆盖 PBF URL + - 默认右侧直接加载 `style.navsea-delivery-full-semantic-extentfix.json` 或 `style.navsea-delivery-full-extentfix.json` + - 也就是说：审计页看到的 style，就是客户端应加载的同一份远程 style + - 只有显式传入 `deliveryTileUrl=` 时，才作为临时实验入口覆盖 tile URL +- 后续全国 full/semantic 审计页固定按这个口径： + - 先改远程 style + - 审计页加载同一份远程 style + - 不用 HTML 默认覆盖 PBF URL 来伪造审计结果 +- 已继续压缩 full-audit 页面 UI： + - 顶部工具栏桌面端压到单行，长版本 chip 省略显示，减少地图空间占用 + - `点击对比` 面板支持展开 / 收起 + - `点击对比` 面板可拖拽到屏幕任意位置，并用 `localStorage` 保留位置与折叠状态 + - 默认落位从右下改为底部居中，避免一开页就压住右侧地图视野 + - 进一步修正了页面横向撑宽问题，保证左右 pane 真正等宽 +- 已更新并部署： + - 源文件：`src/pbf/navsea-compare-full-audit.html` + - 线上文件：`/mnt/sda1/www/newpec/navsea-compare-full-audit.html` + - 线上 style：`/mnt/sda1/www/newpec/domain/style.navsea-delivery-full-semantic-extentfix.json` + - 线上 style：`/mnt/sda1/www/newpec/domain/style.navsea-delivery-full-extentfix.json` + - 页面版本：`full-audit-r6-20260501-1745` +- 当前页面语义版口径： + - 使用 `style.navsea-delivery-full-semantic-extentfix.json` + - sprite / PBF 都以该 style 内部声明为准 + - 当前 style 内部声明的 sprite 是 `/newpec/sprite-semantic/sprite` + - 当前 style 内部声明的 PBF 是 `/pbf-delivery-full-20260418-rebuild` +- 已用 headless Chrome 自检： + - `variant=semantic¢er=130.070228,33.634637&zoom=12` + - 右侧大块空白已消失 + - 1536px 宽桌面下工具栏高度约 `57px` + - `点击对比` 面板展开 / 收起 / 拖拽均可用 + - 默认面板落位已切到底部居中，左右地图 pane 仍保持等宽 + - chip 显示 `HTML: full-audit-r7-20260501-1800` + - chip 显示 `PBF: full-rebuild-pbf-20260418-extentfix` + - 网络请求确认加载的是 `style.navsea-delivery-full-semantic-extentfix.json` + +### 2026-04-25 渔港 20m 续跑改为数据库优先，并新增当前 PRC 状态文件 + +- 已调整 `coastline/build_fish_port_20m_full_mysql_resume.py` 的续跑规则： + - `resume` 现在优先读取 MySQL 里的 `navsea_grid_import_state` 和 `navsea_grid_import_progress` + - 不再把本地 checkpoint 文件当成续跑主来源 + - 已完成的 `PRC_DONE` 会在 resume 时直接跳过 + - 如果用 `--prc` 指定了 PRC，则仍然允许重新跑指定 PRC +- 已新增当前 PRC 状态文件： + - `out/fish_port_20m_full_resume.current_prc.json` + - 每次开始一个 PRC 时都会写入当前 PRC、开始时间、已完成 PRC 列表等信息 + - 每个 PRC 完成时也会更新一次，方便人工查看当前跑到哪一段 +- 本地 checkpoint 文件仍保留作辅助记录： + - `out/fish_port_20m_full_resume.checkpoint.json` +- 已做最小验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + +### 2026-04-24 渔港 20m 海岸线改为 PRC 级缓存 + +- 已修正 `coastline/build_fish_port_20m_full_mysql_resume.py` 的主要慢点： + - 现在海岸线只在每个 `PRC` 开始时读取一次 + - 只在每个 `PRC` 开始时做一次 `unary_union` + - 只在每个 `PRC` 开始时做一次缓冲 `buffer` + - 每个 `part` 不再重复打开海岸线包、重复解析 GML、重复合并整包海岸线 +- 这次修正的目标是： + - 把“一个渔港要跑 30 分钟”的最主要重复成本先砍掉 + - 保留现有的 part 级扫描和 checkpoint 逻辑 +- 已做最小验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` +- 现场诊断补充： + - 当前这次运行没有写入 `out/fish_port_20m_full_resume.log` + - `navsea_grid_import_progress` 里也还没有新进度行 + - 说明它现在仍停留在单个 `part` 的几何扫描阶段，还没跑到 PRC 收口写库 + +### 2026-04-24 渔港 20m 再接 coast_200m 作为粗筛底座 + +- 已进一步把 `coastline/build_fish_port_20m_full_mysql_resume.py` 接到数据库里的 `coast_200m` 底库： + - 启动时先从 MySQL 读取 `coast_200m` + - 先把 `coast_200m` 做一次粗筛缓冲 + - 每个渔港 `part` 先用这层粗筛底座裁切一次 + - 再进入原有的 20m 细扫与海岸线判定 +- 这一步的目的： + - 让 20m 扫描只看“港区 + 海岸 200m 底座”覆盖到的范围 + - 先砍掉明显不相关的空白区域，再做精判 +- 已补充参数： + - `--coast-coarse-margin-m` +- 已做最小验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` + +### 2026-04-24 渔港 20m 内部 TRACE 加密 + +- 已把 `build_cells_for_geometry` 里的内部执行步骤打印出来： + - 准备阶段会打印 `fish_bounds`、扫描范围和 tile 范围 + - tile 命中时会打印 tile bbox + - row 命中时会打印 row y 坐标和横向范围 + - 候选格命中时会打印候选数量 + - 函数结束时会打印累计统计 +- 这样可以直接看出时间花在： + - tile 粗筛 + - row 细筛 + - 候选格精判 + - 还是最后写库 +- 已做最小验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` + +### 2026-04-24 渔港 20m 改为按 coast_200m 小格逐窗扫 + +- 已把 `coastline/build_fish_port_20m_full_mysql_resume.py` 再收紧一层： + - 不再只用整片 `coast_200m` union 去裁 `part` + - 改为直接按命中的 `coast_200m` 小格逐窗处理 + - 每个 200m 窗口内再跑原有 20m 细扫 +- 这一步的目标： + - 让 `candidate_cells` 真正随着海岸 200m 的局部区域缩小 + - 避免一个大 `part` 把几十万甚至上百万个 20m 候选格都扫进去 +- 已做最小验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` + +### 2026-04-24 渔港 20m 打印 coast_200m 命中格数 + +- 已继续给 `coastline/build_fish_port_20m_full_mysql_resume.py` 加日志： + - 每次按当前渔港 bbox 去数据库取 `coast_200m` 时，会打印取回的粗格窗口数 + - 这样可以直接看出这个 `part` 最终拿到了多少个 200m 粗筛格 +- 已做最小验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` + +### 2026-04-24 渔港 20m 对 200m 小窗启用直算模式 + +- 已给 `build_cells_for_geometry` 增加小窗直算快路径： + - 当输入几何本身已经缩到单个 `200m` 粗窗级别时 + - 不再再走 `2000m tile` 外循环 + - 直接按 `20m` 格线性扫描 +- 目的： + - 让 `candidate_cells` 与窗口尺寸严格对齐 + - 避免小窗还反复进入大范围 tile 扫描逻辑 +- 已做最小验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` + +### 2026-04-24 渔港 20m 结果页改为海面 / 陆基双层显示 + +- 已把 `src/pbf/navsea-fukuoka-saga-coast200-fish20-full.html` 改成双层渲染： + - `LAND_BASE` 用红色显示 + - `SEA_SURFACE` 用蓝色半透明显示 +- 已重新导出 MySQL 静态资产： + - `coastline/export_navgrid_mysql_assets.py --db-name navsea_fukuoka_saga_grid --fish-prc 40 41 --include-fish-sea` +- 这样页面不再只显示少量红色陆基块，而是能直接看到完整的渔港 20m 结果分布 + +### 2026-04-24 渔港 20m 红层提到海面层之上 + +- 已继续修正渲染顺序： + - 海面蓝层先画 + - 陆基红层后画并置顶 + - 红层透明度和描边略微加重 +- 这次修正是为了解决蓝色海面层把红色陆基层盖住的问题 +- 页面版本已升到 `v2.4` + +### 2026-04-24 全国 200x200 改成样区默认加载，红格已可见 + +- 已把全国海岸 200m 底库再导出一次，修正为经纬度 GeoJSON 输出： + - `src/pbf/coastline-mysql/japan_coast_200m/` + - `count = 336215` +- 已从全国底库额外导出九州样区： + - `src/pbf/coastline-mysql/japan_coast_200m_hakata/` + - `count = 3088` +- 现在预览页默认加载样区版，而不是全国全量版： + - 源文件：`src/pbf/navsea-coastline-fukuoka-saga-200m.html` + - 当前版本：`v1.5` +- 页面已经确认可以直接看到红色 200x200 格子 +- 全国全量仍保留在“全国概览”按钮里，方便后续切换回去看整体分布 + +### 2026-04-24 全国海岸 200x200 已从数据库导出并接入页面 + +- 已新增全国海岸 200x200 静态导出脚本： + - `coastline/export_japan_coast_grid_mysql_assets.py` +- 已从 MySQL `navsea_japan_coast_grid` 导出全国 200x200 数据： + - 输出目录：`src/pbf/coastline-mysql/japan_coast_200m/` + - 主要文件： + - `coast_200m_grid.geojson` + - `manifest.json` +- 全国导出数量： + - `336215` 个 200x200 硬阻塞格 +- 已把原来的红格预览页切换为全国视图： + - 源文件：`src/pbf/navsea-coastline-fukuoka-saga-200m.html` + - 线上页：`http://192.168.200.184/newpec/navsea-coastline-fukuoka-saga-200m.html` +- 当前页面版本： + - `v1.2` +- 已做一次本地 headless 视觉自检： + - 页面能正常加载全国红格 + - 视觉上可见全国海岸线周边的 200x200 红色覆盖 + - 页面初始视图已自动定位到全国范围 + +### 2026-04-24 全国 200x200 预览改成默认细节视图 + +- 由于全国总览缩放下 200m 格子太小，不容易直接看见，已把预览页改成默认显示九州海岸细节 +- 现在页面保留两个视图按钮： + - `九州细节` + - `全国概览` +- 当前页面版本已升到： + - `v1.3` +- 这样可以先确认红色格网确实存在，再切回全国概览看整体分布 + +### 2026-04-24 全国 200x200 红格再次加强可见性 + +- 已把同一份全国 200x200 数据的显示样式继续加重： + - 红色填充更实 + - 外轮廓更粗 + - 增加白色 halo 辅助识别 +- 当前页面版本已升到： + - `v1.4` +- 这次不改数据，只改视觉权重，目标是让红格在底图上更容易一眼看出来 + +### 2026-04-24 fish port 20m resume 脚本启动修复 + +- 已修复以下启动错误： + - `coastline/build_fish_port_20m_full_mysql_resume.py` + - 问题原因：`@trace_fn("load_xml")` 在 `trace_fn` 定义之前执行，导致模块加载时直接 `NameError` +- 已调整为先定义通用 tracing 装饰器，再装饰 `load_xml` +- 已做最小验证： + - `./.venv/bin/python -m py_compile coastline/build_fish_port_20m_full_mysql_resume.py` + - `./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --help` +- 当前状态： + - 脚本可以正常完成参数解析 + - 后续可以继续按原命令跑断点恢复任务 + +### 2026-04-24 渔港 20m 全国模式与指定 PRC 模式分离 + +- 已把 `coastline/build_fish_port_20m_full_mysql_resume.py` 的运行语义拆成两种： + - 默认模式：全国全量重建，不再把 `PRC` 当作必选条件 + - `--prc` 模式：只重建指定 PRC，并且忽略 `resume` 断点 +- 指定 `--prc` 时会先删除对应 PRC 已写入的旧格网，再按该 PRC 从头完整重建 +- 这样可以避免全国断点续跑和指定片区重建混在一起，导致日志看起来像“总是从 45 开始” + +### 2026-04-24 渔港 20m 格网扫描做了轻量提速 + +- 已优化 `coastline/build_fish_port_20m_full_mysql_resume.py` 的 `build_cells_for_geometry` +- 这次改动重点： + - 把每行新建的 `y_center/y_bottom/y_top_arr` 临时数组去掉了 + - 把 `x_values + half`、`x_values + cell_size_m` 这类重复数组改成按 tile 预先计算 + - 保留原有的 `intersects_xy` 过滤语义，不改判定结果 +- 这属于低风险提速，主要目标是减少大 PRC / 大 cluster 下的临时数组分配和内存抖动 + +### 2026-04-24 渔港 20m 扫描增加 profiling 汇总 + +- 已给 `coastline/build_fish_port_20m_full_mysql_resume.py` 的 `build_cells_for_geometry` 增加汇总 profiling +- 当前会记录并在每个 PRC 完成后输出： + - `tile_checked` + - `tile_hit` + - `row_checked` + - `row_hit` + - `mask_row_hit` + - `candidate_cells` + - `precise_cell_checks` + - `precise_hits` +- 这样可以直接看出一个 PRC 的时间到底花在： + - 粗 tile 过滤 + - 行级过滤 + - 候选格精判 + - 还是最终落库 +- 同时修正了 `profile` 的 PRC 级默认值，避免某些全跳过分支在收口日志处触发 `UnboundLocalError` + +### 2026-04-24 渔港 20m 进一步拆成单港部件处理 + +- 已把 `coastline/build_fish_port_20m_full_mysql_resume.py` 的处理粒度从“cluster 整体”再下沉一层 +- 当前流程变成： + - 先按 `PRC` 收集候选渔港线 + - 再按空间连通性聚成 cluster + - 然后把每个 cluster 的 `fish_geom` 拆成独立几何部件逐个处理 +- 这样单次扫格的 `bounds` 会更接近单个渔港，而不是整个 `PRC` 的大外框 +- 断点也同步细化到 `part` 级，避免中途停掉后还要重算整个 cluster + +### 2026-04-24 海岸线按单港 bbox 裁切 + +- 已将 `coastline/build_fish_port_20m_full_mysql_resume.py` 的海岸线读取进一步改为按当前渔港部件 `bbox` 裁切 +- 当前逻辑： + - 先拿单个渔港部件的 `bbox` + - 再去 `C23-06_{PRC}_GML.zip` 里只保留与该 `bbox` 相交的海岸线 + - 最终缓冲区也会再裁回这个 `bbox` +- 这样海岸线判定和 20m 格网扫描都不会再越过当前渔港的局部范围 +- 这符合“一个渔港算完进数据，然后下一个渔港”的执行口径 +- 另外给海岸线检索窗口加了很小的边界扩展，避免贴边海岸线因为刚好不穿过 `bbox` 而被误判为空；但最终格网仍只在原始 `bbox` 内计算 + +### 2026-04-24 全国海岸 200m 底库现状核对 + +- 已确认仓库里存在全国海岸 200m 的融合层定义： + - `out/navsea_national_navigation_fusion.sqlite` + - `layer_registry` 里已注册 `coastline_200m` +- 但当前这份库里的 `fusion_cell` 仍然是空的，说明全国海岸 200m 真值格还没有实际落库 +- 这意味着： + - 现阶段不能直接拿这份 national fusion 库来替代 `C23-06_*_GML.zip` + - 但后续一旦把全国海岸 200m 真值格填进去，这条 `PRC -> bbox -> 海岸格` 的路径就可以直接复用 + +### 2026-04-24 全国版海岸 200m 生成脚本已补齐 + +- 已新增全国版海岸 200m 生成脚本： + - `coastline/build_japan_coast_grid.py` +- 该脚本默认扫描： + - `coastline/C23-06_*_GML.zip` +- 当前处理方式是按单个 `C23-06` 包逐包生成 200m 底库，避免把全国海岸一次性合成一个巨大的几何对象 +- 默认输出： + - `out/coastline/japan_coast_grid/japan_coast_grid.sqlite` + - `out/coastline/japan_coast_grid/japan_coast_grid.meta.json` +- 可选输出： + - `japan_coastline.geojson` + - `japan_grid.geojson` +- 已做语法验证： + - `./.venv/bin/python -m py_compile coastline/build_japan_coast_grid.py` +- 说明： + - 当前仓库里的 MySQL 样例导入链路仍然是福冈/佐贺配置 + - 全国海岸 200m 底库已经可以先单独生成，后续再把 MySQL / 前端导入入口切到这份全国底库 + +### 2026-04-24 全国海岸 200m 直写 MySQL 脚本已新增 + +- 已新增全国海岸 200m 直接写 MySQL 的脚本： + - `coastline/build_japan_coast_grid_mysql.py` +- 该脚本会： + - 自动扫描 `coastline/C23-06_*_GML.zip` + - 逐包解析海岸线 GML + - 投影到米制坐标 + - 生成 200m 海岸硬阻塞格 + - 直接写入 MySQL +- 默认 MySQL 数据库名： + - `navsea_japan_coast_grid` +- 当前输出的主体表： + - `navsea_grid_layer_meta` + - `navsea_grid_package_stat` + - `navsea_grid_cell` +- 已支持可选的 manifest 输出，便于后续再接前端或导出脚本 + +### 2026-04-24 全国海岸 200m MySQL 口径收紧为只存阻塞格 + +- 已把 `coastline/build_japan_coast_grid_mysql.py` 的写库口径收紧为： + - 只写 `HARD_BLOCKED` + - 不再把 `NAVIGABLE_CANDIDATE` 之类的安全可航格写入 MySQL +- 这样全国海岸底库和之前福冈/佐贺测试版的口径一致： + - 底库只保留阻塞格 + - 候选安全格仅作为统计信息，不落库 +- 这次修正是按用户明确要求对齐，不再混入可航行候选格，避免污染后续全国融合层口径 + +### 2026-04-23 日本海岸线底库确认 + +- 已确认本地 `coastline/` 目录下存在日本海岸线原始数据包： + - `C23-06_*_GML.zip` +- 这些数据包为国土数值信息海岸线数据，包内包含： + - `*_g.xml` 主几何文件 + - `*_Coastline.shp/.dbf/.shx` 辅助矢量文件 + - `KS-META-*.xml` 元数据 +- 抽查结果显示： + - 几何是 `gml:Curve` + - 坐标是经纬度形式的 `lat lon` + - `gml:boundedBy` 覆盖日本海域范围 +- 这说明海岸线任务可以直接进入独立的路由底库预处理线： + - 先做海岸线导入与清洗 + - 再做米制投影 + - 再做 200m fishnet / `hardBlocked` 标记 +- 当前判断： + - 这条线与现有九州可航栅格工程是兼容的 + - 但它更适合独立成“全国海岸线硬阻塞底库”模块，而不是直接塞进现有 MVT / 显示流水线 + +### 2026-04-23 福冈 / 佐贺海岸样例栅格 + +- 已新增样例生成脚本： + - `coastline/build_fukuoka_saga_coast_grid.py` +- 已新增浏览器资产导出脚本： + - `coastline/export_fukuoka_saga_map_assets.py` +- 已新增可视化预览页： + - `src/pbf/navsea-coastline-fukuoka-saga.html` +- 已实际跑通输入包： + - `coastline/C23-06_40_GML.zip` + - `coastline/C23-06_41_GML.zip` +- 这两个包分别对应： + - 福冈县 + - 佐贺县 +- 输出样例库： + - `out/coastline/fukuoka_saga_test/fukuoka_saga_coast_grid.sqlite` +- 当前样例规模： + - 海岸线曲线数：`1067` + - 200m 格子数：`710562` + - `HARD_BLOCKED`：`9743` + - `NAVIGABLE_CANDIDATE`：`700819` +- 当前输出体积： + - 约 `118MB` +- 浏览器 GeoJSON 资产： + - `src/pbf/coastline/fukuoka_saga_blocked_cells.geojson` + - `src/pbf/coastline/fukuoka_saga_coastline_lines.geojson` + - `src/pbf/coastline/fukuoka_saga_manifest.json` +- 已部署到可直接打开的预览地址： + - `http://192.168.200.184/newpec/navsea-coastline-fukuoka-saga/` + - 目录入口：`index.html` +- 当前结论： + - 福冈 / 佐贺这条切片已经能作为 coastline 全国底库的第一版样例基线 +- 下一步可以直接沿着同一脚本扩到其他県包 + - 最近一次重跑已按低资源方式执行：`nice -n 10`，并将外部线程压到 1 + - 新输出目录： + - `out/coastline/fukuoka_saga_test_slow/` + - 结果与原样例一致： + - 海岸线曲线数：`1067` + - 200m 格子数：`710562` + - `HARD_BLOCKED`：`9743` + - `NAVIGABLE_CANDIDATE`：`700819` + +### 2026-04-23 九州港湾 / 渔港叠加预览 + +- 已新增港湾叠加导出脚本： + - `coastline/export_port_overlay_kyushu_assets.py` +- 已新增港湾叠加预览页： + - `src/pbf/navsea-port-overlay-kyushu.html` +- 数据源： + - `out/navsea_kyushu_navigability_full_run.export_snapshot.sqlite` +- 已导出资产： + - `src/pbf/port-overlay/kyushu/port_overlay_kyushu_areas.geojson` + - `src/pbf/port-overlay/kyushu/port_overlay_kyushu_centers.geojson` + - `src/pbf/port-overlay/kyushu/port_overlay_kyushu_manifest.json` +- 当前规模： + - AOI 总数：`5907` + - 漁港：`3169` + - 一般港湾：`1830` + - 小港湾：`908` +- 已部署到可直接打开的预览地址： + - `http://192.168.200.184/newpec/coastline/navsea-port-overlay-kyushu/` +- 当前判断： + - 这版先把港湾 / 渔港 AOI 叠加视图做出来了 + - 后续若要做 200m / 50m 分类格网，可以直接在这层基础上继续切分 + - 但这版目前仍是基于九州工程库里的 `port_aoi` 近似预览，不是正式的 `C02-08.zip / C09-06.zip` 原始港湾/渔港底图 + +### 2026-04-23 港湾数据源纠正 + +- 已确认港湾/渔港正式原始数据应来自： + - `coastline/C02-08.zip`：国土数值信息（港湾）数据 + - `coastline/C09-06.zip`：国土数值信息（漁港）数据 +- 当前九州预览中的圆形 AOI 仅是工程库近似叠加，不应作为正式港湾区域边界 +- 下一步需要把港湾/渔港层从这两个原始包重新生成，再替换现有预览 + +### 2026-04-23 港湾正式源预览已接入 + +- 已新增正式源导出脚本： + - `coastline/export_official_port_overlay_kyushu_assets.py` +- 已将预览页切换为正式源： + - `src/pbf/navsea-port-overlay-kyushu.html` +- 已导出并部署正式资产： + - `src/pbf/port-overlay-official/kyushu/official_port_overlay_kyushu_lines.geojson` + - `src/pbf/port-overlay-official/kyushu/official_port_overlay_kyushu_points.geojson` + - `src/pbf/port-overlay-official/kyushu/official_port_overlay_kyushu_manifest.json` +- 当前规模： + - 线要素：`1965` + - 点要素：`1366` + - C02-08 港湾线：`721` + - C09-06 漁港线：`1244` +- 预览地址： + - `http://192.168.200.184/newpec/coastline/navsea-port-overlay-kyushu/` +- 说明： + - 这版已经不再使用 `port_aoi` 近似圆圈 + - 现在读取的是 `C02-08.zip / C09-06.zip` 正式源 + +### 2026-04-23 coastline-only 港湾预览已切换 + +- 已新增仅依赖 `coastline/` 原始包的干净脚本： + - `coastline/export_coastline_only_port_overlay_kyushu_assets.py` +- 已新增仅依赖 coastline 原始包的新页面： + - `src/pbf/navsea-port-overlay-coastline-only-kyushu.html` +- 已导出正式资产： + - `src/pbf/coastline-only/kyushu/coastline_only_port_overlay_kyushu_lines.geojson` + - `src/pbf/coastline-only/kyushu/coastline_only_port_overlay_kyushu_points.geojson` + - `src/pbf/coastline-only/kyushu/coastline_only_port_overlay_kyushu_manifest.json` +- 当前规模： + - 线要素：`1965` + - 点要素：`1366` +- 已部署到可直接打开的预览地址： + - `http://192.168.200.184/newpec/coastline/navsea-port-overlay-coastline-only-kyushu/` +- 说明： + - 这版不碰数据库 + - 只用 `coastline/C02-08.zip` 和 `coastline/C09-06.zip` + - 作为后续正式港湾 / 渔港分类栅格的干净基线 + +### 2026-04-23 福冈 / 佐贺渔港 50m 栅格预览 + +- 已新增仅使用渔港数据的 50m 栅格导出脚本： + - `coastline/export_fish_port_grid_fukuoka_saga_assets.py` +- 已新增对应浏览器预览页： + - `src/pbf/navsea-fish-port-grid-fukuoka-saga.html` +- 数据源： + - 仅使用 `coastline/C09-06.zip` + - 仅保留 `PRC=40/41` +- 当前导出结果： + - 渔港边界线：`128` + - 50m 栅格：`5772` + - `polygonize` 成功，未走线缓冲回退 +- 已部署到可直接打开的预览地址： + - `http://192.168.200.184/newpec/coastline/navsea-fish-port-grid-fukuoka-saga/` +- 当前判断： + - 这版已经去掉港湾层和中心点，只保留福冈 / 佐贺渔港 + - 适合作为后续渔港 50m 分类栅格的人工核查样例 + +### 2026-04-23 福冈 / 佐贺渔港 海面 / 陆基分层预览 + +- 已新增港内导航分层脚本： + - `coastline/export_fish_port_navsplit_fukuoka_saga_assets.py` +- 已新增对应浏览器预览页： + - `src/pbf/navsea-fish-port-navsplit-fukuoka-saga.html` +- 数据源： + - `coastline/C09-06.zip` 渔港边界 + - `coastline/C23-06_40_GML.zip` + - `coastline/C23-06_41_GML.zip` +- 分类目标： + - `LAND_BASE` + - `SEA_SURFACE` +- 当前导出结果： + - 渔港边界线：`128` + - 海岸线参考：`1067` + - 50m 格网：`5772` + - 其中陆基：`737` + - 海面：`5035` +- 已部署到可直接打开的预览地址： + - `http://192.168.200.184/newpec/coastline/navsea-fish-port-navsplit-fukuoka-saga/` +- 当前判断： + - 这版已经不是单纯边界线，而是能在渔港范围内直接看出海面和陆基的分层 + - 适合作为后续导航路径规划的港内判定样例 + +### 2026-04-23 渔港层口径修正为 20m 不可航行格 + +- 已确认用户最新要求： + - 海岸线继续使用 `C23-06_40/41` 做 `200m` 级别底座 + - 渔港层必须独立于海岸线，只使用 `C09-06.zip` + - 渔港层改为 `20m` 小格，不再用海岸线去推导港内可达性 + - 渔港格子的语义改成“不可航行格”，用于导航底图的高密度港区阻断层 +- 已据此重做渔港三层预览脚本： + - `coastline/render_fukuoka_saga_three_layer_png.py` +- 已导出新版三层 PNG： + - `src/pbf/coastline-only/fukuoka_saga_three_layer/fukuoka_saga_three_layer_v7.png` +- 当前三层口径： + - 海岸线：`200m` 红色 + - 渔港：`20m` 不可航行格 +- PBF 海上障碍：`50m` 黄色 + +### 2026-04-23 渔港 20m 重建脚本改为按 PRC 流式入库 + +- 已调整全国渔港重建脚本： + - `coastline/build_fish_port_20m_full_mysql_resume.py` +- 这次改动的目标： + - 不再一次性保留全部渔港线到内存 + - 改为先发现 PRC 列表，再按 PRC 一段一段收线、计算、入库 + - 每完成一个 PRC / cluster 就立即写 MySQL 并保存 checkpoint +- 当前运行方式： + - 继续使用 `--resume` + - 通过 checkpoint 接着上次 40 / 41 的验证断点往后跑 +- 当前处理策略： + - 渔港线只保留当前 PRC 的数据 + - 海岸线按 PRC 现算现用 + - 目标是把 swap 压低，避免全国版把内存顶满 + +### 2026-04-24 MySQL 栅格导出到 HTML 静态资产 + +- 已新增 MySQL 栅格静态资产导出脚本： + - `coastline/export_navgrid_mysql_assets.py` +- 用途： + - 从 `navsea_fukuoka_saga_grid.navsea_grid_cell` 导出页面使用的 GeoJSON + - 输出到 `src/pbf/coastline-mysql/fukuoka_saga/` + - 再同步到 `newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前福冈 / 佐贺测试页导出口径： + - 海岸 `coast_200m`：全量，`9743` + - 渔港 `fish_port_20m`：只导出 `PRC=40/41` 且 `LAND_BASE`，`2904` + - 海上障碍 `hazard_50m`：全量，`176454` +- 当前测试页： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` + +### 2026-04-24 渔港重建脚本增加心跳输出 + +- 已给 `coastline/build_fish_port_20m_full_mysql_resume.py` 加入运行心跳： + - 每累计固定数量的 cell 就写一次库 + - 每次写库都会打印 `HEARTBEAT` + - 输出包含： + - `PRC` + - `cluster` + - 当前累计写入格数 + - 已用时间 +- 目的： + - 避免长跑任务在中间阶段看起来像“死掉” + - 让你能从日志直接确认它还在持续推进 + +### 2026-04-24 渔港重建状态改为按 PRC 完成输出 + +- 已把 `coastline/build_fish_port_20m_full_mysql_resume.py` 的日志粒度改小： + - 取消批次级心跳输出 + - 每个 PRC 完成后只输出一条状态 + - 这条状态同时写入数据库进度表 +- 新增进度表： + - `navsea_grid_import_progress` +- 写入内容： + - `job_name` + - `layer_name` + - `prc` + - `status_name` + - `cluster_count` + - `inserted_cells` + - `land_cells` + - `sea_cells` + - `elapsed_sec` +- 这样你手工看日志时不会被一堆中间行干扰，也能在数据库里查到每个渔港段的完成状态 + +### 2026-04-24 渔港重建日志改为更密集的屏幕输出 + +- 已再次加密 `coastline/build_fish_port_20m_full_mysql_resume.py` 的屏幕打印： + - `PRC` 开始时打印一次 + - `cluster` 开始时打印一次 + - 每次批量提交时打印一次 + - `PRC` 完成时仍打印一次总结 +- 这样长跑时你能从终端里更容易判断它不是卡住，而是在持续推进 +- 业务写库节奏不变，还是按固定批量提交，并在数据库里保留进度记录 + +### 2026-04-24 渔港重建日志再细分到处理步骤 + +- 已继续把 `coastline/build_fish_port_20m_full_mysql_resume.py` 的屏幕打印拆细： + - `PRC` 开始 + - 读取海岸线 + - 海岸缓冲完成 + - cluster 数量 + - 每个 cluster 开始 + - cluster 合并渔港线 + - 生成格网 + - 写库完成 +- 这样你在终端里可以直接看出卡在哪一步，不会只看到一个“开始”和一个“结束” + +### 2026-04-24 渔港重建日志统一带时间戳 + +- 已给 `coastline/build_fish_port_20m_full_mysql_resume.py` 的状态输出统一加时间戳 +- 当前日志格式： + - `[YYYY-MM-DDTHH:MM:SS] [fish_port_20m] ...` +- checkpoint 语义不变： + - 当前断点仍是 `last_prc = 45` + - 所以下一次 `--resume` 会从 `46` 继续 + +### 2026-04-24 渔港 20m 格网生成提速 + +- 已优化 `coastline/build_fish_port_20m_full_mysql_resume.py` 的格网扫描逻辑： + - 不再纯粹逐格创建 `box()` 后再判断 + - 先按整行做裁剪 + - 再用 `intersects_xy` 批量筛候选格子 + - 最后只对候选格子做精确几何判断 +- 目标： + - 降低 20m 扫格阶段的几何对象创建数量 + - 缩短每个 cluster 的格网生成时间 + +### 2026-04-24 渔港 20m 格网继续按扫描块分段 + +- 已把 `coastline/build_fish_port_20m_full_mysql_resume.py` 的格网扫描再切成更小的块： + - 新增 `--scan-tile-m` + - 默认按 `2000m` 级别分块扫描 + - 每个扫描块内再做 20m 网格判断 +- 这样大 cluster 不会一次把整个外接框都压在同一个扫描循环里 +- 目标是让单次停顿更短，也让日志更容易看出它正在持续推进 + +### 2026-04-24 渔港重建增加函数级耗时日志 + +- 已给 `coastline/build_fish_port_20m_full_mysql_resume.py` 的几个大函数加上统一 trace 日志： + - `load_xml` + - `discover_prcs` + - `collect_fish_lines_for_prc` + - `cluster_line_indices` + - `load_coastlines_for_prc` + - `build_cells_for_geometry` + - `main` +- 日志会记录： + - 函数名 + - 入参摘要 + - 开始时间 + - 完成耗时 + - 异常信息 +- 目的： + - 后面可以直接定位最耗时的步骤 + - 方便分析是解析、聚类、格网扫描还是写库最慢 + +### 2026-04-23 福冈 / 佐贺三层格网改为 MySQL 主存储 + +- 用户明确要求：三层网格不要再用 SQLite，改为 MySQL 主库保存，GeoJSON 仅作为导出和预览产物 +- 已创建 MySQL 数据库：`navsea_fukuoka_saga_grid` +- 已建立两张表： + - `navsea_grid_layer_meta` + - `navsea_grid_cell` +- 已导入三层网格数据： + - 海岸线 `200m`：`9743` + - 渔港 `20m`：`2196` + - 海上危险 `50m`：`183339` +- 海上危险层只取 PBF 中的航行危险语义，不包含鱼礁 +- 已生成并写出 MySQL 导出资产： + - `src/pbf/coastline-mysql/fukuoka_saga/coast_200m_grid.geojson` + - `src/pbf/coastline-mysql/fukuoka_saga/fish_port_20m_grid.geojson` + - `src/pbf/coastline-mysql/fukuoka_saga/hazard_50m_grid.geojson` + - `src/pbf/coastline-mysql/fukuoka_saga/manifest.json` +- 已将预览页切到 MySQL 导出资产： + - `src/pbf/navsea-fukuoka-saga-coast200-fish20-full.html` + - 版本：`v2.0` +- 已部署到预览目录： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前判断： + - MySQL 已经成为这组福冈 / 佐贺三层格网的主存储 + - 后续需要导出或审核时，再从 MySQL 导出 GeoJSON 即可，不再依赖 SQLite 主库 + +### 2026-04-23 修正渔港 20m 未显示的问题 + +- 用户反馈页面上 20m 渔港格没有出来 +- 复查发现页面代码错误地按不存在的 `classification` 字段过滤渔港格 +- 实际导出字段为： + - `state_name` + - `class_name` +- 已修正页面逻辑： + - 只保留 `state_name === "LAND_BASE"` 或 `class_name === "LAND_BASE"` 的 20m 格 +- 已重新部署到： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前判断： + - 20m 渔港格本身是存在的 + - 之前没有显示纯粹是前端筛选字段写错，不是 MySQL 数据缺失 + +### 2026-04-23 神集岛危险物位置复核 + +- 用户指出神集岛附近一处渔网在 PBF 预览中明显偏北，和官方海图不一致 +- 已复核神集岛附近的同一 z12 tile（`12/3526/1642`）在多个交付根中一致： + - `pbf-delivery-kyushu-reencoded` + - `pbf-delivery-full-20260418-rebuild` + - `pbf-delivery-full-semantic-20260416` + - `pbf-delivery-full-20260415` +- 复核到的 `fixed_fishing_gear_area` 位置大致落在： + - 经度 `129.965 ~ 129.985` + - 纬度 `33.543 ~ 33.553` +- 对比神集岛附近参考点（神集岛ヘリポート）约为： + - 经度 `129.96885` + - 纬度 `33.54101` +- 这说明当前 PBF 源里的危险物位置确实位于岛北侧，而不是用户给出的海图中红框的东南侧 +- 当前判断： + - 这不是 MySQL 导入或 HTML 渲染造成的偏移 + - 也不是单个交付根的偶发问题 + - 更像是 PBF 危险物源数据本身与目标海图存在位置偏差 +- 后续如果要继续把危险物层用于导航规划，需要额外做源数据校正或改用更可信的危险物来源 +- 已把新版预览页部署到： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前结论： + - 渔港层已经从“海岸辅助红绿判定”切换为“港区自身 20m 阻断栅格” + - 后续如果要再细分港内通道，需要在这个 20m 纯港区底图之上再叠加更细规则，而不是回头用海岸线推导港内状态 + +### 2026-04-23 福冈 / 佐贺渔港 25m 海面 / 陆基分层预览 + +- 已新增 25m 版本脚本： + - `coastline/export_fish_port_navsplit_fukuoka_saga_assets.py` +- 已新增对应浏览器预览页： + - `src/pbf/navsea-fish-port-navsplit-fukuoka-saga-25m.html` +- 参数： + - `cell_size_m = 25` + - `coast_buffer_m = 15` +- 当前导出结果： + - 渔港边界线：`128` + - 海岸线参考：`1067` + - 50m 预览时的 5772 格，扩细后为 `22535` 格 + - 其中陆基：`1742` + - 海面：`20793` +- 已部署到可直接打开的预览地址： + - `http://192.168.200.184/newpec/coastline/navsea-fish-port-navsplit-fukuoka-saga-25m/` +- 当前判断： + - 25m 明显更适合看窄水道、泊位前缘和贴岸转折 + - 这版比 50m 更重，但对导航底库更有价值 + +### 2026-04-23 福冈 / 佐贺渔港 20m 海面 / 陆基分层预览 + +- 已新增 20m 版本脚本复用： + - `coastline/export_fish_port_navsplit_fukuoka_saga_assets.py` +- 已新增对应浏览器预览页： + - `src/pbf/navsea-fish-port-navsplit-fukuoka-saga-20m.html` +- 参数： + - `cell_size_m = 20` + - `coast_buffer_m = 12` +- 当前导出结果： + - 渔港边界线：`128` + - 海岸线参考：`1067` + - 20m 格网：`35026` + - 其中陆基：`2196` + - 海面：`32830` +- 已部署到可直接打开的预览地址： + - `http://192.168.200.184/newpec/coastline/navsea-fish-port-navsplit-fukuoka-saga-20m/` +- 当前判断： + - 20m 已经足够细，港内边角和窄通道会更清楚 + - 数据量继续上升，但仍然适合样区核查 + +### 2026-04-23 福冈 / 佐贺渔港 港内可通行 / 不可通行预览 + +- 已新增导航语义视图页面： + - `src/pbf/navsea-fish-port-navpass-fukuoka-saga-20m.html` +- 语义映射： + - `SEA_SURFACE -> PORT_PASSABLE` + - `LAND_BASE -> PORT_BLOCKED` +- 仍沿用 20m 分类格网数据，不重新算几何 +- 已部署到可直接打开的预览地址： + - `http://192.168.200.184/newpec/coastline/navsea-fish-port-navpass-fukuoka-saga-20m/` +- 当前判断： + - 这版更贴近导航场景，直接表达“能不能走” + - 后续如果认可这个口径，可以把该语义正式写入导出数据 + +### 2026-04-23 全国导航融合骨架 + +- 已新增全国融合计划文档： + - `tasks/route/NavSea_National_Navigation_Fusion_Plan_v1.md` +- 已新增全国融合骨架脚本： + - `navsea_national_navigation_fusion.py` +- 已写出模板配置： + - `tasks/route/NavSea_National_Navigation_Fusion_Plan_v1.template.json` +- 已初始化融合骨架 SQLite： + - `out/navsea_national_navigation_fusion.sqlite` +- 已写出骨架 manifest： + - `out/navsea_national_navigation_fusion.manifest.json` +- 当前判断： + - 这不是旧样区延伸，而是从零开始的全国融合容器 + - 后续可以按统一 schema 接入海岸线 200m、渔港 20m 和 PBF 障碍层 + - 海岸线输入已明确只使用 `coastline/` 原始包，不再混入 PBF 来源 + +### 2026-04-22 语义小库设计 + +- 明确了海上路径检索和危险物碰撞检测可以抽成一个派生语义小数据库 +- 新增文档： + - `tasks/pbf/NavSea_semantic_query_index_design.md` +- 当前推荐定位： + - PBF 负责完整载体和渲染 + - 小数据库负责语义索引、碰撞检测、路径检索 + - 两者通过 `fid` 和 `render_layer / source_layer_std` 互相回查 +- 进一步补充了移动设备落地口径： + - 推荐单文件 SQLite + - 用 `semantic_object + semantic_geometry + grid_cell` 三层结构 + - 格网层负责快速判定，几何层负责精确验证 + +### 2026-04-22 九州 SQLite 落库 + +- 已生成九州轻量语义 SQLite： + - `out/navsea_kyushu_semantic.sqlite` +- 对应生成脚本： + - `navsea_build_kyushu_sqlite.py` +- 当前内容规模： + - `semantic_object`：`1786316` + - `grid_cell`：`65147` +- 这版 SQLite 重点保留： + - `P陸域` / `P穴` + - `collision` + - `grounding` + - `navigation_mark` + - `boundary_reference` + - `route_reference` +- 这版不再追求完整渲染几何，目标是让移动端可以直接做航线规划和危险检测查询 + +### 2026-04-22 九州核心 SQLite 压缩版 + +- 进一步压缩生成了核心版： + - `out/navsea_kyushu_core.sqlite` +- 当前体积： + - `178MB` +- 当前内容规模： + - `semantic_object`：`1058445` + - `grid_cell`：`26231` +- 核心版只保留： + - `P陸域` + - `P穴` + - `collision` + - 浅水搁浅层 +- 核心版去掉了： + - `navigation_mark` + - `boundary_reference` + - `route_reference` + - 其他非必要辅助对象 +- 适合更紧的移动端包体和更快的本地查询 + +### 2026-04-22 最终发布字段清单 + +- 整理了接近发布版的字段边界说明： + - `tasks/pbf/NavSea_Final_Release_Field_List.md` +- 这份文档把字段分成了四类： + - 必留字段 + - 推荐保留字段 + - 工程版专属字段 + - 旧字段兼容集 +- 当前特别要盯住的点： + - `fid` 和 `render_layer` 不能在最终发布版里被误裁掉 + - `fid_legacy_raw`、`trace_status`、`source_layer_jp` 等应留在工程版或 trace 口径 + - 最终发布版应以 `fid + render_layer + canonical_* + chart_*` 为核心 + +### 2026-04-22 九州可航行性测试版试跑 + +- 新增测试脚本： + - `navsea_build_kyushu_navigability_test.py` +- 这条线当前目标不是“九州全海域一次性全精度栅格化”，而是先验证： + - `z12` 几何恢复到 Web Mercator 的流程可跑 + - `200m` 主层和港口 `50m` 细化层的 SQLite 结构可落 + - `tile default + sparse/refined grid` 的双层查询口径可用 +- 本轮试跑确认了两个关键现实： + - 按 feature 逐个落 `200m` 格子，写放大和内存放大会非常重，不适合九州全量直接跑 + - 更合理的测试口径应先收敛到： + - 九州全域保留 `z12 tile default=open_water` + - 只对港口触发 tile 做 `200m/50m` 细化 +- 这条线后续若要继续推进，优先方向应是： + - tile 级批处理 + - 港口/近岸 AOI 优先 + - 不再对纯开放水域做全量细格预烘焙 + - 现已把脚本收成可切换： + - `full-domain` + - `port-tiles` + +### 2026-04-22 路径规划快照回查核对 + +- 额外核对了 `P陸域` / `P穴` 的回查链路 +- 当前 `pbf_analysis` 中，通过 `feature_geojson_lookup` 关联后： + - `P陸域` 的 `fid` 不是空 + - `P穴` 的 `fid` 也不是空 +- 这说明： + - `landIdentifiers` 里出现大量 `fid=null`，更像是消费端快照组装时没有正确 join 回查表 + - 不是当前 `pbf` 源数据天然缺 `fid` +- 另一个容易混淆的点： + - `land_area` 不是 `feature_geojson_lookup` 里的原始 `vt_layer` 名 + - 它是语义/样式层名，真正应回查的是 `P陸域` / `P穴` 这类对象层 + +### 2026-04-23 移动端发布包口径 + +- 已明确当前 `12G` 级别九州全域库只能作为工程库，不适合直接下发移动端 +- 新增方案文档： + - `tasks/pbf/NavSea_移动端可航栅格压缩方案_2026-04-23.md` +- 已新增移动核心包导出脚本： + - `navsea_export_mobile_core_sqlite.py` +- 当前导出方式： + - 先从工程库做 SQLite 一致性快照 + - 再从快照里导出非 `open_water` 的稀疏核心格 +- 当前推荐发布口径： + - 工程全域库继续保留在服务端 + - 移动端核心包只保留稀疏 `200m` + - `50m` 只作为港口增量包 + +### 2026-04-21 航线规划 Q1 对接答复 + +- 基于当前仓库里的 builder、Domain、样式和对象保全审计，整理了： + - `tasks/route/Q1_answer.md` +- 这份答复明确区分了三类结论： + - 当前仓库已证实 + - 当前只能谨慎判断 + - 仍需上游 PBF 生产侧确认 +- 当前给航线规划的保守口径先定为： + - `z10-z12` 是建议工作区间 + - `z12` 是最保守可靠级别 + - `outline/boundary` 不作为 polygon 重建来源 + - `land_area` 是唯一陆海主判定层 + - 关键规划对象不能只看 style `minzoom`，要区分数据保留与显示起点 + ### 2026-04-16 图标改进与语义迁移 - 完成了新旧 sprite / 图标键的梳理与迁移准备 @@ -178,3 +1395,285 @@ - AOI 高 diff 主要仍是和原始 `newpec` 的整体视觉差异，不是旧坏片那种几何炸裂 - 备注： - 全量 `extent shift` 长扫耗时过大，已暂停，后续如需可改成更小范围抽样或直接复用几何审计结论 + +## 2026-04-23 福冈/佐贺海岸线预览已部署 + +- 已把福冈/佐贺海岸线 200m 栅格预览页部署到： + - `/mnt/sda1/www/newpec/coastline/navsea-coastline-fukuoka-saga/` +- 页面入口： + - `http://192.168.200.184/newpec/coastline/navsea-coastline-fukuoka-saga/` +- 页面资源： + - `index.html` + - `coastline/fukuoka_saga_manifest.json` + - `coastline/fukuoka_saga_blocked_cells.geojson` + - `coastline/fukuoka_saga_coastline_lines.geojson` +- 数据来源： + - `coastline/C23-06_40_GML.zip` + - `coastline/C23-06_41_GML.zip` +- 说明： + - 这版是只用 `coastline/` 原始包做的福冈/佐贺海岸线预览，不使用 PBF 海岸线数据 + +## 2026-04-23 福冈/佐贺海岸 200m + 渔港 20m 合并预览 + +- 新增合并预览页： + - `/mnt/sda1/www/newpec/coastline/navsea-fukuoka-saga-coast200-fish20/` +- 页面入口： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20/` +- 页面内容： + - 海岸线 `200m` 硬阻塞格网 + - 渔港 `20m` 海面/陆基格网 + - 海岸线轮廓与渔港边界同时显示 +- 数据来源： + - 海岸线：`coastline/C23-06_40_GML.zip`、`coastline/C23-06_41_GML.zip` + - 渔港：`coastline/C09-06.zip` + +## 2026-04-23 渔港范围内 20m 覆盖 200m + +- 已将合并预览页调整为： + - 渔港 `20m` 范围内不再显示海岸 `200m` 格子 + - 海岸 `200m` 只保留在渔港范围外作为骨架层 +- 这样更符合后续导航分层： + - 海岸层负责全国骨架 + - 渔港层负责港内细分 + +## 2026-04-23 渔港 20m 改为绿色显示 + +- 已将福冈/佐贺合并预览中的渔港 `20m` 图层改为绿色系 +- 海岸 `200m` 仍保持淡红色 +- 这样更容易在同一张图上区分： + - 红色：海岸骨架 + - 绿色：渔港港内格网 + +## 2026-04-23 福冈 / 佐贺渔港全量 raster 预览 + +- 已把福冈 / 佐贺全量渔港 AOI 渲染成 20m 分辨率的 raster 叠加图 +- 输出文件： + - `src/pbf/coastline-only/fukuoka_saga_fish_port_navsplit_full_20m_aoi/fish_port_navsplit_fukuoka_saga_full_20m_aoi.png` + - `src/pbf/coastline-only/fukuoka_saga_fish_port_navsplit_full_20m_aoi/fish_port_navsplit_fukuoka_saga_full_20m_aoi_raster_manifest.json` +- 预览页： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 说明： + - 这版是全量 AOI raster，不再把渔港 20m 数据铺成超大的 GeoJSON + - 更适合浏览器直接查看全范围覆盖情况 + +## 2026-04-23 福冈 / 佐贺三层预览口径 + +- 已把全量预览页调整为三层语义： + - 红色：海岸线 `200m` 硬阻塞 + - 绿色：渔港内 `20m` 不可航行区域 + - 黄色：海上的渔网 / 标识 / 障碍 +- 黄色层来源： + - `navsea_delivery` 向量瓦片中的海上障碍语义层 +- 当前预览页版本： + - `v1.1` + +## 2026-04-23 渔港重叠区改为红色 20m 小框 + +- 已根据最新确认把渔港层进一步收紧： + - 海岸线仍按 `200m` 红色骨架保留 + - 渔港层只显示与海岸重叠的 `20m` 红色小框 + - 不再使用绿色渔港块表达港内范围 +- 渔港层最新导出结果： + - `src/pbf/coastline-only/fukuoka_saga_three_layer/fukuoka_saga_three_layer_v8.png` + - `src/pbf/coastline-only/fukuoka_saga_three_layer/fukuoka_saga_three_layer_manifest.json` +- 页面已更新为： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前页面版本： + - `v1.5` +- 当前判断： + - 渔港层现在已经和海岸层解耦 + - 重叠区按渔港 `20m` 红色小框表达，更符合导航底图的视觉语义 + +## 2026-04-23 渔港层切回 20m GeoJSON 小格 + +- 已根据最新对照图重新修正渔港层表达： + - 不再使用 `v8.png` 那张 raster 贴图 + - 改为直接读取 `fish_port_navsplit_fukuoka_saga_20m_grid.geojson` + - 渔港层按 `20m` 小格绘制 + - 视觉上只保留红色小框，不再画绿色大块 +- 新页面版本： + - `v1.6` +- 新页面入口： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前判断： + - 这才是用户要的港内细格表达 + - 渔港层应由真正的 20m 格网构成，而不是整块 raster 贴图 + +## 2026-04-23 渔港 20m 范围内裁掉海岸 200m + +- 已继续修正合并预览的显示规则： + - 渔港 `20m` 小格出现的位置，海岸 `200m` 大格不再显示 + - 海岸格网只保留在渔港范围外 +- 页面版本更新为： + - `v1.7` +- 当前入口： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前判断： + - 这条规则更符合用户的直觉要求 + - 渔港层和海岸层在视觉上不再重叠抢位 + +## 2026-04-23 渔港层只保留不可航行红格 + +- 已继续收紧渔港层显示： + - 只保留 `LAND_BASE` 对应的红色 20m 小格 + - `SEA_SURFACE` 不再显示 + - 渔港范围内只看不可航行格 +- 页面版本更新为： + - `v1.8` +- 当前入口： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前判断： + - 这才符合“有小格子的地方只保留红格”的最终口径 + - 渔港层现在只表达阻断，不再表达可通行海面 + +## 2026-04-23 神集岛原始几何复核 + +- 已回到最原始的 `newpec` MVT 瓦片核对神集岛附近的渔具定置区： + - 瓦片路径：`/home/wwwroot/newpec/exported_auto/tile.mapple-on.jp__newpec-mvt-20260106__z___x___y_.pbf/tiles/12/3526/1642.pbf` + - 对应原始层：`P漁具定置箇所` +- 该瓦片里与神集岛相关的两个对象是： + - `fid=18003411` + - 顶点范围约 `lon 129.965292 ~ 129.966978` + - 顶点范围约 `lat 33.543513 ~ 33.545354` + - 落在岛西北侧/北侧近岸 + - `fid=18003509` + - 顶点范围约 `lon 129.976376 ~ 129.985278` + - 顶点范围约 `lat 33.547511 ~ 33.553343` + - 落在岛东北侧外海 +- 结论： + - 这两条原始几何都不在用户手绘红框的东南侧位置 + - 说明我们当前抓到的 `newpec` 原始瓦片，和用户手头“右图”所指的位置并不一致 + - 后续如果要把黄格落到红框处，不能只改栅格逻辑，必须换源、换对象，或者先确认红框对应的具体图层 / fid + +## 2026-04-23 福冈 / 佐贺预览切换为 semantic style 底图 + +- 已把当前三层栅格预览页的底图切换为： + - `http://192.168.200.184/newpec/domain/style.navsea-delivery-full-semantic.json` +- 更新后的页面仍保留三层格网叠加： + - 海岸 `200m` + - 渔港 `20m` + - 海上障碍 `50m` +- 页面版本已更新到： + - `v2.1` +- 当前页面文件： + - `src/pbf/navsea-fukuoka-saga-coast200-fish20-full.html` +- 当前部署入口： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` + +## 2026-04-23 神集岛渔具定置区偏移再次确认 + +- 在切换为 `style.navsea-delivery-full-semantic.json` 底图后，神集岛附近的黄色渔具定置网格仍然明显落在岛北侧/东北侧，而不是用户标出的东南红框位置。 +- 这进一步说明： + - 不是当前底图样式导致偏移 + - 不是栅格绘制逻辑导致偏移 + - 而是原始 `newpec` 瓦片里的渔具定置几何与用户参考图存在位置不一致 + +## 2026-04-23 神集岛坐标原点修正 + +- 重新核对 `newpec` 原始瓦片后发现，之前的判断漏掉了一个关键点： + - 这批瓦片的 tile 内几何 y 轴应按“底部原点”理解 + - 之前在若干转换脚本里误用了“顶部原点”的反向换算 +- 复核示例： + - `fid=18003509` + - 使用底部原点换算后，中心点约为 `129.980826, 33.532361` + - 这已经明显比之前的错误换算更接近用户图上的南侧位置 +- 已修正脚本： + - `coastline/build_fukuoka_saga_navgrid_mysql.py` + - `coastline/render_fukuoka_saga_three_layer_png.py` + - `navsea_mvt_to_geojson.py` +- 已重新导入福冈 / 佐贺 MySQL 三层格网： + - 海岸 `200m`：`9743` + - 渔港 `20m`：`2196` + - 海上障碍 `50m`：`176454` +- 当前页面版本： + - `v2.2` +- 结论： + - 之前把 tile y 轴方向看反了，这是这次渔网位置异常的主要根因 + - 现在应以新导入的三层格网为准，再继续核对神集岛附近的黄格位置 + +## 2026-04-23 全国渔港 20m 全量重建脚本已启动 + +- 新增脚本： + - `coastline/build_fish_port_20m_full_mysql_resume.py` +- 脚本职责： + - 直接读取 `coastline/C09-06.zip` + - 按 `PRC` 分组处理全国渔港 + - 使用对应的 `coastline/C23-06_{PRC}_GML.zip` 做海岸线叠加分类 + - 将 `fish_port_20m` 重新写入 MySQL + - 支持 `--reset-layer` 清空旧数据 + - 支持 `--resume` 断点续跑 +- 已验证样例： + - `PRC=40/41` 试跑成功 + - 福冈 / 佐贺先导结果：`35026` 个 20m 格子 +- 当前全国续跑状态： + - 已从断点 `last_prc=40` 继续执行 + - 运行命令：`env PYTHONUNBUFFERED=1 nice -n 10 ./.venv/bin/python coastline/build_fish_port_20m_full_mysql_resume.py --resume` + - 该任务会从 `PRC=41` 往后继续扫全国渔港 + +## 2026-04-23 内存与交换空间诊断 + +- 当前机器 swap 使用约 `1.4GiB`，但不是单一 Python 进程独占导致。 +- 现阶段主要常驻进程占用包括： + - VSCode / Pylance 相关 Node 进程，RSS 较高 + - MySQL 守护进程 + - 当前全国渔港重建 Python 进程 +- 当前全国渔港脚本自身 RSS 约几百 MB 级，属于可控范围，但如果后续把更多 PRC 的几何一次性装入内存，仍有继续上涨风险。 +- 后续若要进一步压内存，优先优化方向是： + - 按 PRC 分段加载，不要长期保留所有渔港线 + - 海岸线几何按 PRC 现算现丢 + - 更细粒度地分批提交 MySQL + +## 2026-04-24 福冈 / 佐贺三层预览 UI 可视性修正 + +- 已把 `navsea-fukuoka-saga-coast200-fish20-full.html` 的海岸 `200m` 颜色改为绿色，避免和渔港红格混淆。 +- 已给左侧预览面板增加收起 / 展开按钮，折叠后只保留标题和切换按钮，尽量把地图画面让出来。 +- 当前页面版本更新到： + - `v2.5` +- 当前部署入口： + - `http://192.168.200.184/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` + +## 2026-04-24 海岸 200m 遮罩逻辑回退 + +- 复核发现：页面之前把海岸 `200m` 图层按 `fish_port_20m` 的 bbox 做了遮罩，导致用户在当前港区附近看不到该层的大量格子。 +- 这不是坐标系错误，`coast_200m` 和 `fish_port_20m` 的 GeoJSON 坐标都还是经纬度。 +- 已回退前端遮罩逻辑，改成直接显示完整的 `coast_200m`。 +- 页面版本已继续更新到： + - `v2.6` + +## 2026-04-24 渔港 20m 外扩增强 + +- 当前测试里，`fish_port_20m` 还偏“贴线”，港区外围包络不够完整。 +- 已给 `coastline/build_fish_port_20m_full_mysql_resume.py` 增加渔港几何外扩参数： + - `--fish-buffer-m` +- 默认值已提升到 `120m`，用于把 20m 栅格更完整地围住港区，避免只剩边缘带状格子。 + +## 2026-04-24 渔港 20m 改为直接跟随 200m 粗格扫描 + +- 结合呼子港的实际结果，确认之前“先渔港线、再外扩、再切 20m”的办法会把湾内区域漏掉。 +- 现在改成： + - 先抓出当前渔港范围命中的 `coast_200m` 粗格 + - 直接把这些 `200m` 粗格作为 `20m` 细扫底盘 + - 再用同一套海岸线碰撞逻辑判定 `LAND_BASE / SEA_SURFACE` +- 这样 20m 的范围会更接近用户手工圈出的那几块红框，而不是只沿着港岸线出一圈带状格。 + +## 2026-04-24 呼子港重跑结果已更新 + +- 这次重跑后，`fish_port_20m` 总数已经从旧版的 `28726` 提升到 `126888`。 +- 当前按 PRC 统计： + - `PRC=40`：`12798` + - `PRC=41`：`114090` +- 呼子港附近局部 bbox 内也能看到明显的陆 / 海混合格，不再是只剩外围带状格。 +- 已重新导出并同步到 Web 目录： + - `/mnt/sda1/www/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` + +## 2026-04-24 福冈 / 佐贺结果导出与 Web 同步 + +- 已从 MySQL 重新导出三层静态数据： + - `coast_200m: 9743` + - `fish_port_20m: 28726` + - `hazard_50m: 176454` +- 导出时 `fish_port_20m` 采用 `LAND_BASE + SEA_SURFACE` 全量导出，便于一起看港区外围形状。 +- 已同步到 Web 目录： + - `/mnt/sda1/www/newpec/coastline/navsea-fukuoka-saga-coast200-fish20-full/` +- 当前页面版本： + - `v2.6` diff --git a/coastline/build_fish_port_20m_full_mysql_resume.py b/coastline/build_fish_port_20m_full_mysql_resume.py new file mode 100644 index 0000000..523c1db --- /dev/null +++ b/coastline/build_fish_port_20m_full_mysql_resume.py @@ -0,0 +1,1823 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import datetime as dt +import json +import math +import time +import zipfile +import xml.etree.ElementTree as ET +from functools import wraps +from dataclasses import dataclass +from pathlib import Path +from typing import Iterable + +import numpy as np +import pymysql +from shapely import get_parts, intersects_xy +from shapely.geometry import LineString, box +from shapely.ops import polygonize, unary_union +from shapely.prepared import prep + + +DB_NAME = "navsea_japan_coast_grid" +DB_USER = "root" +DB_PASSWORD = "2chi9ks2" +DB_HOST = "localhost" +DB_SOCKET = "/tmp/mysql.sock" + +FISH_SRC = "coastline/C09-06.zip" +COAST_DIR = Path("coastline") +LAYER_NAME = "fish_port_20m" +CELL_SIZE_M = 20.0 +COAST_BUFFER_M = 12.0 +COAST_COARSE_LAYER = "coast_200m" +DEFAULT_COAST_COARSE_MARGIN_M = 0.0 +DEFAULT_CHECKPOINT = "out/fish_port_20m_full_resume.checkpoint.json" +DEFAULT_CURRENT_PRC_FILE = "out/fish_port_20m_full_resume.current_prc.json" +RESUME_JOB_NAME = "fish_port_20m_full_mysql_resume" + +RADIUS = 6378137.0 +MAX_MERCATOR_LAT = 85.0511287798066 + + +@dataclass(frozen=True) +class FishLine: + prc: str + fpc: str | None + coords: list[tuple[float, float]] + bbox: tuple[float, float, float, float] + + +def mercator_x(lon: float) -> float: + return RADIUS * math.radians(lon) + + +def mercator_y(lat: float) -> float: + lat = max(min(lat, MAX_MERCATOR_LAT), -MAX_MERCATOR_LAT) + return RADIUS * math.log(math.tan(math.pi / 4.0 + math.radians(lat) / 2.0)) + + +def lon_from_mercator(x: float) -> float: + return math.degrees(x / RADIUS) + + +def lat_from_mercator(y: float) -> float: + return math.degrees(2.0 * math.atan(math.exp(y / RADIUS)) - math.pi / 2.0) + + +def local_name(tag: str) -> str: + return tag.split("}", 1)[-1] + + +def lonlat_from_text(text: str) -> tuple[float, float]: + parts = [p for p in text.replace(",", " ").split() if p] + if len(parts) < 2: + raise ValueError(f"invalid coordinate text: {text!r}") + lat = float(parts[0]) + lon = float(parts[1]) + return lon, lat + + +def log(message: str) -> None: + print(f"[{dt.datetime.now().isoformat(timespec='seconds')}] {message}") + + +def summarize_value(value): + if isinstance(value, Path): + return str(value) + if isinstance(value, ET.Element): + return f"" + if isinstance(value, dict): + return f"dict(len={len(value)})" + if isinstance(value, (list, tuple, set)): + return f"{type(value).__name__}(len={len(value)})" + if hasattr(value, "geom_type"): + try: + bounds = tuple(round(v, 3) for v in value.bounds) + except Exception: + bounds = None + return f"{getattr(value, 'geom_type', type(value).__name__)}(bounds={bounds})" + return repr(value) + + +def trace_fn(label: str | None = None): + def decorator(func): + @wraps(func) + def wrapper(*args, **kwargs): + name = label or func.__name__ + arg_bits = [summarize_value(arg) for arg in args] + kw_bits = [f"{key}={summarize_value(val)}" for key, val in kwargs.items()] + joined = ", ".join(arg_bits + kw_bits) + start = time.monotonic() + log(f"[TRACE] {name} 开始 {joined}") + try: + result = func(*args, **kwargs) + except Exception as exc: + log(f"[TRACE] {name} 异常 {type(exc).__name__}: {exc}") + raise + elapsed = time.monotonic() - start + log(f"[TRACE] {name} 完成 elapsed={elapsed:.2f}s") + return result + + return wrapper + + return decorator + + +@trace_fn("load_xml") +def load_xml(zip_path: Path) -> ET.Element: + with zipfile.ZipFile(zip_path) as zf: + xml_name = next(name for name in zf.namelist() if name.endswith(".xml") and "META" not in name) + return ET.fromstring(zf.read(xml_name)) + + +def build_point_lookup(root: ET.Element) -> dict[str, tuple[float, float]]: + lookup: dict[str, tuple[float, float]] = {} + for point in root.iter(): + if local_name(point.tag) != "GM_Point": + continue + point_id = point.attrib.get("id") + if not point_id: + continue + coord_text = None + for el in point.iter(): + if local_name(el.tag).endswith("DirectPosition.coordinate") and (el.text or "").strip(): + coord_text = el.text.strip() + break + if coord_text is None: + continue + lookup[point_id] = lonlat_from_text(coord_text) + return lookup + + +def extract_curve_coords(curve: ET.Element, point_lookup: dict[str, tuple[float, float]]) -> list[tuple[float, float]]: + coords: list[tuple[float, float]] = [] + for pos in curve.iter(): + tag = local_name(pos.tag) + if not tag.startswith("GM_Position"): + continue + direct = None + ref = None + for child in list(pos): + child_tag = local_name(child.tag) + if child_tag.endswith("DirectPosition.coordinate") and (child.text or "").strip(): + direct = child.text.strip() + break + if child_tag.endswith("PointRef.point"): + ref = child.attrib.get("idref") + if direct is not None: + coords.append(lonlat_from_text(direct)) + elif ref and ref in point_lookup: + coords.append(point_lookup[ref]) + return coords + + +def bbox_from_coords(coords: list[tuple[float, float]]) -> tuple[float, float, float, float]: + xs = [p[0] for p in coords] + ys = [p[1] for p in coords] + return min(xs), min(ys), max(xs), max(ys) + + +def align_floor(value: float, step: float) -> float: + return math.floor(value / step) * step + + +def align_ceil(value: float, step: float) -> float: + return math.ceil(value / step) * step + + +def expand_bbox(bbox: tuple[float, float, float, float], margin_m: float) -> tuple[float, float, float, float]: + return (bbox[0] - margin_m, bbox[1] - margin_m, bbox[2] + margin_m, bbox[3] + margin_m) + + +def mysql_connect(database: str | None = None): + kwargs = { + "host": DB_HOST, + "user": DB_USER, + "password": DB_PASSWORD, + "charset": "utf8mb4", + "autocommit": False, + "cursorclass": pymysql.cursors.Cursor, + } + if DB_SOCKET: + kwargs["unix_socket"] = DB_SOCKET + if database: + kwargs["database"] = database + return pymysql.connect(**kwargs) + + +def ensure_database() -> None: + conn = mysql_connect() + try: + with conn.cursor() as cur: + cur.execute( + f"CREATE DATABASE IF NOT EXISTS `{DB_NAME}` " + "DEFAULT CHARACTER SET utf8mb4 COLLATE utf8mb4_unicode_ci" + ) + conn.commit() + finally: + conn.close() + + +def ensure_schema(conn) -> None: + stmts = [ + """ + CREATE TABLE IF NOT EXISTS navsea_grid_import_state ( + job_name VARCHAR(64) NOT NULL PRIMARY KEY, + layer_name VARCHAR(32) NOT NULL, + source_name VARCHAR(191) NOT NULL, + checkpoint_json LONGTEXT NOT NULL, + updated_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP ON UPDATE CURRENT_TIMESTAMP + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + """ + CREATE TABLE IF NOT EXISTS navsea_grid_import_progress ( + id BIGINT UNSIGNED NOT NULL AUTO_INCREMENT PRIMARY KEY, + job_name VARCHAR(64) NOT NULL, + layer_name VARCHAR(32) NOT NULL, + prc VARCHAR(8) NOT NULL, + status_name VARCHAR(32) NOT NULL, + cluster_count INT NOT NULL, + inserted_cells BIGINT NOT NULL, + land_cells BIGINT NOT NULL, + sea_cells BIGINT NOT NULL, + elapsed_sec DOUBLE NOT NULL, + created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP, + KEY idx_job_layer_prc (job_name, layer_name, prc), + KEY idx_job_created (job_name, created_at) + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + """ + CREATE TABLE IF NOT EXISTS navsea_grid_layer_meta ( + layer_name VARCHAR(32) NOT NULL PRIMARY KEY, + description VARCHAR(255) NOT NULL, + source_desc TEXT NOT NULL, + cell_size_m DOUBLE NOT NULL, + feature_count BIGINT NOT NULL, + bbox_min_lon DOUBLE NOT NULL, + bbox_min_lat DOUBLE NOT NULL, + bbox_max_lon DOUBLE NOT NULL, + bbox_max_lat DOUBLE NOT NULL, + export_file VARCHAR(255) DEFAULT NULL, + created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + """ + CREATE TABLE IF NOT EXISTS navsea_grid_cell ( + id BIGINT UNSIGNED NOT NULL AUTO_INCREMENT PRIMARY KEY, + layer_name VARCHAR(32) NOT NULL, + cell_id VARCHAR(64) NOT NULL, + row_idx INT NOT NULL, + col_idx INT NOT NULL, + cell_size_m DOUBLE NOT NULL, + state_name VARCHAR(32) NOT NULL, + class_name VARCHAR(32) NOT NULL, + source_name VARCHAR(191) NOT NULL, + min_lon DOUBLE NOT NULL, + min_lat DOUBLE NOT NULL, + max_lon DOUBLE NOT NULL, + max_lat DOUBLE NOT NULL, + created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP, + UNIQUE KEY uniq_layer_cell (layer_name, cell_id), + KEY idx_layer_state (layer_name, state_name), + KEY idx_layer_rowcol (layer_name, row_idx, col_idx) + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + ] + with conn.cursor() as cur: + for stmt in stmts: + cur.execute(stmt) + conn.commit() + + +def upsert_meta( + cur, + layer_name: str, + description: str, + source_desc: str, + cell_size_m: float, + feature_count: int, + bbox: tuple[float, float, float, float], + export_file: str | None, +) -> None: + cur.execute( + """ + INSERT INTO navsea_grid_layer_meta + (layer_name, description, source_desc, cell_size_m, feature_count, + bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat, export_file) + VALUES + (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s) + ON DUPLICATE KEY UPDATE + description=VALUES(description), + source_desc=VALUES(source_desc), + cell_size_m=VALUES(cell_size_m), + feature_count=VALUES(feature_count), + bbox_min_lon=VALUES(bbox_min_lon), + bbox_min_lat=VALUES(bbox_min_lat), + bbox_max_lon=VALUES(bbox_max_lon), + bbox_max_lat=VALUES(bbox_max_lat), + export_file=VALUES(export_file) + """, + ( + layer_name, + description, + source_desc, + cell_size_m, + feature_count, + bbox[0], + bbox[1], + bbox[2], + bbox[3], + export_file, + ), + ) + + +def insert_cells(cur, rows: list[tuple]) -> None: + if not rows: + return + cur.executemany( + """ + INSERT IGNORE INTO navsea_grid_cell + (layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, source_name, + min_lon, min_lat, max_lon, max_lat) + VALUES + (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s) + """, + rows, + ) + + +def save_checkpoint(path: Path, payload: dict) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + tmp_path = path.with_suffix(path.suffix + ".tmp") + tmp_path.write_text(json.dumps(payload, ensure_ascii=False, indent=2), encoding="utf-8") + tmp_path.replace(path) + + +def load_checkpoint(path: Path) -> dict | None: + if not path.exists(): + return None + return json.loads(path.read_text(encoding="utf-8")) + + +def normalize_prc_list(values: Iterable[str]) -> list[str]: + unique = {normalize_prc_token(str(value)) for value in values} + try: + return sorted(unique, key=lambda item: int(item)) + except ValueError: + return sorted(unique) + + +def normalize_prc_token(value: str) -> str: + value = str(value).strip() + if value.isdigit(): + return f"{int(value):02d}" + return value + + +def normalize_fpc_token(value: str) -> str: + return str(value).strip() + + +def load_resume_state_from_db(conn, job_name: str, layer_name: str) -> dict | None: + with conn.cursor() as cur: + cur.execute( + """ + SELECT checkpoint_json + FROM navsea_grid_import_state + WHERE job_name=%s AND layer_name=%s + """, + (job_name, layer_name), + ) + row = cur.fetchone() + if not row or row[0] is None: + return None + state = json.loads(row[0]) + cur.execute( + """ + SELECT prc + FROM navsea_grid_import_progress + WHERE job_name=%s AND layer_name=%s AND status_name='PRC_DONE' + ORDER BY id + """, + (job_name, layer_name), + ) + completed_prcs = [str(prc) for (prc,) in cur.fetchall()] + completed_prcs = normalize_prc_list(completed_prcs) + state["completed_prcs"] = completed_prcs + state.setdefault("job_name", job_name) + state.setdefault("layer_name", layer_name) + state.setdefault("status", "running") + return state + + +def load_resume_state(conn, _path: Path, job_name: str, layer_name: str) -> dict | None: + state = load_resume_state_from_db(conn, job_name, layer_name) + if state is not None: + return state + return None + + +def load_coarse_mask(conn, layer_name: str) -> tuple[object | None, tuple[float, float, float, float] | None, int]: + with conn.cursor() as cur: + cur.execute( + """ + SELECT min_lon, min_lat, max_lon, max_lat + FROM navsea_grid_cell + WHERE layer_name=%s + """, + (layer_name,), + ) + boxes = [box(float(min_lon), float(min_lat), float(max_lon), float(max_lat)) for min_lon, min_lat, max_lon, max_lat in cur] + + if not boxes: + return None, None, 0 + + geom = unary_union(boxes) + return geom, geom.bounds, len(boxes) + + +def load_coarse_windows_for_bbox( + conn, + *, + layer_name: str, + bbox_mercator: tuple[float, float, float, float], + margin_m: float, +) -> list[tuple[float, float, float, float]]: + minx, miny, maxx, maxy = bbox_mercator + query_minx = minx - margin_m + query_miny = miny - margin_m + query_maxx = maxx + margin_m + query_maxy = maxy + margin_m + + with conn.cursor() as cur: + cur.execute( + """ + SELECT min_lon, min_lat, max_lon, max_lat + FROM navsea_grid_cell + WHERE layer_name=%s + AND NOT ( + max_lon < %s OR min_lon > %s OR + max_lat < %s OR min_lat > %s + ) + ORDER BY row_idx, col_idx + """, + (layer_name, query_minx, query_maxx, query_miny, query_maxy), + ) + windows = [ + ( + float(row_min_lon), + float(row_min_lat), + float(row_max_lon), + float(row_max_lat), + ) + for row_min_lon, row_min_lat, row_max_lon, row_max_lat in cur + ] + log( + f"[TRACE] load_coarse_windows_for_bbox 完成 layer={layer_name} " + f"bbox={tuple(round(v, 3) for v in bbox_mercator)} margin={margin_m:.1f}m " + f"windows={len(windows)}" + ) + return windows + + +def delete_cells_by_ids(conn, layer_name: str, cell_ids: list[str]) -> None: + if not cell_ids: + return + with conn.cursor() as cur: + for offset in range(0, len(cell_ids), 1000): + chunk = cell_ids[offset : offset + 1000] + placeholders = ", ".join(["%s"] * len(chunk)) + cur.execute( + f"DELETE FROM navsea_grid_cell WHERE layer_name=%s AND cell_id IN ({placeholders})", + [layer_name, *chunk], + ) + conn.commit() + + +def drop_existing_layer(conn, layer_name: str) -> None: + with conn.cursor() as cur: + cur.execute("DELETE FROM navsea_grid_cell WHERE layer_name=%s", (layer_name,)) + cur.execute("DELETE FROM navsea_grid_layer_meta WHERE layer_name=%s", (layer_name,)) + cur.execute("DELETE FROM navsea_grid_import_state WHERE job_name=%s", (RESUME_JOB_NAME,)) + cur.execute("DELETE FROM navsea_grid_import_progress WHERE job_name=%s", (RESUME_JOB_NAME,)) + conn.commit() + + +def drop_existing_prcs(conn, layer_name: str, prcs: set[str]) -> None: + if not prcs: + return + source_names = [f"C09-06+coastline PRC={prc}" for prc in sorted(normalize_prc_list(prcs))] + placeholders = ", ".join(["%s"] * len(source_names)) + with conn.cursor() as cur: + cur.execute( + f"DELETE FROM navsea_grid_cell WHERE layer_name=%s AND source_name IN ({placeholders})", + [layer_name, *source_names], + ) + cur.execute( + f"DELETE FROM navsea_grid_import_progress WHERE job_name=%s AND prc IN ({placeholders})", + [RESUME_JOB_NAME, *sorted(prcs)], + ) + conn.commit() + + +def write_prc_status(path: Path, payload: dict) -> None: + save_checkpoint(path, payload) + + +def insert_progress( + cur, + *, + job_name: str, + layer_name: str, + prc: str, + status_name: str, + cluster_count: int, + inserted_cells: int, + land_cells: int, + sea_cells: int, + elapsed_sec: float, +) -> None: + cur.execute( + """ + INSERT INTO navsea_grid_import_progress + (job_name, layer_name, prc, status_name, cluster_count, inserted_cells, land_cells, sea_cells, elapsed_sec) + VALUES + (%s, %s, %s, %s, %s, %s, %s, %s, %s) + """, + ( + job_name, + layer_name, + prc, + status_name, + cluster_count, + inserted_cells, + land_cells, + sea_cells, + elapsed_sec, + ), + ) + + +def rect_lonlat(min_lon: float, min_lat: float, max_lon: float, max_lat: float) -> dict: + return { + "type": "Polygon", + "coordinates": [ + [ + [min_lon, min_lat], + [max_lon, min_lat], + [max_lon, max_lat], + [min_lon, max_lat], + [min_lon, min_lat], + ] + ], + } + + +def bbox_intersects(a: tuple[float, float, float, float], b: tuple[float, float, float, float]) -> bool: + return not (a[2] < b[0] or a[0] > b[2] or a[3] < b[1] or a[1] > b[3]) + + +@trace_fn("discover_prcs") +def discover_prcs(root: ET.Element, selected_prc: set[str] | None) -> list[str]: + obj = root.find(".//{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}OBJ") + if obj is None: + raise SystemExit("OBJ block missing in C09-06.zip") + + prcs: set[str] = set() + for feature in obj.findall(".//{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}CB03"): + prc_el = feature.find("{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}PRC") + prc = normalize_prc_token(prc_el.text) if prc_el is not None and (prc_el.text or "").strip() else "__unknown__" + if selected_prc is not None and prc not in selected_prc: + continue + prcs.add(prc) + return sorted(prcs) + + +@trace_fn("discover_target_prcs") +def discover_target_prcs( + root: ET.Element, + selected_prc: set[str] | None, + selected_fpc: str | None, +) -> list[str]: + if selected_fpc is None: + return discover_prcs(root, selected_prc) + + obj = root.find(".//{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}OBJ") + if obj is None: + raise SystemExit("OBJ block missing in C09-06.zip") + + prcs: set[str] = set() + for feature in obj.findall(".//{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}CB03"): + prc_el = feature.find("{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}PRC") + fpc_el = feature.find("{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}FPC") + prc = normalize_prc_token(prc_el.text) if prc_el is not None and (prc_el.text or "").strip() else "__unknown__" + fpc = normalize_fpc_token(fpc_el.text) if fpc_el is not None and (fpc_el.text or "").strip() else None + if selected_prc is not None and prc not in selected_prc: + continue + if fpc != selected_fpc: + continue + prcs.add(prc) + return sorted(prcs) + + +@trace_fn("collect_fish_lines_for_prc") +def collect_fish_lines_for_prc( + root: ET.Element, + point_lookup: dict[str, tuple[float, float]], + selected_prc: str, + selected_fpc: str | None = None, +) -> list[FishLine]: + obj = root.find(".//{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}OBJ") + if obj is None: + raise SystemExit("OBJ block missing in C09-06.zip") + + fish_lines: list[FishLine] = [] + for feature in obj.findall(".//{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}CB03"): + prc_el = feature.find("{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}PRC") + fpc_el = feature.find("{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}FPC") + prc = normalize_prc_token(prc_el.text) if prc_el is not None and (prc_el.text or "").strip() else "__unknown__" + if prc != selected_prc: + continue + fpc = normalize_fpc_token(fpc_el.text) if fpc_el is not None and (fpc_el.text or "").strip() else None + if selected_fpc is not None and fpc != selected_fpc: + continue + loc = feature.find("{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}LOC") + if loc is None: + continue + ref = loc.attrib.get("idref") + if not ref: + continue + curve = obj.find(f".//*[@id='{ref}']") + if curve is None: + continue + coords = extract_curve_coords(curve, point_lookup) + if len(coords) < 2: + continue + bbox = bbox_from_coords(coords) + fish_lines.append(FishLine(prc=prc, fpc=fpc, coords=coords, bbox=bbox)) + return fish_lines + + +@trace_fn("cluster_line_indices") +def cluster_line_indices(bboxes: list[tuple[float, float, float, float]], gap_m: float = 250.0) -> list[list[int]]: + parent = list(range(len(bboxes))) + + def find(x: int) -> int: + while parent[x] != x: + parent[x] = parent[parent[x]] + x = parent[x] + return x + + def union(a: int, b: int) -> None: + ra = find(a) + rb = find(b) + if ra != rb: + parent[rb] = ra + + expanded = [ + (bb[0] - gap_m, bb[1] - gap_m, bb[2] + gap_m, bb[3] + gap_m) + for bb in bboxes + ] + for i in range(len(expanded)): + for j in range(i + 1, len(expanded)): + if bbox_intersects(expanded[i], expanded[j]): + union(i, j) + + groups: dict[int, list[int]] = {} + for idx in range(len(parent)): + root = find(idx) + groups.setdefault(root, []).append(idx) + return list(groups.values()) + + +@trace_fn("load_coastlines_for_prc") +def load_coastlines_for_prc( + project_root: Path, + prc: str, + bbox: tuple[float, float, float, float] | None = None, + search_margin_m: float = 250.0, +) -> tuple[list[LineString], tuple[float, float, float, float]]: + bbox_geom = box(*expand_bbox(bbox, search_margin_m)) if bbox is not None else None + lines: list[LineString] = [] + bbox = [float("inf"), float("inf"), float("-inf"), float("-inf")] + coast_dir = project_root / COAST_DIR + exact_path = coast_dir / f"C23-06_{prc}_GML.zip" + candidate_paths = [exact_path] if exact_path.exists() else [] + fallback_paths = sorted(coast_dir.glob("C23-06_*_GML.zip")) + + if exact_path.exists(): + log(f"[TRACE] load_coastlines_for_prc 使用同号海岸线包 {exact_path.name}") + elif fallback_paths: + log( + f"[TRACE] load_coastlines_for_prc 缺少同号海岸线包 {exact_path.name}，" + f"将按 bbox 退回扫描 {len(fallback_paths)} 个现有海岸线包" + ) + candidate_paths = fallback_paths + else: + raise FileNotFoundError(str(exact_path)) + + tried_paths: list[Path] = [] + + def scan_path(coast_path: Path) -> None: + root = load_xml(coast_path) + point_lookup = build_point_lookup(root) + found_here = 0 + for curve in root.findall(".//{http://www.opengis.net/gml/3.2}Curve"): + coords: list[tuple[float, float]] = [] + for pos_list in curve.findall(".//{http://www.opengis.net/gml/3.2}posList"): + if not pos_list.text: + continue + values = [float(part) for part in pos_list.text.split() if part] + if len(values) < 4 or len(values) % 2 != 0: + continue + for i in range(0, len(values), 2): + lat = values[i] + lon = values[i + 1] + coords.append((mercator_x(lon), mercator_y(lat))) + if len(coords) < 2: + coords = [(mercator_x(lon), mercator_y(lat)) for lon, lat in extract_curve_coords(curve, point_lookup)] + if len(coords) < 2: + continue + line = LineString(coords) + if bbox_geom is not None and not line.intersects(bbox_geom): + continue + lines.append(line) + found_here += 1 + xs = [p[0] for p in coords] + ys = [p[1] for p in coords] + bbox[0] = min(bbox[0], min(xs)) + bbox[1] = min(bbox[1], min(ys)) + bbox[2] = max(bbox[2], max(xs)) + bbox[3] = max(bbox[3], max(ys)) + tried_paths.append(coast_path) + if found_here > 0: + log( + f"[TRACE] load_coastlines_for_prc 海岸线命中 source={coast_path.name} " + f"hits={found_here}" + ) + + for coast_path in candidate_paths: + scan_path(coast_path) + + if not lines and exact_path.exists(): + other_paths = [p for p in fallback_paths if p != exact_path] + if other_paths: + log( + f"[TRACE] load_coastlines_for_prc 同号包 {exact_path.name} 未命中 bbox，" + f"改用其他 {len(other_paths)} 个海岸线包继续搜索" + ) + for coast_path in other_paths: + scan_path(coast_path) + + if not lines: + raise SystemExit( + f"no coastline curves parsed for PRC={prc}; " + f"tried={len(tried_paths)} bbox={tuple(round(v, 3) for v in bbox_geom.bounds) if bbox_geom is not None else None}" + ) + return lines, (bbox[0], bbox[1], bbox[2], bbox[3]) + + +@trace_fn("build_cells_for_geometry") +def build_cells_for_geometry( + fish_geom, + coast_prepared, + *, + prc: str, + cell_size_m: float, + scan_tile_m: float, + source_name: str, + max_scan_cells: int | None, +) -> tuple[list[tuple], int, int, tuple[float, float, float, float], dict[str, float | int]]: + fish_prepared = prep(fish_geom) + minx, miny, maxx, maxy = fish_geom.bounds + start_x = align_floor(minx, cell_size_m) + start_y = align_floor(miny, cell_size_m) + end_x = align_ceil(maxx, cell_size_m) + end_y = align_ceil(maxy, cell_size_m) + estimated_cols = max(0, int(math.ceil((end_x - start_x) / cell_size_m))) + estimated_rows = max(0, int(math.ceil((end_y - start_y) / cell_size_m))) + estimated_cells = estimated_cols * estimated_rows + if max_scan_cells is not None and estimated_cells > max_scan_cells: + raise RuntimeError( + "scan guard triggered: " + f"estimated_cells={estimated_cells} max_scan_cells={max_scan_cells} " + f"grid={estimated_cols}x{estimated_rows} bounds={tuple(round(v, 3) for v in fish_geom.bounds)} " + f"source={source_name} prc={prc}" + ) + log( + f"[TRACE] build_cells_for_geometry 准备 " + f"fish_bounds={tuple(round(v, 3) for v in fish_geom.bounds)} " + f"scan_range=({start_x:.3f}, {start_y:.3f}, {end_x:.3f}, {end_y:.3f}) " + f"grid={estimated_cols}x{estimated_rows} estimated_cells={estimated_cells} " + f"cell={cell_size_m:.1f}m tile={scan_tile_m:.1f}m source={source_name}" + ) + + rows: list[tuple] = [] + land_count = 0 + sea_count = 0 + export_bbox = [float("inf"), float("inf"), float("-inf"), float("-inf")] + profile = { + "tile_checked": 0, + "tile_hit": 0, + "row_checked": 0, + "row_hit": 0, + "mask_row_hit": 0, + "candidate_cells": 0, + "precise_cell_checks": 0, + "precise_hits": 0, + } + + tile_start_x = align_floor(start_x, scan_tile_m) + tile_start_y = align_floor(start_y, scan_tile_m) + tile_end_x = align_ceil(end_x, scan_tile_m) + tile_end_y = align_ceil(end_y, scan_tile_m) + half = cell_size_m / 2.0 + log( + f"[TRACE] build_cells_for_geometry tile_range=" + f"({tile_start_x:.3f}, {tile_start_y:.3f}, {tile_end_x:.3f}, {tile_end_y:.3f})" + ) + + # 对 200m 粗窗这类很小的输入，直接按 20m 网格线性扫描， + # 避免再走一层 2000m tile 外循环。 + if (end_x - start_x) <= scan_tile_m and (end_y - start_y) <= scan_tile_m: + profile["tile_checked"] = 1 + profile["tile_hit"] = 1 + log(f"[TRACE] build_cells_for_geometry 使用小窗直算模式") + + y = start_y + while y < end_y: + profile["row_checked"] += 1 + row_strip = box(start_x, y, end_x, y + cell_size_m) + if not fish_prepared.intersects(row_strip): + y += cell_size_m + continue + profile["row_hit"] += 1 + if profile["row_hit"] <= 20 or profile["row_hit"] % 50 == 0: + log( + f"[TRACE] build_cells_for_geometry row_hit " + f"row_idx={profile['row_checked']} y={y:.3f} x_span=({start_x:.3f}, {end_x:.3f})" + ) + + ix_start = int(round(start_x / cell_size_m)) + ix_end = int(round(end_x / cell_size_m)) + ix_values = np.arange(ix_start, ix_end, dtype=np.int64) + x_values = ix_values.astype(float) * cell_size_m + x_center_values = x_values + half + x_right_values = x_values + cell_size_m + + mask = ( + intersects_xy(fish_geom, x_center_values, y + half) + | intersects_xy(fish_geom, x_values, y) + | intersects_xy(fish_geom, x_right_values, y) + | intersects_xy(fish_geom, x_values, y + cell_size_m) + | intersects_xy(fish_geom, x_right_values, y + cell_size_m) + ) + if not np.any(mask): + y += cell_size_m + continue + profile["mask_row_hit"] += 1 + if profile["mask_row_hit"] <= 20 or profile["mask_row_hit"] % 50 == 0: + log( + f"[TRACE] build_cells_for_geometry mask_hit " + f"row_idx={profile['row_checked']} candidate_mask={int(np.count_nonzero(mask))}" + ) + + mask = mask | np.r_[False, mask[:-1]] | np.r_[mask[1:], False] + candidate_indices = np.flatnonzero(mask) + profile["candidate_cells"] += int(candidate_indices.size) + if candidate_indices.size and (profile["candidate_cells"] <= 50 or profile["candidate_cells"] % 200 == 0): + log( + f"[TRACE] build_cells_for_geometry candidate_cells " + f"count={int(candidate_indices.size)} total={profile['candidate_cells']} " + f"row_y={y:.3f}" + ) + for idx in candidate_indices: + profile["precise_cell_checks"] += 1 + x = float(x_values[idx]) + cell_geom = box(x, y, x + cell_size_m, y + cell_size_m) + if not fish_prepared.intersects(cell_geom): + continue + profile["precise_hits"] += 1 + is_land = coast_prepared.intersects(cell_geom) + state_name = "LAND_BASE" if is_land else "SEA_SURFACE" + class_name = state_name + if is_land: + land_count += 1 + else: + sea_count += 1 + + min_lon = lon_from_mercator(x) + min_lat = lat_from_mercator(y) + max_lon = lon_from_mercator(x + cell_size_m) + max_lat = lat_from_mercator(y + cell_size_m) + export_bbox[0] = min(export_bbox[0], min_lon) + export_bbox[1] = min(export_bbox[1], min_lat) + export_bbox[2] = max(export_bbox[2], max_lon) + export_bbox[3] = max(export_bbox[3], max_lat) + ix = int(round(x / cell_size_m)) + iy = int(round(y / cell_size_m)) + cell_id = f"{ix}:{iy}" + rows.append( + ( + LAYER_NAME, + cell_id, + iy, + ix, + cell_size_m, + state_name, + class_name, + source_name, + min_lon, + min_lat, + max_lon, + max_lat, + ) + ) + y += cell_size_m + + if not rows: + log( + f"[TRACE] build_cells_for_geometry 结束 rows=0 land=0 sea=0 " + f"tile_hit={profile['tile_hit']}/{profile['tile_checked']} " + f"row_hit={profile['row_hit']}/{profile['row_checked']} " + f"candidate={profile['candidate_cells']} precise={profile['precise_hits']}/{profile['precise_cell_checks']}" + ) + return [], 0, 0, (float("inf"), float("inf"), float("-inf"), float("-inf")), profile + log( + f"[TRACE] build_cells_for_geometry 结束 rows={len(rows)} land={land_count} sea={sea_count} " + f"tile_hit={profile['tile_hit']}/{profile['tile_checked']} " + f"row_hit={profile['row_hit']}/{profile['row_checked']} " + f"candidate={profile['candidate_cells']} precise={profile['precise_hits']}/{profile['precise_cell_checks']}" + ) + return rows, land_count, sea_count, (export_bbox[0], export_bbox[1], export_bbox[2], export_bbox[3]), profile + + tile_y = tile_start_y + while tile_y < tile_end_y: + tile_x = tile_start_x + tile_top = min(tile_y + scan_tile_m, end_y) + while tile_x < tile_end_x: + profile["tile_checked"] += 1 + tile_right = min(tile_x + scan_tile_m, end_x) + tile_geom = box(tile_x, tile_y, tile_right, tile_top) + if not fish_prepared.intersects(tile_geom): + tile_x += scan_tile_m + continue + profile["tile_hit"] += 1 + log( + f"[TRACE] build_cells_for_geometry tile_hit " + f"idx={profile['tile_checked']} bbox=({tile_x:.3f}, {tile_y:.3f}, {tile_right:.3f}, {tile_top:.3f})" + ) + + local_start_x = max(start_x, tile_x) + local_end_x = tile_right + local_start_y = max(start_y, tile_y) + local_end_y = tile_top + + ix_start = int(round(local_start_x / cell_size_m)) + ix_end = int(round(local_end_x / cell_size_m)) + ix_values = np.arange(ix_start, ix_end, dtype=np.int64) + cell_count = int(ix_values.shape[0]) + if cell_count == 0: + tile_x += scan_tile_m + continue + x_values = ix_values.astype(float) * cell_size_m + x_center_values = x_values + half + x_right_values = x_values + cell_size_m + + y = local_start_y + while y < local_end_y: + profile["row_checked"] += 1 + row_strip = box(local_start_x, y, local_end_x, y + cell_size_m) + if not fish_prepared.intersects(row_strip): + y += cell_size_m + continue + profile["row_hit"] += 1 + if profile["row_hit"] <= 20 or profile["row_hit"] % 50 == 0: + log( + f"[TRACE] build_cells_for_geometry row_hit " + f"row_idx={profile['row_checked']} y={y:.3f} x_span=({local_start_x:.3f}, {local_end_x:.3f})" + ) + + mask = ( + intersects_xy(fish_geom, x_center_values, y + half) + | intersects_xy(fish_geom, x_values, y) + | intersects_xy(fish_geom, x_right_values, y) + | intersects_xy(fish_geom, x_values, y + cell_size_m) + | intersects_xy(fish_geom, x_right_values, y + cell_size_m) + ) + if not np.any(mask): + y += cell_size_m + continue + profile["mask_row_hit"] += 1 + if profile["mask_row_hit"] <= 20 or profile["mask_row_hit"] % 50 == 0: + log( + f"[TRACE] build_cells_for_geometry mask_hit " + f"row_idx={profile['row_checked']} candidate_mask={int(np.count_nonzero(mask))}" + ) + + mask = mask | np.r_[False, mask[:-1]] | np.r_[mask[1:], False] + candidate_indices = np.flatnonzero(mask) + profile["candidate_cells"] += int(candidate_indices.size) + if candidate_indices.size and (profile["candidate_cells"] <= 50 or profile["candidate_cells"] % 200 == 0): + log( + f"[TRACE] build_cells_for_geometry candidate_cells " + f"count={int(candidate_indices.size)} total={profile['candidate_cells']} " + f"row_y={y:.3f}" + ) + for idx in candidate_indices: + profile["precise_cell_checks"] += 1 + x = float(x_values[idx]) + cell_geom = box(x, y, x + cell_size_m, y + cell_size_m) + if not fish_prepared.intersects(cell_geom): + continue + profile["precise_hits"] += 1 + is_land = coast_prepared.intersects(cell_geom) + state_name = "LAND_BASE" if is_land else "SEA_SURFACE" + class_name = state_name + if is_land: + land_count += 1 + else: + sea_count += 1 + + min_lon = lon_from_mercator(x) + min_lat = lat_from_mercator(y) + max_lon = lon_from_mercator(x + cell_size_m) + max_lat = lat_from_mercator(y + cell_size_m) + export_bbox[0] = min(export_bbox[0], min_lon) + export_bbox[1] = min(export_bbox[1], min_lat) + export_bbox[2] = max(export_bbox[2], max_lon) + export_bbox[3] = max(export_bbox[3], max_lat) + ix = int(round(x / cell_size_m)) + iy = int(round(y / cell_size_m)) + cell_id = f"{ix}:{iy}" + rows.append( + ( + LAYER_NAME, + cell_id, + iy, + ix, + cell_size_m, + state_name, + class_name, + source_name, + min_lon, + min_lat, + max_lon, + max_lat, + ) + ) + y += cell_size_m + tile_x += scan_tile_m + + if not rows: + log( + f"[TRACE] build_cells_for_geometry 结束 rows=0 land=0 sea=0 " + f"tile_hit={profile['tile_hit']}/{profile['tile_checked']} " + f"row_hit={profile['row_hit']}/{profile['row_checked']} " + f"candidate={profile['candidate_cells']} precise={profile['precise_hits']}/{profile['precise_cell_checks']}" + ) + return [], 0, 0, (float("inf"), float("inf"), float("-inf"), float("-inf")), profile + log( + f"[TRACE] build_cells_for_geometry 结束 rows={len(rows)} land={land_count} sea={sea_count} " + f"tile_hit={profile['tile_hit']}/{profile['tile_checked']} " + f"row_hit={profile['row_hit']}/{profile['row_checked']} " + f"candidate={profile['candidate_cells']} precise={profile['precise_hits']}/{profile['precise_cell_checks']}" + ) + return rows, land_count, sea_count, (export_bbox[0], export_bbox[1], export_bbox[2], export_bbox[3]), profile + + +@trace_fn("main") +def main() -> None: + parser = argparse.ArgumentParser(description="Rebuild nationwide fish-port 20m grid into MySQL with resume support") + parser.add_argument( + "--fish-source", + default=FISH_SRC, + help="渔港原始包 C09-06.zip", + ) + parser.add_argument( + "--checkpoint-file", + default=DEFAULT_CHECKPOINT, + help="辅助断点文件（仅作记录）", + ) + parser.add_argument( + "--current-prc-file", + default=DEFAULT_CURRENT_PRC_FILE, + help="当前运行 PRC 的状态文件", + ) + parser.add_argument( + "--resume", + action="store_true", + default=True, + help="按数据库里的断点继续，默认开启；与 --prc 同用时会被忽略", + ) + parser.add_argument( + "--no-resume", + dest="resume", + action="store_false", + help="忽略数据库断点，从头规划但不清表", + ) + parser.add_argument( + "--reset-layer", + action="store_true", + help="先清空 fish_port_20m 的旧数据并删除断点", + ) + parser.add_argument( + "--prc", + nargs="+", + default=None, + help="只处理指定 PRC；指定时会按该 PRC 完整重建并忽略 resume 断点", + ) + parser.add_argument( + "--fpc", + default=None, + help="只处理指定 FPC（单个渔港）；指定时会先反查所在 PRC，再只重算该港口", + ) + parser.add_argument( + "--cluster-gap-m", + type=float, + default=250.0, + help="聚类连通判定的间隔（米）", + ) + parser.add_argument( + "--cell-size-m", + type=float, + default=CELL_SIZE_M, + help="格网大小（米）", + ) + parser.add_argument( + "--coast-buffer-m", + type=float, + default=COAST_BUFFER_M, + help="海岸缓冲距离（米）", + ) + parser.add_argument( + "--fish-buffer-m", + type=float, + default=120.0, + help="渔港几何外扩缓冲（米），保留兼容参数；当前默认扫描以 coast_200m 粗格为底盘", + ) + parser.add_argument( + "--commit-every", + type=int, + default=5000, + help="每累计多少条 cell 写一次数据库并打印状态", + ) + parser.add_argument( + "--scan-tile-m", + type=float, + default=2000.0, + help="格网扫描分块大小（米）", + ) + parser.add_argument( + "--coast-coarse-margin-m", + type=float, + default=DEFAULT_COAST_COARSE_MARGIN_M, + help="coast_200m 粗筛底座外扩缓冲（米），默认 0m（精确按同区域粗格筛选）", + ) + parser.add_argument( + "--allow-unbounded-fallback", + action="store_true", + help="允许在没有 coast_200m 小窗时回退扫描原始渔港 part；默认关闭以避免超大范围拖僵系统", + ) + parser.add_argument( + "--skip-missing-coarse", + action="store_true", + help="遇到没有有效 coast_200m 交集的 part 时跳过；默认作为异常中止，避免静默漏算", + ) + parser.add_argument( + "--max-scan-cells", + type=int, + default=200000, + help="单次 build_cells_for_geometry 允许估算扫描的最大 20m 格数，超过即中止；设为 0 表示不限制", + ) + args = parser.parse_args() + max_scan_cells = None if args.max_scan_cells <= 0 else args.max_scan_cells + + project_root = Path(__file__).resolve().parent.parent + fish_source = Path(args.fish_source) + if not fish_source.is_absolute(): + fish_source = project_root / fish_source + if not fish_source.exists(): + raise SystemExit(f"missing input: {fish_source}") + + checkpoint_path = Path(args.checkpoint_file) + if not checkpoint_path.is_absolute(): + checkpoint_path = project_root / checkpoint_path + current_prc_path = Path(args.current_prc_file) + if not current_prc_path.is_absolute(): + current_prc_path = project_root / current_prc_path + + ensure_database() + conn = mysql_connect(DB_NAME) + try: + ensure_schema(conn) + coast_coarse_geom, coast_coarse_bbox, coast_coarse_count = load_coarse_mask(conn, COAST_COARSE_LAYER) + if coast_coarse_geom is not None and coast_coarse_count > 0: + if args.coast_coarse_margin_m > 0: + coast_coarse_geom = coast_coarse_geom.buffer( + args.coast_coarse_margin_m, + cap_style=2, + join_style=2, + ) + coast_coarse_prepared = prep(coast_coarse_geom) + log( + f"[{LAYER_NAME}] coast_200m 粗筛底座就绪，cells={coast_coarse_count} " + f"bbox={coast_coarse_bbox} margin={args.coast_coarse_margin_m:.1f}m" + ) + else: + coast_coarse_geom = None + coast_coarse_prepared = None + if args.allow_unbounded_fallback: + log(f"[{LAYER_NAME}] coast_200m 粗筛底座为空，已允许回退按原始渔港范围扫描") + elif args.skip_missing_coarse: + log(f"[{LAYER_NAME}] coast_200m 粗筛底座为空，已允许跳过无粗筛窗口的 part") + else: + log(f"[{LAYER_NAME}] coast_200m 粗筛底座为空，后续遇到无粗筛窗口的 part 将作为异常中止") + selected_prc = None if args.prc is None else {normalize_prc_token(str(x)) for x in args.prc} + selected_fpc = normalize_fpc_token(args.fpc) if args.fpc else None + if selected_fpc is not None and args.resume: + log(f"[{LAYER_NAME}] 指定 FPC 模式，忽略 checkpoint 续跑，改为按指定渔港完整重建") + args.resume = False + if selected_prc is not None: + if args.resume: + log(f"[{LAYER_NAME}] 指定 PRC 模式，忽略 checkpoint 续跑，改为按指定 PRC 完整重建") + args.resume = False + log(f"[{LAYER_NAME}] 指定 PRC：{', '.join(sorted(selected_prc))}") + drop_existing_prcs(conn, LAYER_NAME, selected_prc) + if selected_fpc is not None: + log(f"[{LAYER_NAME}] 指定 FPC：{selected_fpc}") + if args.reset_layer: + log(f"[{LAYER_NAME}] 清空旧数据并删除 checkpoint ...") + drop_existing_layer(conn, LAYER_NAME) + if checkpoint_path.exists(): + checkpoint_path.unlink() + if current_prc_path.exists(): + current_prc_path.unlink() + + checkpoint = load_resume_state(conn, checkpoint_path, RESUME_JOB_NAME, LAYER_NAME) if args.resume else None + completed_prcs = set(normalize_prc_list(checkpoint.get("completed_prcs", []))) if checkpoint else set() + if checkpoint: + log( + f"[{LAYER_NAME}] 读取断点：" + f"last_prc={checkpoint.get('last_prc')} " + f"completed={len(completed_prcs)} " + f"current_prc={checkpoint.get('current_prc')}" + ) + + log(f"[{LAYER_NAME}] 开始解析渔港源：{fish_source}") + root = load_xml(fish_source) + point_lookup = build_point_lookup(root) + + prc_list = discover_target_prcs(root, selected_prc, selected_fpc) + if not prc_list: + if selected_fpc is not None: + raise SystemExit(f"no fish-port line features selected for FPC={selected_fpc}") + raise SystemExit("no fish-port line features selected") + if selected_fpc is not None: + log(f"[{LAYER_NAME}] 选中 PRC 数量：{len(prc_list)}（由 FPC 反查得到）") + else: + log(f"[{LAYER_NAME}] 选中 PRC 数量：{len(prc_list)}") + + total_inserted = 0 + total_land = 0 + total_sea = 0 + global_cluster_index = 0 + started_at = time.monotonic() + checkpoint_state = checkpoint or { + "job_name": RESUME_JOB_NAME, + "layer_name": LAYER_NAME, + "source_name": str(fish_source), + "status": "running", + "last_prc": None, + "last_cluster_index": -1, + "last_part_index": -1, + "processed_clusters": 0, + "inserted_cells": 0, + "land_cells": 0, + "sea_cells": 0, + "current_prc": None, + "current_prc_started_at": None, + "completed_prcs": [], + "updated_at": None, + } + + resume_prc = checkpoint.get("current_prc") if checkpoint else None + if resume_prc in completed_prcs: + resume_prc = None + resume_prc_index = prc_list.index(resume_prc) if resume_prc in prc_list else None + if args.resume and selected_prc is None: + pending_prcs = [ + prc + for prc in prc_list + if prc not in completed_prcs or prc == resume_prc + ] + if not pending_prcs: + log(f"[{LAYER_NAME}] 所有 PRC 都已完成，resume 不需要继续") + return + prc_list = pending_prcs + resume_prc_index = prc_list.index(resume_prc) if resume_prc in prc_list else None + + with conn.cursor() as cur: + for prc_index, prc in enumerate(prc_list): + if selected_prc is not None and prc not in selected_prc: + continue + if args.resume and selected_prc is None and prc in completed_prcs and prc != resume_prc: + continue + if resume_prc_index is not None and prc_index < resume_prc_index: + continue + if resume_prc is not None and prc == resume_prc and checkpoint and checkpoint.get("last_cluster_index", -1) >= 0: + log( + f"[{LAYER_NAME}] PRC={prc} 从断点恢复，" + f"last_cluster_index={checkpoint.get('last_cluster_index')}, " + f"last_part_index={checkpoint.get('last_part_index', -1)}" + ) + prc_lines = collect_fish_lines_for_prc(root, point_lookup, prc, selected_fpc=selected_fpc) + if not prc_lines: + if selected_fpc is not None: + log(f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} 没有渔港线要素，跳过") + else: + log(f"[{LAYER_NAME}] PRC={prc} 没有渔港线要素，跳过") + continue + if selected_fpc is not None: + log( + f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} 开始，" + f"线要素={len(prc_lines)}，已累计={checkpoint_state['inserted_cells']}" + ) + else: + log(f"[{LAYER_NAME}] PRC={prc} 开始，线要素={len(prc_lines)}，已累计={checkpoint_state['inserted_cells']}") + checkpoint_state.update( + { + "status": "running", + "current_prc": prc, + "current_prc_started_at": dt.datetime.now().isoformat(timespec="seconds"), + "completed_prcs": normalize_prc_list(completed_prcs), + "updated_at": dt.datetime.now().isoformat(timespec="seconds"), + } + ) + save_checkpoint(checkpoint_path, checkpoint_state) + write_prc_status( + current_prc_path, + { + "job_name": RESUME_JOB_NAME, + "layer_name": LAYER_NAME, + "status": "running", + "current_prc": prc, + "current_prc_started_at": checkpoint_state["current_prc_started_at"], + "last_prc": checkpoint_state.get("last_prc"), + "completed_prcs": normalize_prc_list(completed_prcs), + "updated_at": checkpoint_state["updated_at"], + **({"fpc": selected_fpc} if selected_fpc is not None else {}), + }, + ) + if selected_fpc is not None: + log(f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} 准备按单港 bbox 裁切海岸线并做缓冲 ...") + else: + log(f"[{LAYER_NAME}] PRC={prc} 准备按单港 bbox 裁切海岸线并做缓冲 ...") + + prc_bboxes = [line.bbox for line in prc_lines] + prc_bbox = ( + min(bb[0] for bb in prc_bboxes), + min(bb[1] for bb in prc_bboxes), + max(bb[2] for bb in prc_bboxes), + max(bb[3] for bb in prc_bboxes), + ) + cluster_groups = cluster_line_indices(prc_bboxes, gap_m=args.cluster_gap_m) + if selected_fpc is not None: + log(f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} cluster 数={len(cluster_groups)}") + else: + log(f"[{LAYER_NAME}] PRC={prc} cluster 数={len(cluster_groups)}") + + if selected_fpc is not None: + log(f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} 正在缓存整港海岸线并生成缓冲区 ...") + else: + log(f"[{LAYER_NAME}] PRC={prc} 正在缓存整县海岸线并生成缓冲区 ...") + prc_coast_lines, prc_coast_bbox = load_coastlines_for_prc( + project_root, + prc, + bbox=( + mercator_x(prc_bbox[0]), + mercator_y(prc_bbox[1]), + mercator_x(prc_bbox[2]), + mercator_y(prc_bbox[3]), + ), + search_margin_m=max(args.coast_buffer_m, args.coast_coarse_margin_m), + ) + if selected_fpc is not None: + log( + f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} 海岸线缓存完成，" + f"曲线数={len(prc_coast_lines)} bbox={prc_coast_bbox}" + ) + else: + log( + f"[{LAYER_NAME}] PRC={prc} 海岸线缓存完成，" + f"曲线数={len(prc_coast_lines)} bbox={prc_coast_bbox}" + ) + prc_coast_geom = unary_union(prc_coast_lines) + prc_coast_buffer = prc_coast_geom.buffer(args.coast_buffer_m, cap_style=2, join_style=2) + prc_coast_prepared = prep(prc_coast_buffer) + if selected_fpc is not None: + log(f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} 海岸缓冲缓存完成，开始处理各 part ...") + else: + log(f"[{LAYER_NAME}] PRC={prc} 海岸缓冲缓存完成，开始处理各 part ...") + + resume_cluster = 0 + if checkpoint and checkpoint.get("last_prc") == prc: + resume_cluster = int(checkpoint.get("last_cluster_index", -1)) + 1 + + prc_started_at = time.monotonic() + prc_inserted_before = checkpoint_state["inserted_cells"] + prc_land_before = checkpoint_state["land_cells"] + prc_sea_before = checkpoint_state["sea_cells"] + selected_port_rows: dict[str, tuple] | None = {} if selected_fpc is not None else None + selected_port_coarse_windows: set[tuple[float, float, float, float]] | None = set() if selected_fpc is not None else None + prc_profile = { + "tile_checked": 0, + "tile_hit": 0, + "row_checked": 0, + "row_hit": 0, + "mask_row_hit": 0, + "candidate_cells": 0, + "precise_cell_checks": 0, + "precise_hits": 0, + } + + for cluster_idx, cluster in enumerate(cluster_groups): + global_cluster_index += 1 + if checkpoint and checkpoint.get("last_prc") == prc and cluster_idx < resume_cluster: + continue + log(f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} 开始，成员={len(cluster)}") + log(f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1} 正在合并渔港线 ...") + + cluster_lines = [ + LineString([(mercator_x(lon), mercator_y(lat)) for lon, lat in prc_lines[i].coords]) + for i in cluster + ] + merged = unary_union(cluster_lines) + polygons = list(polygonize(merged)) + if polygons: + fish_geom = unary_union(polygons) + geom_mode = "polygonize" + else: + fish_geom = merged.buffer(args.cell_size_m, cap_style=2, join_style=2) + geom_mode = "buffer" + if fish_geom.is_empty: + log(f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx} 为空，跳过") + continue + log(f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1} 几何模式={geom_mode}，正在生成格网 ...") + + fish_parts = [part for part in get_parts(fish_geom) if not part.is_empty] + if not fish_parts: + log(f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1} 没有可拆分部件，跳过") + continue + + resume_part = 0 + if checkpoint and checkpoint.get("last_prc") == prc and checkpoint.get("last_cluster_index") == cluster_idx: + resume_part = int(checkpoint.get("last_part_index", -1)) + 1 + + for part_idx, fish_part in enumerate(fish_parts): + if checkpoint and checkpoint.get("last_prc") == prc and checkpoint.get("last_cluster_index") == cluster_idx and part_idx < resume_part: + continue + fish_bbox = fish_part.bounds + coarse_windows = load_coarse_windows_for_bbox( + conn, + layer_name=COAST_COARSE_LAYER, + bbox_mercator=fish_bbox, + margin_m=args.coast_coarse_margin_m, + ) + scan_windows: list[object] = [] + if coarse_windows: + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} 取回 coast_200m 格子={len(coarse_windows)}" + ) + for coarse_window in coarse_windows: + window_geom = box(*coarse_window) + scan_geom = fish_part.intersection(window_geom) + if not scan_geom.is_empty: + scan_windows.append(scan_geom) + if selected_port_coarse_windows is not None: + selected_port_coarse_windows.add(tuple(float(v) for v in coarse_window)) + if not scan_windows: + if args.allow_unbounded_fallback: + scan_windows = [fish_part] + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} 没有有效 coast_200m 交集，已允许回退原始范围" + ) + else: + message = ( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} 没有有效 coast_200m 交集；" + f"fish_bbox={tuple(round(v, 3) for v in fish_bbox)} " + f"coarse_windows={len(coarse_windows)}。" + ) + if args.skip_missing_coarse: + log(message + "已按 --skip-missing-coarse 跳过") + continue + raise RuntimeError( + message + + "按当前设计这属于异常：渔港 20m 必须以数据库 coast_200m 为底盘。" + + "如确需旧式整片扫描请显式加 --allow-unbounded-fallback；" + + "如确认允许漏算该 part，请显式加 --skip-missing-coarse。" + ) + else: + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} 实际扫描交集窗口={len(scan_windows)}" + ) + if selected_fpc is not None: + log( + f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} bbox={fish_bbox} 使用海岸线缓存 ..." + ) + else: + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} bbox={fish_bbox} 使用 PRC 级海岸线缓存 ..." + ) + coast_prepared = prc_coast_prepared + log(f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1} part={part_idx + 1} 开始格网 ...") + land_count = 0 + sea_count = 0 + bbox_lonlat = (float("inf"), float("inf"), float("-inf"), float("-inf")) + profile = {key: 0 for key in prc_profile} + part_inserted = 0 + batch: list[tuple] = [] + + for window_idx, scan_geom in enumerate(scan_windows): + if scan_geom.is_empty: + continue + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} window={window_idx + 1}/{len(scan_windows)} " + f"正在生成格网 ..." + ) + window_rows, window_land, window_sea, window_bbox, window_profile = build_cells_for_geometry( + scan_geom, + coast_prepared, + prc=prc, + cell_size_m=args.cell_size_m, + scan_tile_m=args.cell_size_m * 10.0, + source_name=( + f"C09-06+coastline PRC={prc}" + + (f" FPC={selected_fpc}" if selected_fpc is not None else "") + ), + max_scan_cells=max_scan_cells, + ) + land_count += window_land + sea_count += window_sea + bbox_lonlat = ( + min(bbox_lonlat[0], window_bbox[0]), + min(bbox_lonlat[1], window_bbox[1]), + max(bbox_lonlat[2], window_bbox[2]), + max(bbox_lonlat[3], window_bbox[3]), + ) + for key in prc_profile: + prc_profile[key] += window_profile[key] + profile[key] += window_profile[key] + + for row in window_rows: + if selected_port_rows is not None: + selected_port_rows[row[1]] = row + batch.append(row) + if len(batch) >= args.commit_every: + if selected_fpc is not None: + delete_cells_by_ids(conn, LAYER_NAME, [item[1] for item in batch]) + insert_cells(cur, batch) + conn.commit() + part_inserted += len(batch) + total_inserted += len(batch) + checkpoint_state["inserted_cells"] += len(batch) + checkpoint_state["updated_at"] = dt.datetime.now().isoformat(timespec="seconds") + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} " + f"提交={part_inserted} 累计={checkpoint_state['inserted_cells']} " + f"land={checkpoint_state['land_cells']} sea={checkpoint_state['sea_cells']}" + ) + batch.clear() + + if batch: + if selected_fpc is not None: + delete_cells_by_ids(conn, LAYER_NAME, [item[1] for item in batch]) + insert_cells(cur, batch) + conn.commit() + part_inserted += len(batch) + total_inserted += len(batch) + checkpoint_state["inserted_cells"] += len(batch) + checkpoint_state["updated_at"] = dt.datetime.now().isoformat(timespec="seconds") + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1}/{len(cluster_groups)} " + f"part={part_idx + 1}/{len(fish_parts)} " + f"提交={part_inserted} 累计={checkpoint_state['inserted_cells']} " + f"land={checkpoint_state['land_cells']} sea={checkpoint_state['sea_cells']}" + ) + + if part_inserted == 0: + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1} " + f"part={part_idx + 1} 没有生成 cell" + ) + checkpoint_state.update( + { + "last_prc": prc, + "last_cluster_index": cluster_idx, + "last_part_index": part_idx, + "processed_clusters": checkpoint_state["processed_clusters"] + 1, + "updated_at": dt.datetime.now().isoformat(timespec="seconds"), + } + ) + save_checkpoint(checkpoint_path, checkpoint_state) + continue + if selected_fpc is not None: + log( + f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} cluster={cluster_idx + 1} " + f"part={part_idx + 1} 格网完成，cell={part_inserted} land={land_count} sea={sea_count} bbox={bbox_lonlat}" + ) + log( + f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} cluster={cluster_idx + 1} " + f"part={part_idx + 1} 写库完成" + ) + else: + log( + f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1} " + f"part={part_idx + 1} 格网完成，cell={part_inserted} land={land_count} sea={sea_count} bbox={bbox_lonlat}" + ) + log(f"[{LAYER_NAME}] PRC={prc} cluster={cluster_idx + 1} part={part_idx + 1} 写库完成") + checkpoint_state.update( + { + "last_prc": prc, + "last_cluster_index": cluster_idx, + "last_part_index": part_idx, + "processed_clusters": checkpoint_state["processed_clusters"] + 1, + "land_cells": checkpoint_state["land_cells"] + land_count, + "sea_cells": checkpoint_state["sea_cells"] + sea_count, + "current_prc": prc, + "completed_prcs": normalize_prc_list(completed_prcs | {prc}), + "updated_at": dt.datetime.now().isoformat(timespec="seconds"), + } + ) + save_checkpoint(checkpoint_path, checkpoint_state) + write_prc_status( + current_prc_path, + { + "job_name": RESUME_JOB_NAME, + "layer_name": LAYER_NAME, + "status": "running", + "current_prc": prc, + "current_prc_started_at": checkpoint_state.get("current_prc_started_at"), + "last_prc": prc, + "last_cluster_index": cluster_idx, + "last_part_index": part_idx, + "completed_prcs": normalize_prc_list(completed_prcs | {prc}), + "inserted_cells": checkpoint_state["inserted_cells"], + "land_cells": checkpoint_state["land_cells"], + "sea_cells": checkpoint_state["sea_cells"], + "updated_at": checkpoint_state["updated_at"], + **({"fpc": selected_fpc} if selected_fpc is not None else {}), + }, + ) + total_land += land_count + total_sea += sea_count + + prc_elapsed = time.monotonic() - prc_started_at + prc_inserted_now = checkpoint_state["inserted_cells"] - prc_inserted_before + prc_land_now = checkpoint_state["land_cells"] - prc_land_before + prc_sea_now = checkpoint_state["sea_cells"] - prc_sea_before + summary_inserted = len(selected_port_rows) if selected_port_rows is not None else prc_inserted_now + summary_coarse = len(selected_port_coarse_windows) if selected_port_coarse_windows is not None else 0 + insert_progress( + cur, + job_name="fish_port_20m_full_mysql_resume", + layer_name=LAYER_NAME, + prc=prc, + status_name="FPC_DONE" if selected_fpc is not None else "PRC_DONE", + cluster_count=len(cluster_groups), + inserted_cells=summary_inserted, + land_cells=prc_land_now, + sea_cells=prc_sea_now, + elapsed_sec=prc_elapsed, + ) + conn.commit() + cur.execute( + """ + INSERT INTO navsea_grid_import_state + (job_name, layer_name, source_name, checkpoint_json) + VALUES + (%s, %s, %s, %s) + ON DUPLICATE KEY UPDATE + layer_name=VALUES(layer_name), + source_name=VALUES(source_name), + checkpoint_json=VALUES(checkpoint_json) + """, + ( + RESUME_JOB_NAME, + LAYER_NAME, + str(fish_source), + json.dumps(checkpoint_state, ensure_ascii=False), + ), + ) + conn.commit() + if selected_fpc is not None: + log( + f"[{LAYER_NAME}] PRC={prc} / FPC={selected_fpc} 完成 cluster={len(cluster_groups)} " + f"200x200={summary_coarse} final_20x20={summary_inserted} " + f"land={prc_land_now} sea={prc_sea_now} elapsed={prc_elapsed:.1f}s" + ) + else: + log( + f"[{LAYER_NAME}] PRC={prc} 完成 cluster={len(cluster_groups)} " + f"本段写入={prc_inserted_now} land={prc_land_now} sea={prc_sea_now} " + f"elapsed={prc_elapsed:.1f}s" + ) + log( + f"[{LAYER_NAME}] PRC={prc} profiling " + f"tile={prc_profile['tile_hit']}/{prc_profile['tile_checked']} " + f"row={prc_profile['row_hit']}/{prc_profile['row_checked']} " + f"mask_row={prc_profile['mask_row_hit']} " + f"candidate={prc_profile['candidate_cells']} " + f"precise={prc_profile['precise_hits']}/{prc_profile['precise_cell_checks']}" + ) + + completed_prcs.add(prc) + checkpoint_state["completed_prcs"] = normalize_prc_list(completed_prcs) + checkpoint_state["status"] = "idle" + checkpoint_state["current_prc"] = prc + checkpoint_state["current_prc_started_at"] = checkpoint_state.get("current_prc_started_at") + checkpoint_state["updated_at"] = dt.datetime.now().isoformat(timespec="seconds") + save_checkpoint(checkpoint_path, checkpoint_state) + write_prc_status( + current_prc_path, + { + "job_name": RESUME_JOB_NAME, + "layer_name": LAYER_NAME, + "status": "done", + "current_prc": prc, + "current_prc_started_at": checkpoint_state.get("current_prc_started_at"), + "completed_prcs": normalize_prc_list(completed_prcs), + "finished_at": checkpoint_state["updated_at"], + "updated_at": checkpoint_state["updated_at"], + **({"fpc": selected_fpc, "coarse_200m_intersections": summary_coarse, "final_20x20_cells": summary_inserted} if selected_fpc is not None else {}), + }, + ) + + del prc_lines + + cur.execute("SELECT COUNT(*) FROM navsea_grid_cell WHERE layer_name=%s", (LAYER_NAME,)) + final_count = int(cur.fetchone()[0]) + cur.execute( + """ + SELECT MIN(min_lon), MIN(min_lat), MAX(max_lon), MAX(max_lat) + FROM navsea_grid_cell + WHERE layer_name=%s + """, + (LAYER_NAME,), + ) + bbox_row = cur.fetchone() + bbox = ( + float(bbox_row[0]) if bbox_row and bbox_row[0] is not None else float("inf"), + float(bbox_row[1]) if bbox_row and bbox_row[1] is not None else float("inf"), + float(bbox_row[2]) if bbox_row and bbox_row[2] is not None else float("-inf"), + float(bbox_row[3]) if bbox_row and bbox_row[3] is not None else float("-inf"), + ) + upsert_meta( + cur, + LAYER_NAME, + "全国渔港 20m 可航/不可航格", + "coastline/C09-06.zip + navsea_japan_coast_grid.coast_200m", + args.cell_size_m, + final_count, + bbox, + None, + ) + cur.execute( + """ + INSERT INTO navsea_grid_import_state + (job_name, layer_name, source_name, checkpoint_json) + VALUES + (%s, %s, %s, %s) + ON DUPLICATE KEY UPDATE + layer_name=VALUES(layer_name), + source_name=VALUES(source_name), + checkpoint_json=VALUES(checkpoint_json) + """, + ( + RESUME_JOB_NAME, + LAYER_NAME, + str(fish_source), + json.dumps(checkpoint_state, ensure_ascii=False), + ), + ) + conn.commit() + + log(f"[{LAYER_NAME}] 完成，checkpoint={checkpoint_path}") + finally: + conn.close() + + +if __name__ == "__main__": + main() diff --git a/coastline/build_japan_coast_grid_mysql.py b/coastline/build_japan_coast_grid_mysql.py new file mode 100644 index 0000000..2ad8acb --- /dev/null +++ b/coastline/build_japan_coast_grid_mysql.py @@ -0,0 +1,565 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import datetime as dt +import glob +import json +import math +import zipfile +from dataclasses import dataclass +from pathlib import Path +from typing import Iterable +import xml.etree.ElementTree as ET + +import pymysql +from shapely.geometry import LineString, MultiLineString, box, mapping +from shapely.prepared import prep + + +DB_NAME = "navsea_japan_coast_grid" +DB_USER = "root" +DB_PASSWORD = "2chi9ks2" +DB_HOST = "localhost" +DB_SOCKET = "/tmp/mysql.sock" + +DEFAULT_INPUT_GLOB = "coastline/C23-06_*_GML.zip" +DEFAULT_OUT_DIR = "out/coastline/japan_coast_grid_mysql" +LAYER_NAME = "coast_200m" +CELL_SIZE_M = 200.0 + +RADIUS = 6378137.0 +MAX_MERCATOR_LAT = 85.0511287798066 +NS = {"gml": "http://www.opengis.net/gml/3.2"} + + +@dataclass(frozen=True) +class SourcePackage: + path: Path + prefecture: str + code: str + + +def mercator_x(lon: float) -> float: + return RADIUS * math.radians(lon) + + +def mercator_y(lat: float) -> float: + lat = max(min(lat, MAX_MERCATOR_LAT), -MAX_MERCATOR_LAT) + return RADIUS * math.log(math.tan(math.pi / 4.0 + math.radians(lat) / 2.0)) + + +def lon_from_mercator(x: float) -> float: + return math.degrees(x / RADIUS) + + +def lat_from_mercator(y: float) -> float: + return math.degrees(2.0 * math.atan(math.exp(y / RADIUS)) - math.pi / 2.0) + + +def mercator_bbox_to_lonlat(minx: float, miny: float, maxx: float, maxy: float) -> tuple[float, float, float, float]: + return ( + lon_from_mercator(minx), + lat_from_mercator(miny), + lon_from_mercator(maxx), + lat_from_mercator(maxy), + ) + + +def align_floor(value: float, step: float) -> float: + return math.floor(value / step) * step + + +def align_ceil(value: float, step: float) -> float: + return math.ceil(value / step) * step + + +def parse_text_list(text: str) -> list[float]: + return [float(part) for part in text.split() if part] + + +def detect_prefecture(meta_xml: str) -> str: + title_start = meta_xml.find("") + if title_start == -1: + return "unknown" + title_end = meta_xml.find("", title_start) + if title_end == -1: + return "unknown" + return meta_xml[title_start + 7:title_end].strip() + + +def load_package(path: Path) -> SourcePackage: + with zipfile.ZipFile(path) as zf: + meta_name = next(name for name in zf.namelist() if "META" in name and name.endswith(".xml")) + meta_xml = zf.read(meta_name).decode("shift_jis", errors="replace") + prefecture = detect_prefecture(meta_xml) + code = path.stem.replace("_GML", "") + return SourcePackage(path=path, prefecture=prefecture, code=code) + + +def iter_coastline_lines(zip_path: Path) -> Iterable[LineString]: + with zipfile.ZipFile(zip_path) as zf: + xml_name = next(name for name in zf.namelist() if name.endswith(".xml") and "META" not in name) + root = ET.fromstring(zf.read(xml_name)) + + for curve in root.findall(".//gml:Curve", NS): + coords: list[tuple[float, float]] = [] + for pos_list in curve.findall(".//gml:posList", NS): + if not pos_list.text: + continue + values = parse_text_list(pos_list.text) + if len(values) < 4 or len(values) % 2 != 0: + continue + for i in range(0, len(values), 2): + lat = values[i] + lon = values[i + 1] + coords.append((mercator_x(lon), mercator_y(lat))) + if len(coords) >= 2: + yield LineString(coords) + + +def mysql_connect(database: str | None = None): + kwargs = { + "host": DB_HOST, + "user": DB_USER, + "password": DB_PASSWORD, + "charset": "utf8mb4", + "autocommit": False, + "cursorclass": pymysql.cursors.Cursor, + } + if DB_SOCKET: + kwargs["unix_socket"] = DB_SOCKET + if database: + kwargs["database"] = database + return pymysql.connect(**kwargs) + + +def ensure_database() -> None: + conn = mysql_connect() + try: + with conn.cursor() as cur: + cur.execute( + f"CREATE DATABASE IF NOT EXISTS `{DB_NAME}` " + "DEFAULT CHARACTER SET utf8mb4 COLLATE utf8mb4_unicode_ci" + ) + conn.commit() + finally: + conn.close() + + +def ensure_schema(conn) -> None: + ddl = [ + "DROP TABLE IF EXISTS navsea_grid_cell", + "DROP TABLE IF EXISTS navsea_grid_package_stat", + "DROP TABLE IF EXISTS navsea_grid_layer_meta", + """ + CREATE TABLE navsea_grid_layer_meta ( + layer_name VARCHAR(32) NOT NULL PRIMARY KEY, + description VARCHAR(255) NOT NULL, + source_desc TEXT NOT NULL, + cell_size_m DOUBLE NOT NULL, + feature_count BIGINT NOT NULL, + bbox_min_lon DOUBLE NOT NULL, + bbox_min_lat DOUBLE NOT NULL, + bbox_max_lon DOUBLE NOT NULL, + bbox_max_lat DOUBLE NOT NULL, + export_file VARCHAR(255) DEFAULT NULL, + created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + """ + CREATE TABLE navsea_grid_package_stat ( + source_code VARCHAR(64) NOT NULL, + prefecture VARCHAR(128) NOT NULL, + zip_path TEXT NOT NULL, + line_count BIGINT NOT NULL, + point_count BIGINT NOT NULL, + bbox_min_lon DOUBLE NOT NULL, + bbox_min_lat DOUBLE NOT NULL, + bbox_max_lon DOUBLE NOT NULL, + bbox_max_lat DOUBLE NOT NULL, + grid_min_lon DOUBLE NOT NULL, + grid_min_lat DOUBLE NOT NULL, + grid_max_lon DOUBLE NOT NULL, + grid_max_lat DOUBLE NOT NULL, + blocked_count BIGINT NOT NULL, + candidate_count BIGINT NOT NULL + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + """ + CREATE TABLE navsea_grid_cell ( + id BIGINT UNSIGNED NOT NULL AUTO_INCREMENT PRIMARY KEY, + layer_name VARCHAR(32) NOT NULL, + cell_id VARCHAR(64) NOT NULL, + row_idx INT NOT NULL, + col_idx INT NOT NULL, + cell_size_m DOUBLE NOT NULL, + state_name VARCHAR(32) NOT NULL, + class_name VARCHAR(32) NOT NULL, + source_name VARCHAR(191) NOT NULL, + min_lon DOUBLE NOT NULL, + min_lat DOUBLE NOT NULL, + max_lon DOUBLE NOT NULL, + max_lat DOUBLE NOT NULL, + created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP, + UNIQUE KEY uniq_layer_cell (layer_name, cell_id), + KEY idx_layer_state (layer_name, state_name), + KEY idx_layer_rowcol (layer_name, row_idx, col_idx) + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + ] + with conn.cursor() as cur: + for stmt in ddl: + cur.execute(stmt) + conn.commit() + + +def insert_cells(cur, rows: list[tuple]) -> None: + cur.executemany( + """ + INSERT INTO navsea_grid_cell + (layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, source_name, + min_lon, min_lat, max_lon, max_lat) + VALUES + (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s) + ON DUPLICATE KEY UPDATE + row_idx = row_idx, + col_idx = col_idx, + cell_size_m = cell_size_m, + state_name = IF(VALUES(state_name) = 'HARD_BLOCKED' AND state_name <> 'HARD_BLOCKED', VALUES(state_name), state_name), + class_name = IF(VALUES(state_name) = 'HARD_BLOCKED' AND state_name <> 'HARD_BLOCKED', VALUES(class_name), class_name), + source_name = IF(VALUES(state_name) = 'HARD_BLOCKED' AND state_name <> 'HARD_BLOCKED', VALUES(source_name), source_name), + min_lon = min_lon, + min_lat = min_lat, + max_lon = max_lon, + max_lat = max_lat + """, + rows, + ) + + +def upsert_meta( + cur, + *, + layer_name: str, + description: str, + source_desc: str, + cell_size_m: float, + feature_count: int, + bbox: tuple[float, float, float, float], + export_file: str | None, +) -> None: + cur.execute( + """ + INSERT INTO navsea_grid_layer_meta + (layer_name, description, source_desc, cell_size_m, feature_count, + bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat, export_file) + VALUES + (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s) + ON DUPLICATE KEY UPDATE + description=VALUES(description), + source_desc=VALUES(source_desc), + cell_size_m=VALUES(cell_size_m), + feature_count=VALUES(feature_count), + bbox_min_lon=VALUES(bbox_min_lon), + bbox_min_lat=VALUES(bbox_min_lat), + bbox_max_lon=VALUES(bbox_max_lon), + bbox_max_lat=VALUES(bbox_max_lat), + export_file=VALUES(export_file) + """, + ( + layer_name, + description, + source_desc, + cell_size_m, + feature_count, + bbox[0], + bbox[1], + bbox[2], + bbox[3], + export_file, + ), + ) + + +def main() -> None: + parser = argparse.ArgumentParser(description="Build Japan coastline 200m grid directly into MySQL") + parser.add_argument( + "--input", + dest="inputs", + action="append", + default=[], + help="海岸线 GML zip。可重复指定；默认自动扫描 coastline/C23-06_*_GML.zip", + ) + parser.add_argument( + "--db-name", + default=DB_NAME, + help="MySQL 数据库名", + ) + parser.add_argument( + "--out-dir", + default=DEFAULT_OUT_DIR, + help="输出目录，用于可选 GeoJSON 和元数据", + ) + parser.add_argument("--buffer-m", type=float, default=50.0, help="海岸线缓冲距离，默认 50m") + parser.add_argument( + "--margin-cells", + type=int, + default=2, + help="输出范围外扩的格子数，默认 2 个格子", + ) + parser.add_argument( + "--write-geojson", + action="store_true", + help="同时输出 coastline.geojson 和 grid.geojson", + ) + args = parser.parse_args() + + project_root = Path(__file__).resolve().parent.parent + if args.inputs: + input_paths: list[Path] = [] + for item in args.inputs: + if any(ch in item for ch in "*?[]"): + input_paths.extend(Path(path) for path in sorted(glob.glob(item))) + else: + input_paths.append(Path(item)) + else: + input_paths = [Path(path) for path in sorted(glob.glob(str(project_root / DEFAULT_INPUT_GLOB)))] + + if not input_paths: + raise SystemExit("no coastline packages found") + + packages: list[SourcePackage] = [] + for path in input_paths: + if not path.is_absolute(): + path = project_root / path + if not path.exists(): + raise SystemExit(f"missing input: {path}") + packages.append(load_package(path)) + + out_dir = Path(args.out_dir) + if not out_dir.is_absolute(): + out_dir = project_root / out_dir + out_dir.mkdir(parents=True, exist_ok=True) + + ensure_database() + conn = mysql_connect(args.db_name) + try: + print(f"数据库 {args.db_name} 已连接，准备初始化表结构 ...") + ensure_schema(conn) + + total_lines = 0 + total_points = 0 + package_stats: list[dict] = [] + + for package in packages: + print(f"[{package.code}] 开始处理 {package.path.name}") + lines = list(iter_coastline_lines(package.path)) + if not lines: + print(f"[{package.code}] 没有可用海岸线，跳过") + continue + + coastline_geom = MultiLineString([list(line.coords) for line in lines]) + coastline_buffer = coastline_geom.buffer(args.buffer_m, cap_style=2, join_style=2) + prepared_buffer = prep(coastline_buffer) + coast_bounds = coastline_geom.bounds + + minx = align_floor(coast_bounds[0] - args.buffer_m - args.margin_cells * CELL_SIZE_M, CELL_SIZE_M) + miny = align_floor(coast_bounds[1] - args.buffer_m - args.margin_cells * CELL_SIZE_M, CELL_SIZE_M) + maxx = align_ceil(coast_bounds[2] + args.buffer_m + args.margin_cells * CELL_SIZE_M, CELL_SIZE_M) + maxy = align_ceil(coast_bounds[3] + args.buffer_m + args.margin_cells * CELL_SIZE_M, CELL_SIZE_M) + + line_count = len(lines) + point_count = sum(len(line.coords) for line in lines) + total_lines += line_count + total_points += point_count + + package_blocked = 0 + package_candidate = 0 + package_rows: list[tuple] = [] + + ix0 = int(round(minx / CELL_SIZE_M)) + iy0 = int(round(miny / CELL_SIZE_M)) + ix1 = int(round(maxx / CELL_SIZE_M)) + iy1 = int(round(maxy / CELL_SIZE_M)) + + for ix in range(ix0, ix1): + cell_minx = ix * CELL_SIZE_M + cell_maxx = cell_minx + CELL_SIZE_M + for iy in range(iy0, iy1): + cell_miny = iy * CELL_SIZE_M + cell_maxy = cell_miny + CELL_SIZE_M + cell = box(cell_minx, cell_miny, cell_maxx, cell_maxy) + intersects = prepared_buffer.intersects(cell) + if intersects: + cell_id = f"{ix}:{iy}" + package_rows.append( + ( + LAYER_NAME, + cell_id, + iy, + ix, + CELL_SIZE_M, + "HARD_BLOCKED", + "HARD_BLOCKED", + package.code, + cell_minx, + cell_miny, + cell_maxx, + cell_maxy, + ) + ) + package_blocked += 1 + else: + package_candidate += 1 + + if len(package_rows) >= 5000: + with conn.cursor() as cur: + insert_cells(cur, package_rows) + conn.commit() + print( + f"[{package.code}] 已写入 {len(package_rows)} 行，" + f"blocked={package_blocked} candidate={package_candidate}" + ) + package_rows.clear() + + if package_rows: + with conn.cursor() as cur: + insert_cells(cur, package_rows) + conn.commit() + + package_stats.append( + { + "source_code": package.code, + "prefecture": package.prefecture, + "zip_path": str(package.path), + "line_count": line_count, + "point_count": point_count, + "bbox_mercator": [coast_bounds[0], coast_bounds[1], coast_bounds[2], coast_bounds[3]], + "grid_bounds_mercator": [minx, miny, maxx, maxy], + "blocked_count": package_blocked, + "candidate_count": package_candidate, + } + ) + + with conn.cursor() as cur: + cur.execute( + """ + INSERT INTO navsea_grid_package_stat( + source_code, prefecture, zip_path, line_count, point_count, + bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat, + grid_min_lon, grid_min_lat, grid_max_lon, grid_max_lat, + blocked_count, candidate_count + ) VALUES ( + %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s + ) + """, + ( + package.code, + package.prefecture, + str(package.path), + line_count, + point_count, + coast_bounds[0], + coast_bounds[1], + coast_bounds[2], + coast_bounds[3], + minx, + miny, + maxx, + maxy, + package_blocked, + package_candidate, + ), + ) + conn.commit() + print( + f"[{package.code}] 完成 line={line_count} point={point_count} " + f"blocked={package_blocked} candidate={package_candidate}" + ) + + with conn.cursor() as cur: + cur.execute("ANALYZE TABLE navsea_grid_cell") + cur.execute("ANALYZE TABLE navsea_grid_package_stat") + cur.execute("ANALYZE TABLE navsea_grid_layer_meta") + conn.commit() + + with conn.cursor() as cur: + cur.execute(f"SELECT COUNT(*) FROM navsea_grid_cell WHERE layer_name=%s", (LAYER_NAME,)) + final_total = int(cur.fetchone()[0] or 0) + cur.execute( + f"SELECT COUNT(*) FROM navsea_grid_cell WHERE layer_name=%s AND state_name=%s", + (LAYER_NAME, "HARD_BLOCKED"), + ) + blocked_count = int(cur.fetchone()[0] or 0) + cur.execute( + f"SELECT COUNT(*) FROM navsea_grid_cell WHERE layer_name=%s AND state_name=%s", + (LAYER_NAME, "NAVIGABLE_CANDIDATE"), + ) + candidate_count = int(cur.fetchone()[0] or 0) + cur.execute( + """ + SELECT MIN(min_lon), MIN(min_lat), MAX(max_lon), MAX(max_lat) + FROM navsea_grid_cell + WHERE layer_name=%s + """, + (LAYER_NAME,), + ) + bbox_row = cur.fetchone() + + grid_bbox = ( + float(bbox_row[0]) if bbox_row and bbox_row[0] is not None else float("inf"), + float(bbox_row[1]) if bbox_row and bbox_row[1] is not None else float("inf"), + float(bbox_row[2]) if bbox_row and bbox_row[2] is not None else float("-inf"), + float(bbox_row[3]) if bbox_row and bbox_row[3] is not None else float("-inf"), + ) + + build_time = dt.datetime.now().isoformat(timespec="seconds") + grid_bbox_lonlat = mercator_bbox_to_lonlat(*grid_bbox) + + manifest = { + "build_time": build_time, + "db_name": args.db_name, + "source_scope": "japan", + "source_code": "C23-06_*", + "source_count": len(packages), + "source_packages": [str(pkg.path) for pkg in packages], + "coastline_line_count": total_lines, + "coastline_point_count": total_points, + "grid_cell_count": final_total, + "blocked_count": blocked_count, + "candidate_count": candidate_count, + "grid_bounds_lonlat": list(grid_bbox_lonlat), + "package_stats": package_stats, + "notes": "全国海岸线 200m 硬阻塞栅格，直接写入 MySQL", + } + manifest_path = out_dir / "japan_coast_grid_mysql.manifest.json" + manifest_path.write_text(json.dumps(manifest, ensure_ascii=False, indent=2), encoding="utf-8") + + with conn.cursor() as cur: + upsert_meta( + cur, + layer_name=LAYER_NAME, + description="全国海岸线 200m 硬阻塞格", + source_desc="coastline/C23-06_*_GML.zip", + cell_size_m=CELL_SIZE_M, + feature_count=final_total, + bbox=grid_bbox_lonlat, + export_file=str(manifest_path.relative_to(project_root)), + ) + conn.commit() + + print("完成") + print(f" MySQL: {args.db_name}") + print(f" packages={len(packages)}") + print(f" coastline_lines={total_lines}") + print(f" grid_cells={final_total}") + print(f" blocked={blocked_count}") + print(f" candidate={candidate_count}") + print(f" manifest: {manifest_path}") + finally: + conn.close() + + +if __name__ == "__main__": + main() diff --git a/coastline/build_japan_hazard_50m_mysql.py b/coastline/build_japan_hazard_50m_mysql.py new file mode 100644 index 0000000..0f24950 --- /dev/null +++ b/coastline/build_japan_hazard_50m_mysql.py @@ -0,0 +1,516 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import datetime as dt +import json +import math +from pathlib import Path +from typing import Iterable + +import mapbox_vector_tile +import pymysql +from shapely.geometry import box, shape +from shapely.prepared import prep + + +DB_NAME = "navsea_japan_coast_grid" +DB_USER = "root" +DB_PASSWORD = "2chi9ks2" +DB_HOST = "localhost" +DB_SOCKET = "/tmp/mysql.sock" + +DEFAULT_TILE_ROOT = Path("/home/wwwroot/pbf-delivery-full-20260418-rebuild") +DEFAULT_TILE_Z = 12 +DEFAULT_OUT_DIR = "src/pbf/coastline-mysql/japan_national" +DEFAULT_SOURCE_DESC = "pbf-delivery-full-20260418-rebuild z12" + +HAZARD_CELL_M = 50.0 +HAZARD_LAYER_NAME = "hazard_50m" +HAZARD_DESCRIPTION = "全国 PBF 海上障碍 50m 黄格" +HAZARD_LAYERS = ( + "navigation_hazard_area", + "fixed_fishing_gear_area", + "anchor_caution_hazard_area", + "navigation_hazard_point", + "anchor_caution_hazard_point", + "navigation_marks", +) +BREAKWATER_LAYER = "baseline_area" +POINT_LAYERS = { + "navigation_hazard_point", + "anchor_caution_hazard_point", + "navigation_marks", +} +HAZARD_CANONICAL_OBJECT_TYPES = { + "breakwater", +} + +RADIUS = 6378137.0 +MAX_MERCATOR_LAT = 85.0511287798066 + + +def mercator_x(lon: float) -> float: + return RADIUS * math.radians(lon) + + +def mercator_y(lat: float) -> float: + lat = max(min(lat, MAX_MERCATOR_LAT), -MAX_MERCATOR_LAT) + return RADIUS * math.log(math.tan(math.pi / 4.0 + math.radians(lat) / 2.0)) + + +def lon_from_mercator(x: float) -> float: + return math.degrees(x / RADIUS) + + +def lat_from_mercator(y: float) -> float: + return math.degrees(2.0 * math.atan(math.exp(y / RADIUS)) - math.pi / 2.0) + + +def rect_geojson(min_lon: float, min_lat: float, max_lon: float, max_lat: float) -> dict: + return { + "type": "Feature", + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [min_lon, min_lat], + [max_lon, min_lat], + [max_lon, max_lat], + [min_lon, max_lat], + [min_lon, min_lat], + ] + ], + }, + } + + +def local_name(tag: str) -> str: + return tag.split("}", 1)[-1] + + +def load_json(path: Path) -> dict: + return json.loads(path.read_text(encoding="utf-8")) + + +def stream_geojson(path: Path, features: Iterable[dict]) -> int: + path.parent.mkdir(parents=True, exist_ok=True) + count = 0 + with path.open("w", encoding="utf-8") as fh: + fh.write('{"type":"FeatureCollection","features":[\n') + first = True + for feature in features: + if not first: + fh.write(",\n") + fh.write(json.dumps(feature, ensure_ascii=False)) + first = False + count += 1 + fh.write("\n]}\n") + return count + + +def mysql_connect(database: str | None = None): + kwargs = { + "host": DB_HOST, + "user": DB_USER, + "password": DB_PASSWORD, + "charset": "utf8mb4", + "autocommit": False, + "cursorclass": pymysql.cursors.Cursor, + } + if DB_SOCKET: + kwargs["unix_socket"] = DB_SOCKET + if database: + kwargs["database"] = database + return pymysql.connect(**kwargs) + + +def ensure_database(name: str) -> None: + conn = mysql_connect() + try: + with conn.cursor() as cur: + cur.execute( + f"CREATE DATABASE IF NOT EXISTS `{name}` " + "DEFAULT CHARACTER SET utf8mb4 COLLATE utf8mb4_unicode_ci" + ) + conn.commit() + finally: + conn.close() + + +def ensure_schema(conn) -> None: + stmts = [ + """ + CREATE TABLE IF NOT EXISTS navsea_grid_layer_meta ( + layer_name VARCHAR(32) NOT NULL PRIMARY KEY, + description VARCHAR(255) NOT NULL, + source_desc TEXT NOT NULL, + cell_size_m DOUBLE NOT NULL, + feature_count BIGINT NOT NULL, + bbox_min_lon DOUBLE NOT NULL, + bbox_min_lat DOUBLE NOT NULL, + bbox_max_lon DOUBLE NOT NULL, + bbox_max_lat DOUBLE NOT NULL, + export_file VARCHAR(255) DEFAULT NULL, + created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + """ + CREATE TABLE IF NOT EXISTS navsea_grid_cell ( + id BIGINT UNSIGNED NOT NULL AUTO_INCREMENT PRIMARY KEY, + layer_name VARCHAR(32) NOT NULL, + cell_id VARCHAR(64) NOT NULL, + row_idx INT NOT NULL, + col_idx INT NOT NULL, + cell_size_m DOUBLE NOT NULL, + state_name VARCHAR(32) NOT NULL, + class_name VARCHAR(32) NOT NULL, + source_name VARCHAR(191) NOT NULL, + min_lon DOUBLE NOT NULL, + min_lat DOUBLE NOT NULL, + max_lon DOUBLE NOT NULL, + max_lat DOUBLE NOT NULL, + created_at TIMESTAMP NOT NULL DEFAULT CURRENT_TIMESTAMP, + UNIQUE KEY uniq_layer_cell (layer_name, cell_id), + KEY idx_layer_state (layer_name, state_name), + KEY idx_layer_rowcol (layer_name, row_idx, col_idx) + ) ENGINE=InnoDB DEFAULT CHARSET=utf8mb4 + """, + ] + with conn.cursor() as cur: + for stmt in stmts: + cur.execute(stmt) + conn.commit() + + +def replace_layer_rows(conn, layer_name: str) -> None: + with conn.cursor() as cur: + cur.execute("DELETE FROM navsea_grid_cell WHERE layer_name=%s", (layer_name,)) + cur.execute("DELETE FROM navsea_grid_layer_meta WHERE layer_name=%s", (layer_name,)) + conn.commit() + + +def insert_cells(cur, rows: list[tuple]) -> None: + cur.executemany( + """ + INSERT INTO navsea_grid_cell + (layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, source_name, + min_lon, min_lat, max_lon, max_lat) + VALUES + (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s) + """, + rows, + ) + + +def upsert_meta( + cur, + *, + layer_name: str, + description: str, + source_desc: str, + cell_size_m: float, + feature_count: int, + bbox: tuple[float, float, float, float], + export_file: str, +) -> None: + cur.execute( + """ + INSERT INTO navsea_grid_layer_meta + (layer_name, description, source_desc, cell_size_m, feature_count, + bbox_min_lon, bbox_min_lat, bbox_max_lon, bbox_max_lat, export_file) + VALUES + (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s) + ON DUPLICATE KEY UPDATE + description=VALUES(description), + source_desc=VALUES(source_desc), + cell_size_m=VALUES(cell_size_m), + feature_count=VALUES(feature_count), + bbox_min_lon=VALUES(bbox_min_lon), + bbox_min_lat=VALUES(bbox_min_lat), + bbox_max_lon=VALUES(bbox_max_lon), + bbox_max_lat=VALUES(bbox_max_lat), + export_file=VALUES(export_file) + """, + ( + layer_name, + description, + source_desc, + cell_size_m, + feature_count, + bbox[0], + bbox[1], + bbox[2], + bbox[3], + export_file, + ), + ) + + +def tile_mercator_bounds(z: int, x: int, y: int) -> tuple[float, float, float, float]: + n = 2**z + lon_left = x / n * 360.0 - 180.0 + lon_right = (x + 1) / n * 360.0 - 180.0 + lat_top = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * y / n)))) + lat_bottom = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * (y + 1) / n)))) + return mercator_x(lon_left), mercator_y(lat_bottom), mercator_x(lon_right), mercator_y(lat_top) + + +def transform_geom_from_tile(geom, bounds_m: tuple[float, float, float, float], extent: int): + minx, miny, maxx, maxy = bounds_m + dx = maxx - minx + dy = maxy - miny + + def coord(x: float, y: float): + mx = minx + (x / extent) * dx + my = miny + (y / extent) * dy + return mx, my + + def walk(obj): + if isinstance(obj[0], (int, float)): + return coord(obj[0], obj[1]) + return [walk(item) for item in obj] + + return walk(geom) + + +def iter_source_tiles(tile_root: Path, zoom: int) -> list[Path]: + zoom_root = tile_root / str(zoom) + if not zoom_root.exists(): + raise SystemExit(f"missing tile zoom root: {zoom_root}") + return sorted(zoom_root.glob("*/*.pbf")) + + +def iter_hazard_cells(tile_root: Path, zoom: int) -> Iterable[tuple[int, int]]: + hazard_cells: set[tuple[int, int]] = set() + tiles = iter_source_tiles(tile_root, zoom) + for index, p in enumerate(tiles, start=1): + try: + tx = int(p.parent.name) + ty = int(p.stem) + except ValueError: + continue + try: + tile = mapbox_vector_tile.decode(p.read_bytes()) + except Exception as exc: + print(f"[hazard_50m] 跳过无法解析的 tile: {p} ({exc.__class__.__name__})") + continue + + if index % 5000 == 0: + print(f"[hazard_50m] 已扫描 tile {index}/{len(tiles)}: {p.parent.parent.name}/{tx}/{ty}.pbf") + + for layer_name in HAZARD_LAYERS + (BREAKWATER_LAYER,): + layer = tile.get(layer_name) + if not layer: + continue + extent = int(layer.get("extent") or 1048576) + minx, miny, maxx, maxy = tile_mercator_bounds(zoom, tx, ty) + for feature in layer.get("features", []): + props = feature.get("properties") or {} + if layer_name == BREAKWATER_LAYER and props.get("canonical_object_type") not in HAZARD_CANONICAL_OBJECT_TYPES: + continue + if ( + props.get("canonical_object_type") == "fish_reef" + or props.get("chart_symbol_code") == "fish_reef" + or props.get("class_name") == "魚礁" + ): + continue + geom = feature.get("geometry") or {} + coords = geom.get("coordinates") + if not coords: + continue + shp_coords = transform_geom_from_tile(coords, (minx, miny, maxx, maxy), extent) + shp = shape({"type": geom.get("type"), "coordinates": shp_coords}) + if shp.is_empty: + continue + if layer_name in POINT_LAYERS: + shp = shp.buffer(25.0) + minx2, miny2, maxx2, maxy2 = shp.bounds + start_x = math.floor(minx2 / HAZARD_CELL_M) * HAZARD_CELL_M + start_y = math.floor(miny2 / HAZARD_CELL_M) * HAZARD_CELL_M + end_x = math.ceil(maxx2 / HAZARD_CELL_M) * HAZARD_CELL_M + end_y = math.ceil(maxy2 / HAZARD_CELL_M) * HAZARD_CELL_M + prep_geom = prep(shp) + cy = start_y + while cy < end_y: + cx = start_x + while cx < end_x: + cell = box(cx, cy, cx + HAZARD_CELL_M, cy + HAZARD_CELL_M) + if prep_geom.intersects(cell): + hazard_cells.add((int(round(cx / HAZARD_CELL_M)), int(round(cy / HAZARD_CELL_M)))) + cx += HAZARD_CELL_M + cy += HAZARD_CELL_M + return sorted(hazard_cells, key=lambda item: (item[1], item[0])) + + +def cell_bbox_from_mercator(ix: int, iy: int, cell_m: float) -> tuple[float, float, float, float]: + minx = ix * cell_m + miny = iy * cell_m + maxx = minx + cell_m + maxy = miny + cell_m + return ( + lon_from_mercator(minx), + lat_from_mercator(miny), + lon_from_mercator(maxx), + lat_from_mercator(maxy), + ) + + +def import_hazard_layer( + conn, + *, + layer_name: str, + description: str, + source_desc: str, + cell_size_m: float, + tile_root: Path, + zoom: int, + export_path: Path, +) -> tuple[int, tuple[float, float, float, float]]: + cells = iter_hazard_cells(tile_root, zoom) + batch: list[tuple] = [] + count = 0 + export_bbox = [float("inf"), float("inf"), float("-inf"), float("-inf")] + + with conn.cursor() as cur: + print(f"[{layer_name}] 开始扫描全国 PBF 危险层：{tile_root} / z{zoom}") + with export_path.open("w", encoding="utf-8") as fh: + fh.write('{"type":"FeatureCollection","features":[\n') + first = True + for ix, iy in cells: + min_lon, min_lat, max_lon, max_lat = cell_bbox_from_mercator(ix, iy, cell_size_m) + export_bbox[0] = min(export_bbox[0], min_lon) + export_bbox[1] = min(export_bbox[1], min_lat) + export_bbox[2] = max(export_bbox[2], max_lon) + export_bbox[3] = max(export_bbox[3], max_lat) + cell_id = f"{ix}:{iy}" + batch.append( + ( + layer_name, + cell_id, + iy, + ix, + cell_size_m, + "HAZARD_50M", + "HAZARD", + source_desc, + min_lon, + min_lat, + max_lon, + max_lat, + ) + ) + feature = { + "type": "Feature", + "properties": { + "layer_name": layer_name, + "cell_id": cell_id, + "row": iy, + "col": ix, + "cell_size_m": cell_size_m, + "state_name": "HAZARD_50M", + "class_name": "HAZARD", + "source_name": source_desc, + }, + "geometry": rect_geojson(min_lon, min_lat, max_lon, max_lat)["geometry"], + } + if not first: + fh.write(",\n") + fh.write(json.dumps(feature, ensure_ascii=False)) + first = False + count += 1 + if len(batch) >= 5000: + insert_cells(cur, batch) + conn.commit() + batch.clear() + if count % 20000 == 0: + print(f"[{layer_name}] 已处理 {count} 个格子 ...") + if batch: + insert_cells(cur, batch) + conn.commit() + fh.write("\n]}\n") + + if count == 0: + raise RuntimeError(f"{layer_name} 扫描结果为空，请检查 tile_root={tile_root} zoom={zoom}") + + print(f"[{layer_name}] 导入完成：{count} 条，导出 {export_path}") + return count, (export_bbox[0], export_bbox[1], export_bbox[2], export_bbox[3]) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description="重算全国 50x50 海上障碍格并写入统一 MySQL 库。") + parser.add_argument("--db-name", default=DB_NAME, help="MySQL 数据库名") + parser.add_argument("--tile-root", type=Path, default=DEFAULT_TILE_ROOT, help="全国 PBF 根目录") + parser.add_argument("--zoom", type=int, default=DEFAULT_TILE_Z, help="PBF 瓦片 zoom，默认 12") + parser.add_argument("--out-dir", default=DEFAULT_OUT_DIR, help="GeoJSON 输出目录") + parser.add_argument("--source-desc", default=DEFAULT_SOURCE_DESC, help="写入元数据的来源说明") + return parser.parse_args() + + +def main() -> None: + args = parse_args() + + project_root = Path(__file__).resolve().parent.parent + out_dir = Path(args.out_dir) + if not out_dir.is_absolute(): + out_dir = project_root / out_dir + out_dir.mkdir(parents=True, exist_ok=True) + + ensure_database(args.db_name) + conn = mysql_connect(args.db_name) + try: + print(f"数据库 {args.db_name} 已连接，准备初始化表结构 ...") + ensure_schema(conn) + replace_layer_rows(conn, HAZARD_LAYER_NAME) + + hazard_export = out_dir / "hazard_50m_grid.geojson" + hazard_count, hazard_bbox = import_hazard_layer( + conn, + layer_name=HAZARD_LAYER_NAME, + description=HAZARD_DESCRIPTION, + source_desc=args.source_desc, + cell_size_m=HAZARD_CELL_M, + tile_root=args.tile_root, + zoom=args.zoom, + export_path=hazard_export, + ) + + with conn.cursor() as cur: + upsert_meta( + cur, + layer_name=HAZARD_LAYER_NAME, + description=HAZARD_DESCRIPTION, + source_desc=args.source_desc, + cell_size_m=HAZARD_CELL_M, + feature_count=hazard_count, + bbox=hazard_bbox, + export_file=str(hazard_export.relative_to(project_root)), + ) + conn.commit() + + manifest = { + "database": args.db_name, + "generated_at": dt.datetime.now().isoformat(timespec="seconds"), + "layer": { + "name": HAZARD_LAYER_NAME, + "count": hazard_count, + "bbox_lonlat": list(hazard_bbox), + "cell_size_m": HAZARD_CELL_M, + "export_file": str(hazard_export.relative_to(project_root)), + "source_desc": args.source_desc, + "tile_root": str(args.tile_root), + "zoom": args.zoom, + }, + } + manifest_path = out_dir / "manifest.json" + manifest_path.write_text(json.dumps(manifest, ensure_ascii=False, indent=2), encoding="utf-8") + + print("已导入 MySQL 数据库:", args.db_name) + print(" - hazard_50m:", hazard_count, hazard_export) + print(" - manifest:", manifest_path) + finally: + conn.close() + + +if __name__ == "__main__": + main() diff --git a/coastline/export_navgrid_mysql_assets.py b/coastline/export_navgrid_mysql_assets.py new file mode 100644 index 0000000..2cbcf1c --- /dev/null +++ b/coastline/export_navgrid_mysql_assets.py @@ -0,0 +1,485 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import datetime as dt +import json +import math +from pathlib import Path + +import pymysql + + +DB_NAME = "navsea_japan_coast_grid" +DB_USER = "root" +DB_PASSWORD = "2chi9ks2" +DB_HOST = "localhost" +DB_SOCKET = "/tmp/mysql.sock" + +DEFAULT_OUT_DIR = "src/pbf/coastline-mysql/japan_national" + +DENSITY_OVERVIEW_FILE = "density_overview_grid.geojson" +CELL_SIZES = { + "coast_200m": 200.0, + "fish_port_20m": 20.0, + "hazard_50m": 50.0, +} + +RADIUS = 6378137.0 +MAX_MERCATOR_LAT = 85.0511287798066 + + +def mysql_connect(database: str, *, cursorclass=pymysql.cursors.Cursor): + kwargs = { + "host": DB_HOST, + "user": DB_USER, + "password": DB_PASSWORD, + "database": database, + "charset": "utf8mb4", + "cursorclass": cursorclass, + } + if DB_SOCKET: + kwargs["unix_socket"] = DB_SOCKET + return pymysql.connect(**kwargs) + + +def rect_geometry(min_lon: float, min_lat: float, max_lon: float, max_lat: float) -> dict: + return { + "type": "Polygon", + "coordinates": [ + [ + [min_lon, min_lat], + [max_lon, min_lat], + [max_lon, max_lat], + [min_lon, max_lat], + [min_lon, min_lat], + ] + ], + } + + +def lon_from_mercator(x: float) -> float: + return math.degrees(x / RADIUS) + + +def lat_from_mercator(y: float) -> float: + return math.degrees(2.0 * math.atan(math.exp(y / RADIUS)) - math.pi / 2.0) + + +def mercator_bbox_to_lonlat(minx: float, miny: float, maxx: float, maxy: float) -> tuple[float, float, float, float]: + return ( + lon_from_mercator(minx), + lat_from_mercator(miny), + lon_from_mercator(maxx), + lat_from_mercator(maxy), + ) + + +def normalize_bbox(layer_name: str, min_lon: float, min_lat: float, max_lon: float, max_lat: float) -> tuple[float, float, float, float]: + if layer_name == "coast_200m" and ( + abs(min_lon) > 180.0 or abs(max_lon) > 180.0 or abs(min_lat) > 90.0 or abs(max_lat) > 90.0 + ): + return mercator_bbox_to_lonlat(min_lon, min_lat, max_lon, max_lat) + return min_lon, min_lat, max_lon, max_lat + + +def fetch_stats(conn, where_sql: str, params: tuple) -> tuple[int, list[float] | None]: + sql = f""" + SELECT COUNT(*), MIN(min_lon), MIN(min_lat), MAX(max_lon), MAX(max_lat) + FROM navsea_grid_cell + WHERE {where_sql} + """ + with conn.cursor() as cur: + cur.execute(sql, params) + row = cur.fetchone() + count = int(row[0] or 0) + if count == 0: + return 0, None + min_lon, min_lat, max_lon, max_lat = normalize_bbox( + params[0] if params else "", + float(row[1]), + float(row[2]), + float(row[3]), + float(row[4]), + ) + return count, [min_lon, min_lat, max_lon, max_lat] + + +def export_layer( + conn, + *, + path: Path, + where_sql: str, + params: tuple, +) -> int: + path.parent.mkdir(parents=True, exist_ok=True) + sql = f""" + SELECT layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, + source_name, min_lon, min_lat, max_lon, max_lat + FROM navsea_grid_cell + WHERE {where_sql} + ORDER BY row_idx, col_idx + """ + count = 0 + with conn.cursor(pymysql.cursors.SSCursor) as cur, path.open("w", encoding="utf-8") as fh: + cur.execute(sql, params) + fh.write('{"type":"FeatureCollection","features":[\n') + first = True + for row in cur: + ( + layer_name, + cell_id, + row_idx, + col_idx, + cell_size_m, + state_name, + class_name, + source_name, + min_lon, + min_lat, + max_lon, + max_lat, + ) = row + min_lon, min_lat, max_lon, max_lat = normalize_bbox( + row[0], + float(min_lon), + float(min_lat), + float(max_lon), + float(max_lat), + ) + feature = { + "type": "Feature", + "properties": { + "layer_name": layer_name, + "cell_id": cell_id, + "row": int(row_idx), + "col": int(col_idx), + "cell_size_m": float(cell_size_m), + "state_name": state_name, + "class_name": class_name, + "source_name": source_name, + }, + "geometry": rect_geometry(min_lon, min_lat, max_lon, max_lat), + } + if not first: + fh.write(",\n") + fh.write(json.dumps(feature, ensure_ascii=False)) + first = False + count += 1 + fh.write("\n]}\n") + return count + + +def _range_for_overlap(min_value: float, max_value: float, fine_size_m: float) -> tuple[int, int]: + start = int(math.floor(min_value / fine_size_m)) + end = int(math.floor((max_value - 1e-9) / fine_size_m)) + return start, end + + +def _overlaps_finer_cells(row_idx: int, col_idx: int, coarse_size_m: float, fine_size_m: float, fine_keys: set[tuple[int, int]]) -> bool: + min_x = col_idx * coarse_size_m + min_y = row_idx * coarse_size_m + max_x = min_x + coarse_size_m + max_y = min_y + coarse_size_m + row_start, row_end = _range_for_overlap(min_y, max_y, fine_size_m) + col_start, col_end = _range_for_overlap(min_x, max_x, fine_size_m) + for fine_row in range(row_start, row_end + 1): + for fine_col in range(col_start, col_end + 1): + if (fine_row, fine_col) in fine_keys: + return True + return False + + +def export_density_overview(conn, *, path: Path) -> tuple[int, dict]: + path.parent.mkdir(parents=True, exist_ok=True) + + kept_counts = { + "20": 0, + "50": 0, + "200": 0, + } + kept_bbox = [float("inf"), float("inf"), float("-inf"), float("-inf")] + keys_20: set[tuple[int, int]] = set() + keys_50: set[tuple[int, int]] = set() + + def write_feature(fh, feature: dict, first: bool) -> bool: + if not first: + fh.write(",\n") + fh.write(json.dumps(feature, ensure_ascii=False)) + return False + + with conn.cursor(pymysql.cursors.SSCursor) as cur, path.open("w", encoding="utf-8") as fh: + fh.write('{"type":"FeatureCollection","features":[\n') + first = True + + # 20m: always keep + cur.execute( + """ + SELECT layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, + source_name, min_lon, min_lat, max_lon, max_lat + FROM navsea_grid_cell + WHERE layer_name=%s + ORDER BY row_idx, col_idx + """, + ("fish_port_20m",), + ) + for row in cur: + ( + layer_name, + cell_id, + row_idx, + col_idx, + cell_size_m, + state_name, + class_name, + source_name, + min_lon, + min_lat, + max_lon, + max_lat, + ) = row + keys_20.add((int(row_idx), int(col_idx))) + kept_counts["20"] += 1 + min_lon, min_lat, max_lon, max_lat = normalize_bbox( + layer_name, + float(min_lon), + float(min_lat), + float(max_lon), + float(max_lat), + ) + kept_bbox[0] = min(kept_bbox[0], min_lon) + kept_bbox[1] = min(kept_bbox[1], min_lat) + kept_bbox[2] = max(kept_bbox[2], max_lon) + kept_bbox[3] = max(kept_bbox[3], max_lat) + feature = { + "type": "Feature", + "properties": { + "layer_name": layer_name, + "cell_id": cell_id, + "row": int(row_idx), + "col": int(col_idx), + "cell_size_m": float(cell_size_m), + "state_name": state_name, + "class_name": class_name, + "source_name": source_name, + "density_key": "20", + "density_level": 3, + "density_name": "高密度", + }, + "geometry": rect_geometry(min_lon, min_lat, max_lon, max_lat), + } + first = write_feature(fh, feature, first) + + # 50m: keep only when no 20m cell overlaps this area + cur.execute( + """ + SELECT layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, + source_name, min_lon, min_lat, max_lon, max_lat + FROM navsea_grid_cell + WHERE layer_name=%s + ORDER BY row_idx, col_idx + """, + ("hazard_50m",), + ) + for row in cur: + ( + layer_name, + cell_id, + row_idx, + col_idx, + cell_size_m, + state_name, + class_name, + source_name, + min_lon, + min_lat, + max_lon, + max_lat, + ) = row + row_idx_i = int(row_idx) + col_idx_i = int(col_idx) + if _overlaps_finer_cells(row_idx_i, col_idx_i, CELL_SIZES["hazard_50m"], CELL_SIZES["fish_port_20m"], keys_20): + continue + keys_50.add((row_idx_i, col_idx_i)) + kept_counts["50"] += 1 + min_lon, min_lat, max_lon, max_lat = normalize_bbox( + layer_name, + float(min_lon), + float(min_lat), + float(max_lon), + float(max_lat), + ) + kept_bbox[0] = min(kept_bbox[0], min_lon) + kept_bbox[1] = min(kept_bbox[1], min_lat) + kept_bbox[2] = max(kept_bbox[2], max_lon) + kept_bbox[3] = max(kept_bbox[3], max_lat) + feature = { + "type": "Feature", + "properties": { + "layer_name": layer_name, + "cell_id": cell_id, + "row": row_idx_i, + "col": col_idx_i, + "cell_size_m": float(cell_size_m), + "state_name": state_name, + "class_name": class_name, + "source_name": source_name, + "density_key": "50", + "density_level": 2, + "density_name": "中密度", + }, + "geometry": rect_geometry(min_lon, min_lat, max_lon, max_lat), + } + first = write_feature(fh, feature, first) + + # 200m: keep only when neither 20m nor 50m cell overlaps this area + cur.execute( + """ + SELECT layer_name, cell_id, row_idx, col_idx, cell_size_m, state_name, class_name, + source_name, min_lon, min_lat, max_lon, max_lat + FROM navsea_grid_cell + WHERE layer_name=%s + ORDER BY row_idx, col_idx + """, + ("coast_200m",), + ) + for row in cur: + ( + layer_name, + cell_id, + row_idx, + col_idx, + cell_size_m, + state_name, + class_name, + source_name, + min_lon, + min_lat, + max_lon, + max_lat, + ) = row + row_idx_i = int(row_idx) + col_idx_i = int(col_idx) + if _overlaps_finer_cells(row_idx_i, col_idx_i, CELL_SIZES["coast_200m"], CELL_SIZES["fish_port_20m"], keys_20): + continue + if _overlaps_finer_cells(row_idx_i, col_idx_i, CELL_SIZES["coast_200m"], CELL_SIZES["hazard_50m"], keys_50): + continue + kept_counts["200"] += 1 + min_lon, min_lat, max_lon, max_lat = normalize_bbox( + layer_name, + float(min_lon), + float(min_lat), + float(max_lon), + float(max_lat), + ) + kept_bbox[0] = min(kept_bbox[0], min_lon) + kept_bbox[1] = min(kept_bbox[1], min_lat) + kept_bbox[2] = max(kept_bbox[2], max_lon) + kept_bbox[3] = max(kept_bbox[3], max_lat) + feature = { + "type": "Feature", + "properties": { + "layer_name": layer_name, + "cell_id": cell_id, + "row": row_idx_i, + "col": col_idx_i, + "cell_size_m": float(cell_size_m), + "state_name": state_name, + "class_name": class_name, + "source_name": source_name, + "density_key": "200", + "density_level": 1, + "density_name": "低密度", + }, + "geometry": rect_geometry(min_lon, min_lat, max_lon, max_lat), + } + first = write_feature(fh, feature, first) + + fh.write("\n]}\n") + + total = kept_counts["20"] + kept_counts["50"] + kept_counts["200"] + if total == 0: + raise RuntimeError("density overview export produced no features") + summary = { + "count": total, + "by_density": kept_counts, + "bbox": None if kept_bbox[0] == float("inf") else kept_bbox, + } + return total, summary + + +def main() -> None: + parser = argparse.ArgumentParser(description="Export NavSea MySQL grid layers to static GeoJSON assets") + parser.add_argument("--db-name", default=DB_NAME) + parser.add_argument("--out-dir", default=DEFAULT_OUT_DIR) + args = parser.parse_args() + + project_root = Path(__file__).resolve().parent.parent + out_dir = Path(args.out_dir) + if not out_dir.is_absolute(): + out_dir = project_root / out_dir + out_dir.mkdir(parents=True, exist_ok=True) + + layers = { + "coast_200m": { + "file": "coast_200m_grid.geojson", + "where": "layer_name=%s", + "params": ("coast_200m",), + "cell_size_m": 200.0, + }, + "fish_port_20m": { + "file": "fish_port_20m_grid.geojson", + "where": "layer_name=%s", + "params": ("fish_port_20m",), + "cell_size_m": 20.0, + }, + "hazard_50m": { + "file": "hazard_50m_grid.geojson", + "where": "layer_name=%s", + "params": ("hazard_50m",), + "cell_size_m": 50.0, + }, + } + + manifest_layers: dict[str, dict] = {} + conn = mysql_connect(args.db_name) + try: + for layer_name, spec in layers.items(): + out_path = out_dir / spec["file"] + count = export_layer(conn, path=out_path, where_sql=spec["where"], params=spec["params"]) + stat_count, bbox = fetch_stats(conn, spec["where"], spec["params"]) + if count != stat_count: + raise RuntimeError(f"{layer_name} export count mismatch: {count} != {stat_count}") + manifest_layers[layer_name] = { + "count": count, + "bbox_lonlat": bbox, + "cell_size_m": spec["cell_size_m"], + "export_file": str(out_path.relative_to(project_root)), + } + print(f"{layer_name}: {count} -> {out_path}") + + density_path = out_dir / DENSITY_OVERVIEW_FILE + density_count, density_summary = export_density_overview(conn, path=density_path) + manifest_layers["density_overview"] = { + "count": density_count, + "bbox_lonlat": density_summary["bbox"], + "cell_size_m": None, + "export_file": str(density_path.relative_to(project_root)), + "by_density": density_summary["by_density"], + } + print(f"density_overview: {density_count} -> {density_path}") + finally: + conn.close() + + manifest = { + "database": args.db_name, + "generated_at": dt.datetime.now().isoformat(timespec="seconds"), + "layers": manifest_layers, + } + manifest_path = out_dir / "manifest.json" + manifest_path.write_text(json.dumps(manifest, ensure_ascii=False, indent=2), encoding="utf-8") + print(f"manifest: {manifest_path}") + + +if __name__ == "__main__": + main() diff --git a/find_fpc_by_port_name.py b/find_fpc_by_port_name.py new file mode 100644 index 0000000..d47d8c1 --- /dev/null +++ b/find_fpc_by_port_name.py @@ -0,0 +1,282 @@ +#!/usr/bin/env python3 +from __future__ import annotations + +import argparse +import csv +import json +import re +import zipfile +import xml.etree.ElementTree as ET +from dataclasses import dataclass, asdict +from pathlib import Path + + +FISH_SRC_DEFAULT = "coastline/C09-06.zip" + + +@dataclass(frozen=True) +class MatchRow: + fpc: str + prc: str | None + primary_name: str | None + secondary_name: str | None + source_id: str | None + score: tuple[int, int, str] + + +def local_name(tag: str) -> str: + return tag.split("}", 1)[-1] + + +def normalize_text(value: str) -> str: + value = value.strip() + value = value.replace("　", " ") + value = re.sub(r"\s+", "", value) + return value.lower() + + +def load_xml(zip_path: Path) -> ET.Element: + with zipfile.ZipFile(zip_path) as zf: + xml_name = next(name for name in zf.namelist() if name.endswith(".xml") and "META" not in name) + return ET.fromstring(zf.read(xml_name)) + + +def get_text(feature: ET.Element, key: str) -> str | None: + el = feature.find(f"{{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}}{key}") + if el is None or not (el.text or "").strip(): + return None + return el.text.strip() + + +def feature_name_tokens(feature: ET.Element) -> list[str]: + tokens: list[str] = [] + for key in ("NA2", "NA4", "FCF", "AAC", "CFP", "FPA"): + value = get_text(feature, key) + if value: + tokens.append(value) + return tokens + + +def feature_display_name(feature: ET.Element) -> tuple[str | None, str | None]: + na2 = get_text(feature, "NA2") + na4 = get_text(feature, "NA4") + fcf = get_text(feature, "FCF") + primary = na2 or na4 or fcf + secondary = None + if primary == na2: + secondary = fcf or na4 + elif primary == na4: + secondary = na2 or fcf + else: + secondary = na2 or na4 + return primary, secondary + + +def catalog_display_name(row: dict[str, str]) -> str | None: + for key in ("name", "name_ja", "port_name", "港名", "漁港名", "名称"): + value = row.get(key) + if value and value.strip(): + return value.strip() + return None + + +def catalog_fpc(row: dict[str, str]) -> str | None: + for key in ("fpc", "FPC"): + value = row.get(key) + if value and value.strip(): + return value.strip() + return None + + +def catalog_prc(row: dict[str, str]) -> str | None: + for key in ("prc", "PRC"): + value = row.get(key) + if value and value.strip(): + return value.strip() + return None + + +def score_match(query: str, tokens: list[str], primary: str | None) -> tuple[int, int, str] | None: + q = normalize_text(query) + token_norms = [normalize_text(token) for token in tokens if token] + primary_norm = normalize_text(primary) if primary else "" + + if primary_norm == q: + return (0, 0, primary_norm) + if q in token_norms: + return (1, 0, primary_norm) + if primary_norm and q in primary_norm: + return (2, len(primary_norm), primary_norm) + for token in token_norms: + if q in token: + return (3, len(token), primary_norm or token) + return None + + +def load_catalog_matches(catalog_path: Path, query: str) -> list[MatchRow]: + suffix = catalog_path.suffix.lower() + rows: list[dict[str, str]] = [] + if suffix in (".csv", ".tsv"): + delimiter = "\t" if suffix == ".tsv" else "," + with catalog_path.open("r", encoding="utf-8-sig", newline="") as fh: + reader = csv.DictReader(fh, delimiter=delimiter) + rows.extend(dict(row) for row in reader) + elif suffix in (".json", ".jsonl"): + text = catalog_path.read_text(encoding="utf-8") + if suffix == ".jsonl": + rows.extend(json.loads(line) for line in text.splitlines() if line.strip()) + else: + payload = json.loads(text) + if isinstance(payload, list): + rows.extend(dict(row) for row in payload) + elif isinstance(payload, dict) and isinstance(payload.get("items"), list): + rows.extend(dict(row) for row in payload["items"]) + else: + raise SystemExit("catalog json must be a list or {items:[...]}") + else: + raise SystemExit("catalog must be .csv, .tsv, .json or .jsonl") + + matches: list[MatchRow] = [] + for row in rows: + name = catalog_display_name(row) + fpc = catalog_fpc(row) + if not name or not fpc: + continue + score = score_match(query, [name], name) + if score is None: + continue + matches.append( + MatchRow( + fpc=fpc, + prc=catalog_prc(row), + primary_name=name, + secondary_name=None, + source_id=row.get("source_id") or row.get("id") or row.get("sid"), + score=score, + ) + ) + matches.sort(key=lambda row: (row.score, row.fpc)) + return matches + + +def collect_matches(zip_path: Path, query: str) -> list[MatchRow]: + root = load_xml(zip_path) + obj = root.find(".//{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}OBJ") + if obj is None: + raise SystemExit("OBJ block missing in input package") + + matches: list[MatchRow] = [] + for feature in obj.findall(".//{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}CB03"): + fpc = get_text(feature, "FPC") + prc = get_text(feature, "PRC") + if not fpc: + continue + tokens = feature_name_tokens(feature) + primary, secondary = feature_display_name(feature) + score = score_match(query, tokens, primary) + if score is None: + continue + source_id = None + loc = feature.find("{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}LOC") + if loc is not None: + source_id = loc.attrib.get("idref") + matches.append( + MatchRow( + fpc=fpc, + prc=prc, + primary_name=primary, + secondary_name=secondary, + source_id=source_id, + score=score, + ) + ) + + matches.sort(key=lambda row: (row.score, row.fpc)) + return matches + + +def main() -> None: + parser = argparse.ArgumentParser(description="按渔港名字查找 FPC") + parser.add_argument("name", help="渔港名字，支持包含匹配") + parser.add_argument( + "--source", + default=FISH_SRC_DEFAULT, + help="输入渔港原始包，默认 coastline/C09-06.zip", + ) + parser.add_argument( + "--catalog", + default=None, + help="可选的港名对照表（CSV/TSV/JSON/JSONL）；如果原始包里找不到名字，建议用它", + ) + parser.add_argument( + "--exact", + action="store_true", + help="只保留精确匹配", + ) + parser.add_argument( + "--json", + action="store_true", + help="以 JSON 输出", + ) + parser.add_argument( + "--max-results", + type=int, + default=20, + help="最多输出多少条候选结果", + ) + args = parser.parse_args() + + source_path = Path(args.source) + if not source_path.is_absolute(): + source_path = Path(__file__).resolve().parent / source_path + if not source_path.exists(): + raise SystemExit(f"missing input: {source_path}") + + matches: list[MatchRow] = [] + if args.catalog: + catalog_path = Path(args.catalog) + if not catalog_path.is_absolute(): + catalog_path = Path(__file__).resolve().parent / catalog_path + if not catalog_path.exists(): + raise SystemExit(f"missing catalog: {catalog_path}") + matches = load_catalog_matches(catalog_path, args.name) + + if not matches: + matches = collect_matches(source_path, args.name) + if args.exact: + q = normalize_text(args.name) + matches = [row for row in matches if normalize_text(row.primary_name or "") == q] + + matches = matches[: args.max_results] + if not matches: + if args.catalog: + raise SystemExit( + f"no FPC matched for name={args.name!r} in catalog={args.catalog!r}" + ) + raise SystemExit( + f"no FPC matched for name={args.name!r}; " + "this source package may not contain human-readable port names, " + "try --catalog with a port-name mapping file" + ) + + if args.json: + print( + json.dumps( + [asdict(row) | {"score": list(row.score)} for row in matches], + ensure_ascii=False, + indent=2, + ) + ) + return + + for idx, row in enumerate(matches, start=1): + name_bits = [bit for bit in [row.primary_name, row.secondary_name] if bit] + print( + f"{idx}. FPC={row.fpc} PRC={row.prc or '-'} " + f"name={' / '.join(name_bits) if name_bits else '-'} " + f"source_id={row.source_id or '-'}" + ) + + +if __name__ == "__main__": + main() diff --git a/src/pbf/navsea-coastline-fukuoka-saga-200m.html b/src/pbf/navsea-coastline-fukuoka-saga-200m.html new file mode 100644 index 0000000..82181cd --- /dev/null +++ b/src/pbf/navsea-coastline-fukuoka-saga-200m.html @@ -0,0 +1,459 @@ + + + + + + NavSea 全国海岸 / 渔港 / 障碍红格预览 v2.1 + + + + +
+
+
NavSea 全国海岸 / 渔港 / 障碍红格预览
+
+ 版本：v2.1
+ 数据：全国统一库导出
+ 底盘：navsea_japan_coast_grid
+ 默认视图：全国 200x200
+ 底图：GSI 标准地图（mapple） +
+
+ + + + + +
+
+
当前图层采用红色填充和深色边线显示
+
+
等待加载中...
+
+
当前加载：200x200 全国海岸格
+ + + + +