基于TIFF深度数据生成独立陆地多边形的技术咨询

核心思路 #

你当前仅提取了2-2.2米水深的边缘点，这种方式无法区分独立岛屿。正确做法是先生成完整的不可航行区域掩码（包含陆地和水深<2米的所有区域），再通过矢量化工具将掩码转换为独立多边形要素，自动识别各个岛屿/陆地区域。

实现步骤与代码 #

1. 生成不可航行区域掩码 #

不可航行区域包含两部分：

• 陆地：对应TIFF中的nodata值
• 水深小于2米的区域：由于TIFF中水深以负数表示（例如-2代表水深2米），因此筛选条件为depth_data >= -2（水深≤2米）

2. 矢量化掩码为独立多边形 #

使用rasterio.features.shapes函数，自动识别掩码中的连通区域，生成每个区域的多边形几何，同时通过TIFF的transform参数直接转换为地理坐标。

完整代码 #

import rasterio
import numpy as np
import geopandas as gpd
from shapely.geometry import shape

# 加载GeoTIFF文件
tif_file_path = "ddm_50m.dybde.tiff"
with rasterio.open(tif_file_path, 'r') as r:
    depth_data = r.read(1)
    transform = r.transform  # 直接读取文件自带的坐标转换参数
    nodata = r.nodata
    crs = r.crs

# 生成不可航行区域掩码：陆地(nodata) + 水深≤2米（depth >= -2）
non_nav_mask = (depth_data >= -2) | (depth_data == nodata)

# 矢量化掩码，提取每个连通区域的多边形
shapes = rasterio.features.shapes(
    non_nav_mask.astype(np.uint8),
    mask=non_nav_mask,
    transform=transform
)

# 将多边形转换为GeoDataFrame，方便后续处理和保存
geoms = []
for geom, _ in shapes:
    geoms.append(shape(geom))

gdf = gpd.GeoDataFrame(geometry=geoms, crs=crs)

# 保存为Shapefile或GeoJSON（按需选择）
gdf.to_file("non_navigable_areas.shp")
# gdf.to_file("non_navigable_areas.geojson", driver="GeoJSON")

关键优化说明 #

• 无需手动定义metadata：直接从rasterio读取的文件对象中获取transform、nodata、crs，避免手动输入错误
• 自动区分独立区域：rasterio.features.shapes会自动识别掩码中的连通组件，每个岛屿/陆地会生成独立的多边形
• 坐标自动转换：通过transform参数，直接将像素坐标转换为TIFF对应的地理坐标（EPSG:3034），无需手动计算

可选调整 #

• 过滤小面积噪点：如果需要忽略微小岛屿或噪点，可添加面积过滤步骤：

  # 过滤面积小于2500平方米的区域（50m分辨率下1个像素为2500㎡）
  gdf = gdf[gdf.geometry.area > 2500]
  

• 提取边缘多边形：若确实需要仅保留2-2.2米水深的边缘，可先对不可航行区域做边缘检测（如用skimage.feature.canny），再矢量化边缘点，但直接矢量化完整掩码的实用性更强。

内容的提问来源于stack exchange，提问作者Lukas Eliassen