简介一份面向遥感影像处理与地理信息开发者的技术文档围绕Python中GDAL与skimage库的结合应用系统展示传统图像分割、图像块合并和矢量化实现思路。文档基于作者此前的分割算法实验把felzenszwalb、slic、quickshift等算法迁移到遥感影像场景并给出read_img、write_img、矢量驱动识别等可直接参考的代码函数帮助读者完成栅格数据读取、数据写回和矢量层构建。资源为单一PDF文档大小808KB内容以代码实现为主虽然只有1个文件但包含了从影像读取到分割、合并再到矢量化的完整方法记录并附有作者在实验过程中遇到的棘手问题描述与求助信息。目前该资源已有1713人学习下载适合具备一定Python基础、希望快速上手GDAL与skimage进行遥感影像分割及矢量化的开发者和研究人员参考包体不大但核心步骤集中可作为项目起步时的速查资料。1. 遥感影像分割到矢量化这套 Python GDAL skimage 流程能直接把栅格变成面要素很多做遥感的人卡在同一个地方分割结果在 Python 里怎么看都漂亮但业务上要的是能落库的矢量边界。这套代码把遥感影像的传统图像分割、区域合并和矢量化串成了一条完整链路——读入一张带坐标的 tif用 quickshift 做超像素粗分割再通过 RAG 区域邻接图合并碎块最后用 gdal.Polygonize 输出 shapefile 面要素。整个过程不依赖任何商业 GIS 软件只要 Python 环境里装好 GDAL 和 skimage 就能跑。适合做样本标注、地物粗提取、分割算法对比验证的从业者。我完整跑通后最大的感触是算法本身并不难真正耗时间的是波段顺序、数据类型映射和栅格转矢量这三类“非算法”环节。2. GDAL 影像读写read_img 与 write_img 里的坐标、dtype 和波段顺序这一章解决全流程的地基怎么把一张 tif 读成 numpy 数组处理完再带上地理坐标写回磁盘。代码里的 read_img 和 write_img 是两个可以直接复用的函数但里面有几个参数值得拆开讲清楚否则后续所有步骤都会在数据形状上翻车。2.1 read_img为什么返回宽高、投影和仿射变换六个值def read_img(filename): dataset gdal.Open(filename) im_width dataset.RasterXSize im_height dataset.RasterYSize im_geotrans dataset.GetGeoTransform() im_proj dataset.GetProjection() im_data dataset.ReadAsArray(0, 0, im_width, im_height) del dataset return im_width, im_height, im_proj, im_geotrans, im_datagdal.Open打开的是栅格数据集RasterXSize和RasterYSize分别是列数和行数注意不要和 numpy 的行列习惯搞混。GetGeoTransform()返回一个长度为 6 的元组六个值依次是左上角 X 坐标、像元宽度、X 方向旋转、左上角 Y 坐标、Y 方向旋转、像元高度。对绝大多数标准正射影像旋转项都是 0所以你要记住的实际上是四个值左上角坐标、像元宽、像元高。这六个值后续在write_img里必须原样写回否则输出的 tif 在 GIS 里会丢失地理参考或者位置偏移。ReadAsArray(0, 0, im_width, im_height)读出来的数组形状是(bands, rows, cols)波段在最前面。这个顺序对后续用 skimage 是个关键差异因为 skimage 的分割函数普遍期望输入是(rows, cols, channels)的形状。代码里后面用temp im_data.transpose((2, 1, 0))就是在做这个轴交换。最后del dataset不是可有可无的操作。GDAL 在 Windows 下对文件句柄的释放比较敏感尤其是后面还要写回同一个目录下的 shapefile如果不释放轻则报权限错误重则整个脚本卡住不动。我一般会在read_img里强制写删除因为这个代价几乎为零。2.2 write_img 的 dtype 映射int8、int16 和浮点的坑def write_img(filename, im_proj, im_geotrans, im_data): if int8 in im_data.dtype.name: datatype gdal.GDT_Byte elif int16 in im_data.dtype.name: datatype gdal.GDT_UInt16 else: datatype gdal.GDT_Float32 ...这段 dtype 判断代码有一个隐藏细节判断条件是int8 in im_data.dtype.name而 numpy 的uint8类型名字是uint8字符串子串int8恰好匹配所以能进入GDT_Byte分支。如果你的数组是int16但值是负数GDAL 的GDT_UInt16是存不进去的写入时会被截断或报错。常见做法是在写之前统一做类型转换比如im_data.astype(np.uint8)或按百分比拉伸到 0-255 再写。另一个常见的坑是分割标签图labels2的数据类型是int64直接传给write_img时会落入GDT_Float32分支。这种 tif 能打开但标签值在 GIS 里会被当作浮点显示后续做分类统计时容易出现精度问题。代码里在保存qs_label.tif前没有显式转类型我建议你在自己的流程中补一步labels2.astype(np.int32)避免不必要的隐患。如果你发现写出的 tif 在 QGIS 里打开是全黑或全白不要急着怀疑分割算法先检查数组是不是float32再看最小值是不是 0、最大值是不是远大于 255。很多影像读出来是 16 位反射率数据值域在 0-10000直接写成 Byte 就会溢出成噪声。处理手段是在写之前做线性拉伸或直方图拉伸。2.3 二维/三维分支与多波段写回的细节if len(im_data.shape) 3: im_bands, im_height, im_width im_data.shape else: im_bands, (im_height, im_width) 1, im_data.shape driver gdal.GetDriverByName(GTiff) dataset driver.Create(filename, im_width, im_height, im_bands, datatype) dataset.SetGeoTransform(im_geotrans) dataset.SetProjection(im_proj) if im_bands 1: dataset.GetRasterBand(1).WriteArray(im_data) else: for i in range(im_bands): dataset.GetRasterBand(i 1).WriteArray(im_data[i])这个分支逻辑本身不复杂但有一个容易误用的点else分支里的im_bands, (im_height, im_width) 1, im_data.shape实际上就是把im_data.shape解包成(im_height, im_width)。如果传入的数组是三维且形状是(rows, cols, bands)这个分支会直接误判为单波段导致后面driver.Create创建出的数据集波段数错误写入时也会报维度不匹配。正确的做法是写回前对图像数组统一调用np.transpose((2, 0, 1))把波段维度换回第一位再交给write_img。代码里作者在保存mark0时用的是mark0.transpose((2, 1, 0))恰好和这里的预期一致。如果是自己拼接的 RGB 数组务必验证shape是不是(bands, rows, cols)后再调用。提示多波段写回时GDAL 波段编号从 1 开始不要写成dataset.GetRasterBand(i)那会越界。3. quickshift 粗分割超像素参数与格网骨架的一次成型分割这层是整条流程的核心。代码导入了felzenszwalb、slic、quickshift三个算法但主流程用的是 quickshift这里结合参数说明为什么选它以及骨架提取为什么要绕一道弯。3.1 quickshift 与 slic、felzenszwalb三者如何取舍from skimage.segmentation import felzenszwalb, slic, quickshift from skimage.segmentation import mark_boundaries三个算法里slic生成的是规则且紧致的超像素区域数量能通过n_segments精确控制适合视觉上比较规整的样本标注但对地物边界复杂的山体、水体它容易把细长的河流或田埂切断。felzenszwalb在超大规模影像上速度优势明显但它的scale参数和最终区域尺度的对应关系不直观调参基本靠试。quickshift 属于基于密度峰值的聚类参数ratio可以直观控制颜色和空间距离的权重对遥感影像这种颜色分布相对连续的场景更友好。我一般会先用 quickshift 跑一个初始结果如果超像素过于细碎再用graph.rag_mean_color做区域合并而不是换算法重跑。原因是 quickshift 之后的合并是基于颜色统计的和 RAG 的邻接关系天然匹配而 slic 的超像素形状太规整合并后边界容易显得锯齿感强。3.2 kernel_size、max_dist、ratio 分别控制什么temp im_data.transpose((2, 1, 0)) segments_quick quickshift(temp, kernel_size3, max_dist6, ratio0.5) mark0 mark_boundaries(temp, segments_quick)kernel_size是高斯核的大小控制局部密度估计的窗口。它越大分割结果越平滑细节越少对高分辨率遥感影像我建议从 5 起步如果目标地物很细碎再用 3。max_dist表示两个像素在特征空间中允许连接的最大距离直接决定了最终超像素的尺度。它越大越多的像素会被归到同一个区域代码里的 6 属于比较保守的取值得到的碎块数量多为后续合并留了空间。ratio是颜色距离和空间距离的权重比。0.5 表示两者各占一半如果目标地物颜色差异明显但空间上分散可以把 ratio 调小到 0.2-0.3让算法更偏重颜色聚合如果影像中同类地物颜色变化大、但空间上连续ratio 可以放大到 0.8 左右。mark_boundaries返回的是在原图上叠加了区域边界的高亮图它是三维数组最后一位是 RGB 三通道后面保存边界格网时就用到了这一点。3.3 格网骨架提取为什么要先转单像素grid_path temp_path qs_grid0.tif grid0 np.uint8(re0[0, ...]) write_img(grid_path, im_proj, im_geotrans, grid0) skeleton morphology.skeletonize(grid0) border0 np.multiply(grid0, skeleton) ret, border0 cv2.threshold(border0, 0, 1, cv2.THRESH_BINARY cv2.THRESH_OTSU)mark0是三通道图作者直接取re0[0, ...]也就是红色通道得到的是所有分割边界组成的格网。但这一步得到的边界线通常有 2-3 个像素宽如果直接拿去做矢量化输出的线永远不会贴合真实边界而是产生双线。morphology.skeletonize的作用是把宽度不均匀的线条细化为单像素骨架返回的布尔数组表示每个位置是否属于骨架线。np.multiply(grid0, skeleton)是把原始格网值和骨架位置相乘相当于只保留骨架位置上的像素值然后通过 OTSU 自动阈值把非零值归一化到 0 和 1。这个二值图才是后面矢量化能用的输入。很多人在这一步图省事直接用grid0 0当二值图结果多边形边缘会有大量锯齿和小孔隙后处理成本更高。注意如果影像尺寸很大超过 8000×8000建议先切块跑 quickshift再把分割标签偏移后拼回否则内存很容易吃满。常见做法是每块保留 32 像素重叠做完分割后只取中心区域能有效避免块接缝处的错位。4. RAG 区域合并用 thresh 回调控制超像素合并粒度粗分割得到的超像素通常碎得没法直接用需要按颜色相似度合并相邻区域。这里的做法是构建区域邻接图RAG然后做分层合并。这章重点讲节点权重、合并回调和阈值的选择。4.1 rag_mean_color 节点与边图里到底存的什么g graph.rag_mean_color(temp, segments_quick) labels2 graph.merge_hierarchical( segments_quick, g, thresh5, rag_copyFalse, in_place_mergeTrue, merge_funcmerge_mean_color, weight_func_weight_mean_color )rag_mean_color会为每一个超像素建立一个节点节点属性包括平均颜色mean color、像素总数pixel count和颜色总和total color。边连接的是空间上相邻的节点初始边权是两个节点平均颜色差的欧氏距离。这里的“相邻”不是图像坐标上的上下左右而是基于超像素边界的共边关系所以即使两个区域在颜色上完全一样只要空间上不相邻也不会被合并。merge_hierarchical的核心机制是在图上反复找边权最小的边如果它小于thresh就把两个节点合并然后更新所有邻接边的权重。thresh5意味着允许合并颜色距离在 5 以内的相邻区域。对 0-255 的 RGB 图像来说这个值属于中等偏保守能保留较多细节。4.2 merge_mean_color 回调与 in_place_merge 的坑def merge_mean_color(graph, src, dst): graph.nodes[dst][total color] graph.nodes[src][total color] graph.nodes[dst][pixel count] graph.nodes[src][pixel count] graph.nodes[dst][mean color] (graph.nodes[dst][total color] / graph.nodes[dst][pixel count])这个回调在每次合并前被调用负责把src节点的颜色统计合并进dst节点。关键点在于它维护的是total color和pixel count然后重新计算mean color而不是直接对两个平均颜色做算术平均。如果直接平均面积大的区域会被小区域的颜色拖累导致合并偏向小区域最终结果会出现“细小地物吞噬大片地物”的怪象。def _weight_mean_color(graph, src, dst, n): diff graph.nodes[dst][mean color] - graph.nodes[n][mean color] diff np.linalg.norm(diff) return {weight: diff}weight_func在节点合并后重新计算dst和邻居n之间的边权。这里有个容易忽视的细节in_place_mergeTrue会直接修改传入的图对象g如果后面还需要原始 RAG 做其他分析记得改成rag_copyTrue否则原图会被破坏。我最初跑的时候没注意这个参数后续想复用g做可视化时发现节点数已经变了还以为是算法 bug。4.3 合并阈值怎么调从碎块到面要素的粒度控制合并后输出前需要先检查合并效果是否合理。一个实用的调试手段是输出合并前后的区域数量region_num_before len(np.unique(segments_quick)) labels2 graph.merge_hierarchical( segments_quick, g, thresh5, rag_copyFalse, in_place_mergeTrue, merge_funcmerge_mean_color, weight_func_weight_mean_color ) region_num_after len(np.unique(labels2)) print(before:, region_num_before, after:, region_num_after)thresh的选择和影像地物尺度强相关。做精细地物边界提取比如房屋轮廓thresh1-3足够做林地、农田等大面状地物粗提取thresh10-20能把大块超像素合并成完整地块。代码里的thresh5属于通用初始值适合城市和郊区混合的影像。合并后的标签图直接写盘和可视化label_rgb2 color.label2rgb(labels2, temp, kindavg) rgb_path temp_path qs_label.tif lb labels2.transpose((1, 0)) write_img(rgb_path, im_proj, im_geotrans, lb) mark mark_boundaries(label_rgb2, labels2) save_path temp_path qs_seg.tif re mark.transpose((2, 1, 0)) write_img(save_path, im_proj, im_geotrans, re)label2rgb的kindavg表示每个分割区域用该区域的平均颜色填充这样保存出来的图可以直接当分类示意图用。注意lb labels2.transpose((1, 0))把行和列调换了一下再写盘这一步是因为labels2的形状是(rows, cols)但 GDAL 期望输入的二维数组按(rows, cols)理解时逻辑上和原图是一致的这里转置会让输出和边界标记图在 GIS 中显示的轴向对齐。实际使用之前建议先用小图验证轴向是否正确。5. 避坑排查栅格转矢量的五个典型问题从栅格到矢量是整个流程里最容易“翻车”的一段作者在代码里也专门留了一段注释说明自己遇到的棘手问题。这里我把跑通流程时踩过的五个坑按现象、原因、解决的思路写清楚。5.1 栅格线直接转矢量线为什么只能得到双线现象把分割边界栅格直接转矢量线得到的每条边界都是两条平行线且都不在真实边界上。原因mark_boundaries画出的边界本身是 2-3 个像素宽的“带”栅格矢量化时每个像素面按八邻域规则连接带的左右边缘各自生成一条线中心线反而丢失。GDAL 里没有现成的栅格线转矢量线函数直接用边界带做 Polygonize 只会得到沿着带宽两侧的曲面。解决先对边界栅格做skeletonize细化到单像素再转矢量实在不行就折中先 Polygonize 成面后面再用面转线处理。代码采用的就是后者这也是为什么最终输出的是面要素而不是线要素。5.2 Polygonize 的 mask 波段面要素里出现大量细碎条带现象Polygonize 之后的面要素数量非常多很多是沿骨架线生成的细长条带或小碎面。原因skeletonize 后的二值图存在少量孤立像素和分支Polygonize 会把每个连通区都转成一个面要素没有做面积过滤。解决在 Polygonize 之前先对骨架二值图做一次形态学闭运算或按连通区面积过滤掉小于阈值的区域from skimage import morphology cleaned morphology.remove_small_objects(skeleton.astype(bool), min_size50)min_size按像素数控制具体取多少取决于影像分辨率一般取 30-100 之间。这一步能显著减少矢量面的碎块数量但要注意不要设太大否则细长地物边界会被整个删除。5.3 write_img 遇到 float 数组输出 tif 在 GIS 里全黑现象分割标签图写入后打开是纯黑或花屏。原因标签数组是int64或float64落入写函数的GDT_Float32分支而标签值域可能只有 0-几十GIS 默认拉伸范围不匹配显示出来就是黑的。解决写遥感分割结果前对标签图统一做astype(np.uint8)或astype(np.int32)并且先打印数组的dtype和min/max确认值域合理再写盘。5.4 波段顺序错位分割结果出现横向条纹现象quickshift 的分割结果沿着影像行方向出现明显的条带状看起来像压缩损坏。原因ReadAsArray读出来的数组是(bands, rows, cols)如果没有transpose((2, 1, 0))就直接喂给 skimage算法会把波段当作行来处理导致分割目标完全错误。解决入口处统一做一次形状断言assert temp.shape[2] im_data.shape[0], 波段顺序需要先转置常见的错误是把转置写成transpose((1, 2, 0))得到的是(cols, rows, bands)虽然能运行但分割结果的走向会偏 90 度这类错误特别隐蔽。5.5 shapefile 输出被占用第二次运行直接报错现象第一次运行成功第二次运行时CreateDataSource失败提示无法创建、文件被占用。原因上次运行的.shp和.dbf等文件被 GIS 软件或代码里的数据源对象占用GDAL 不会自己删除覆盖。解决在创建数据源前先尝试删除旧的输出文件import os for ext in [.shp, .shx, .dbf, .prj]: f dst_filename.replace(.shp, ext) if os.path.exists(f): os.remove(f)如果还删不掉检查代码是否有dst_ds、src_ds没有关闭调试时确保在脚本末尾调用del dst_ds和del src_ds。6. 面转线与后处理让分割结果变成真正能用的边界线上一章提到栅格线直接转矢量线会在 GDAL 里得到双线作者最终用 Polygonize 先得到面再留了面转线、线简化、线平滑的待办。这一章给这几个步骤一个可执行的落地方案也是这类分割矢量结果落到业务前的最后一步。from osgeo import ogr src_ds ogr.Open(polygon_shp, 0) src_lyr src_ds.GetLayer(0) line_ds ogr.GetDriverByName(ESRI Shapefile).CreateDataSource(line_shp) line_lyr line_ds.CreateLayer(boundary, geom_typeogr.wkbLineString) for feat in src_lyr: geom feat.geometry() if geom is None: continue boundary_geom geom.GetBoundary() if boundary_geom is not None: out_feat ogr.Feature(line_lyr.GetLayerDefn()) out_feat.SetGeometry(boundary_geom) line_lyr.CreateFeature(out_feat)这段代码的思路是遍历每个面要素调用GetBoundary()把多边形的外环提取出来作为线要素。对带洞的面还需遍历GetInteriorRingCount()把内环也用GetInteriorRing(i)取出来否则会把洞的边界漏掉。转成线之后最直接的优化是线简化。GDAL 的Geometry.Simplify(tolerance)用的就是 Douglas-Peucker 算法tolerance 建议取当前影像像元尺寸的 1-2 倍。比如影像分辨率是 2 米tolerance 取 2-4 米能去掉大部分锯齿又不会把地物拐角拉平。先简化再平滑顺序不要反因为先简化可以大幅减少平滑阶段的顶点数运行更快拟合也更稳。我自己的习惯是保存最终结果前把简化前后的线要素数量和各线长度打印出来对比一次如果总长度波动超过 5%说明 tolerance 设得过大需要回退。从那以后我每次做分割到矢量化实验都强制先拿一小块 1000×1000 的影像跑通全链路再上全图这个习惯帮我躲过了不少内存和坐标轴方向的坑。希望帮到你。本文还有配套的精品资源点击获取