海产养殖遥感提取（探索版）

筏式养殖区从影像到矢量面积的完整方法链：纹理能量 + 水域掩膜 → 形态学净化 → 高斯投影下真算面积。附桑沟湾实测样张与 80 点核查方案。

---

一、什么时候用

适用：

• 想知道「遥感能不能做海产养殖」——本技能给出完整方法链与可跑的脚本
• 需要养殖用海范围 / 筏式养殖区 / 网箱的空间分布与面积统计
• 做养殖险的承保支持：哪些海域在养、养了多少、年际有没有变化
• 投标 / 方案演示需要养殖遥感业绩雏形：本技能带桑沟湾实测样张（含分辨率对照与核查样张）

不适用（先明说，别硬套）：

• 单体权属边界：影像只能给"这一片在养"，给不到"这块属于谁"——那要靠确权数据
• 水下生物量 / 存栏量：光学影像看不到水下，估产是另一条科研路线
• 灾损程度判定：台风过后筏架散架在影像上有信号，但"损失几成"需要现场或多次影像时序

这一条要记住：养殖险与种植险的客户群是同一波人（同一家保险公司、同一个农险部门）， 不是新客户。所以海产养殖不是"转行"，是"在同一个客户身上多接一类活"。

二、这块活的钱在哪

市场位置：平台上的农业/保险类付费技能里，没有一个做海产养殖遥感。 这条线现在基本是空的——但"空"的另一面是"没有验证过的付费需求"，所以本技能免费， 先用它占搜索词、导流到有实战背书的种植险场景件。

能拿出去的东西（本技能自带的实测样张）：

| 样张 | 内容 | |---|---| | 成果图 | 桑沟湾筏式养殖区，5 个图斑 / 2,333.1 ha / 34,996 亩（A4 竖版，含经纬网/比例尺/指北针） | | 分辨率选型对照 | 同一片海在 z16(2.4 m) / z17(1.2 m) / z18(0.6 m) 下分别能看出什么 | | 精度核查样张 | 1.0 km × 1.0 km 局部放大 + 提取边界 + 80 个分层随机核查点（判读列留空） |

★ 一条诚实边界：核查点表里的「判读_是否筏式养殖」列是空的。 本技能交付的是抽样核查方案 + 样张，不是精度指标—— 正式精度数字必须由人工逐点判读后统计（界内正确率 / 界外误提率 / 漏提率）。 任何把"提取面积"直接说成"实测养殖面积"的表述都是不实的。

三、步骤 0：口径识别（★ 不许跳过）

开工前三件必锁。

① 目标物口径 —— 你到底要提哪一种？

| 类型 | 影像特征 | 提取难度 | |---|---|---| | 筏式养殖 | 规则条带/栅格纹理，浮筏有周期性 | ★ 最易（纹理就能抓） | | 网箱养殖 | 圆/方小格阵，密集排列 | 中（需形态 + 密度） | | 滩涂/底播 | 边界弱、与滩面反差小 | 难（多光谱或时序） |

不同目标物用完全不同的判据。这一条不锁死，后面参数全都白调。

② 分辨率口径 —— 分不出条带的影像，做不出区块

要提取的最小单元（比如一个养殖区块、或者单体筏架）在影像上至少要有 几十个像素跨度。 实测对照（同一片海）：

| 级别 | 地面分辨率 | 能看出什么 | |---|---|---| | z16 | 2.4 m | 只看得出「这一片有养殖」，分不出条带、更分不出区块 | | z17 | 1.2 m | 能看出筏架条带走向，区块边界仍模糊 | | z18 | 0.6 m | 条带清晰可辨、可计数，能分出不同区块（品种/权属） |

→ 本技能默认 z18（地面 ≈ 0.48 m）。 ★ 分辨率不够是源头缺陷，后处理救不回来——别拿 z16 硬提然后再做形态学。

③ 面积口径 —— 墨卡托面积不守恒

瓦片拼出来的 GeoTIFF 天然是 EPSG:3857（Web Mercator）。 3857 的面积不守恒：高纬度虚高 1/cos²φ。37°N 处虚高约 1.57 倍。

→ 面积必须在真实米制下算：to_crs(高斯 3°带) 后再算，除以 666.6667 得亩。 → 带号按经度选：122°E → EPSG:4529（CM 123°E）；132°E → 4553。 → 不要用图层自带的 CRS 量面积，也不要用计算几何工具默认口径。

识别不出时的三级兜底

1. 从用户原话提（要提什么、哪个海域、什么用途）
2. 从影像推（纹理特征像不像筏式条带、水域掩膜占比多少）
3. 仍不确定 → 同时出三种参数的对比图，让用户挑，不要静默选一组

四、输入契约

必备

| # | 输入 | 要求 | |---|---|---| | 1 | 影像 | 网络可达的在线瓦片（Esri World Imagery / 高德卫星）或自备的高分影像 GeoTIFF | | 2 | AOI 范围 | WGS84 bbox（W S E N），或一份海域/岸线矢量 |

可选（有则更好）

| # | 输入 | 用途 | |---|---|---| | 3 | 行政区 / 海域权属界 | 统计按权属归并 | | 4 | 历史影像 | 做年际变化监测 | | 5 | 现场照片 / 台账 | 精度核查的独立证据 |

环境依赖

rasterio, geopandas, shapely, pyproj, opencv-python, numpy, Pillow, matplotlib

字段口径

成果图斑只需三个面积字段，且必须标明算法：

area_m2   = 几何面积（高斯投影下）
area_ha   = area_m2 / 10000
area_mu   = area_m2 / 666.6667

五、核心计算

2.1 抓图（grab_tiles.py）—— 拼成一张 GeoTIFF

python grab_tiles.py --src esri --zoom 18 --bbox 122.54 37.08 122.62 37.12 --out area_z18

要点：

• 瓦片编号换算：px = 2 * πR / (2^z) / 256（每像素的墨卡托米）
• ★ 方位不要交叉：西 → 最小列、北 → 最小行。 反过来算会得出负的瓦片总数（实测算出 -2024 片）
• 缺片兜底：抓失败的格子留黑，不要静默跳过（会在地面上留空洞却看不出来）
• 并发 24 线程足够；抓完打印「缺失瓦片 N 片」并给出前几个坐标

2.2 提取（extract_raft.py）—— 纹理能量 + 水域掩膜

python extract_raft.py area_z18.tif --out area --k 9 --tex-q 68 --min-ha 0.5

四步：

1. 水域掩膜：(蓝 − 红) > 5。筏式养殖一定在水面上，这一条先砍掉岸上所有干扰 （房屋、道路、农田纹理都很强，不先砍会全进来）
2. 纹理能量：局部标准差 std = sqrt(boxFilter(g²) − boxFilter(g)²)，窗口 k=9 —— 筏架的周期性条带让局部方差显著抬高
3. 阈值 + 形态学：
   • 阈值在水域子集内取分位数（默认 q68），不在全图取
   • 闭运算 25 px 把条带连成片 → 开运算 9 px 去掉孤立噪点
   • 连通域面积筛选：小于 min_ha（默认 0.5 ha）剔除
4. 矢量化 + 面积：rasterio.features.shapes → geopandas → to_crs(4529) 后算面积（不是 3857）

参数怎么定：在候选网格里选（k ∈ {9,13,17} × tex_q ∈ {60,68,76}）， 看哪个组合的图斑边界最贴合影像上的实际条带范围。别拍脑袋定一个。

2.3 出图与核查

python make_map.py      area_z18.tif area_rafts.gpkg 成果图.png
python make_res_compare.py                                  # 分辨率选型对照
python make_check.py                                        # 精度核查样张 + 点表

• 成果图按 A4 竖版：影像底图 + 图斑半透明填充 + 边界线 + 经纬网 + 图内比例尺 + 指北针 + 无黑框图例
• 核查样张：1 km × 1 km 局部放大，分层随机 80 点（保证界内界外都有样本）， 红圈 = 界内、青圈 = 界外，逐个编号
• 点表导出后判读列留空，交人工逐点在原分辨率影像上判读

六、机器自证（不合格不许交付）

① 面积口径自证（最关键）

g = gpd.read_file("area_rafts.gpkg")
assert g.crs.to_epsg() == 4529, f"CRS 不是 4529，是 {g.crs}"
assert abs(g.geometry.area.sum() / g["area_m2"].sum() - 1) < 1e-6, "面积字段与几何不一致"

别信元数据标着 4529 就完事——set_geometry(crs=) / gdf.crs= 只改标签不改坐标。 唯一真转坐标的是 to_crs()。

② 分辨率自证

图斑的等效半径 ÷ 地面 GSD 应 ≥ 10 像素。 若最小图斑只有几个像素宽 → 该图斑是噪声，不是养殖区（要么调大 min_ha，要么换更高分辨率）。

③ 影像几何自证

拼图的 GeoTIFF 在 GIS 里打开，与已知岸线/岛礁套一下。 实测踩过：transform 少写 256 倍 → 图幅地面尺寸错成 5,772 km（一个省那么大）。

④ 抓图完整性自证

• 打印缺失瓦片数；缺失 > 0 时在成果图上标注"本区域有数据缺口"
• 别把缺片当成"这里没有养殖"

⑤ 顺序自证

把参数故意调到明显错（如 tex_q=10），看输出是不是变成一大片—— 如果调了没反应，说明阈值根本没生效在正确子集上。

七、输出物

<项目>\
  area_z18.tif               # 拼接影像（EPSG:3857）
  area_rafts.gpkg            # 图斑矢量（面积字段在高斯投影下算好）
  area_rafts.csv             # id / area_m2 / area_ha / area_mu
  area_mask.npy              # 二值掩膜（复算用）
  成果图.png                 # A4 竖版专题图
  分辨率选型对照.png          # z16 / z17 / z18 三栏对照
  精度核查样张.png            # 局部放大 + 边界 + 80 点
  精度核查点表.csv            # id / in_extract / lon / lat / 判读列(留空)

面积字段口径：area_m2 为高斯投影下几何面积，area_ha = /10000， area_mu = /666.6667。引用时一律取 area_mu 字段，勿用计算几何。

八、熔断与错误码

| # | 情形 | 动作 | |---|---|---| | 1 | 最小图斑只有几像素宽 | 停：分辨率不够，换 z18 或更高分影像，别用形态学硬凑 | | 2 | 抓图缺失 > 5% 且集中在边缘 | 补抓或在图幅上标注缺口，不许当成"无养殖" | | 3 | 提取结果连成一大片、分不出区块 | 阈值取得太低或形态学闭运算过大，回到 2.2 重调参数 | | 4 | 岸上区域也被提取出来 | 水域掩膜没生效——检查波段顺序（B/R 有没有弄反） | | 5 | 有人要求把提取面积直接当作养殖面积上报 | 拒绝：须先人工判读、统计精度，否则是不实交付 |

精度核查的诚实说明

• 本技能的核查点表判读列是空的——这是故意的，不是遗漏。
• 正式项目必须：人工逐点判读 → 算 界内正确率 / 界外误提率 / 漏提率 → 写进报告。
• 抽样必须分层（界内 + 界外都抽），只抽界内点会系统性高估正确率。

九、常见坑

| # | 症状 | 根因 | 动作 | |---|---|---|---| | 1 | 图幅地面尺寸算成 5,772 km | GeoTIFF transform 少乘了 256 倍 | px = tile_m / 256，不是 tile_m | | 2 | 瓦片总数算成负数（如 -2024） | 经纬度与瓦片行列交叉配对 | 西 → 最小列、北 → 最小行 | | 3 | 面积比预期大 1.5~1.6 倍 | 在 EPSG:3857 上量面积（面积不守恒） | to_crs(4529) 后再量 | | 4 | 改了 CRS 但数值没变 | 用了 gdf.crs = ...（只改标签） | 必须用 to_crs() | | 5 | 岸上房屋/农田全被提取 | 没做水域掩膜，或 B/R 波段弄反 | 先 (B−R) > 5 砍掉陆地 | | 6 | 阈值调了没反应 | 分位数取在了全图而非水域子集 | np.percentile(tex[water], q) | | 7 | 条带被连成一坨、区块分不开 | 闭运算核太大 | 25 px 起步，按 GSD 折算而不是固定像素 | | 8 | 图幅上有看不见的空洞 | 缺片被静默跳过 | 打印缺失清单；缺片留黑不填充 | | 9 | 经纬网标签压在图上挡内容 | 网格标注没加半透明底 | 给标签加 bbox=dict(fc='#000000', alpha=0.55) | | 10 | matplotlib 出图只剩一小块有影像 | imshow(extent=子区) 后被矢量图层撑开 | 画完图层再 set_xlim/set_ylim | | 11 | 出图字号/图例忽大忽小 | 版面参数按像素写死 | 一律按图幅比例定尺寸（Wd*0.013 之类） |

十、免费试算

本技能免费，所以这一节讲的是「你能白拿到什么、以及哪些事别指望它」：

白拿的部分：九节完整方法链 + 5 个可直接改路径运行的脚本 + 三张实测样张 （成果图 / 分辨率对照 / 核查样张）+ 11 条带症状的坑表。跑通一次的成本是磁盘和带宽，不是钱。

别指望它做的：

• 给出合规的精度指标——那需要人工判读，技能里只给了抽样方案和点表
• 判断权属归属——影像给不了"这块属于谁"
• 估水下产量——光学影像看不到水下

如果你要的是能写进报告、能过甲方验收的东西： 本技能是探索版，不作交付承诺。真正有实战背书、带熔断清单和自证断言的是 承保验标 / 长势监测 / 受灾定损那几个场景件——那些是从几十个已交付项目里长出来的。

十一、隐私与数据

• 本技能不收集、不上传、不存储任何数据。抓图脚本走的是公开在线瓦片服务， 不涉及任何客户或承保信息。
• 演示案例为公开海域（桑沟湾），不含任何客户信息、保单信息或个人身份。
• 抓取在线瓦片请遵守对应服务的使用条款；商用交付应改用自购影像或正规数据源， 不要长期依赖公开瓦片做商业成果。
• 成果图上不要叠加养殖户姓名/权属编号等敏感属性——本技能的成果图只画边界与面积。