企业官网建设流程全解析

简介：基于Python的GDAL与skimage库，围绕遥感影像的传统图像分割、影像块合并及结果矢量化展开完整实验，适合具有一定图像处理基础、希望将算法迁移到遥感数据的开发者。资料为PDF文档，共1个文件，大小808KB，已吸引1713人学习。文档从read_img读取影像、write_img写回影像入手，依次演示felzenszwalb、slic、quickshift等分割算子，并说明如何借助区域邻接图（RAG）完成分割块合并，最后通过osgeo的ogr模块将分割结果转为矢量要素。作者还专门记录了实验中遇到的棘手问题，附有完整代码与过程分析，便于读者复现并绕开同类坑点。整体篇幅紧凑、代码可读性强，可作为遥感图像分割与矢量化实验的快速参考。

1. 遥感影像分割合并矢量化，这条传统路线依然能打

很多人一提到遥感影像分割，第一反应就是深度学习。但当你手上一没有标注样本、二没有GPU、三需要今天就要出矢量结果时，传统图像分割反而是最快的路。GDAL负责读写和矢量化，skimage 负责把像素变成斑块、再把斑块合并成地物，这套组合对水体、耕地、林地这类光谱相对均质的目标，半天就能走通全流程。这篇笔记就是沿着“读数据 → 分割 → 合并 → 矢量化”的顺序，把参数和坑一次讲透。适合有基础 GIS 概念、想快速交付成果的从业者，也适合刚入门想搞懂分割原理的开发者。

2. GDAL 与 skimage 的桥接：从地理坐标到内存数组

2.1 为什么是这两个库：职责划分与常见弯路

遥感影像和普通图片最大的区别在于它自带地理坐标。skimage 只认 numpy 数组，它不知道影像的投影和仿射参数；GDAL 恰恰弥补了这一点。常见做法是让 GDAL 负责打开文件、读取波段、裁剪、输出矢量，把图像数据转成 numpy 数组后交给 skimage 处理，处理完再把结果写回带地理信息的栅格，最后用 GDAL 的 Polygonize 转成 Shapefile。

很多人第一次做都会踩一个维度坑：GDAL 的ReadAsArray()对多波段影像返回的数组形状是(波段, 行, 列)，也就是 channel-first；而 skimage 几乎所有函数默认输入都是(行, 列, 波段)，也就是 channel-last。如果不做转置，SLIC 会把波段当成空间维度去计算，结果完全不能看。我习惯在读取后立刻统一转成(H, W, C)，后面的代码就干净很多。

另一个常见弯路是试图用 OpenCV 替代 skimage。OpenCV 确实能做分割，而且性能不错，但它没有现成的区域邻接图合并接口，地理信息也得自己手动带回，做遥感分割总感觉差一口气。skimage 自带的slic、watershed、rag_mean_color、merge_hierarchical正好覆盖“过分割 + 合并”的完整链路，少写很多胶水代码。

2.2 数据准备：裁剪、波段取舍、值域归一化

遥感数据基本没有“小图”。一个 20000 x 20000 的影像直接进内存，float32 下三个波段就是 4.8 GB，还没开始分割电脑已经卡死。所以我一般先用gdal.Warp按研究区范围裁剪，或者做一次降采样跑通参数，确定参数后再在原分辨率上处理。

裁剪用到的代码不长，但坐标顺序很容易写反：

from osgeo import gdal src_ds = gdal.Open(r"E:/data/tiff/raw_landsat.tif") # outputBounds 顺序是 (minX, minY, maxX, maxY)，不是 (xmin, xmax, ymin, ymax) warp_opt = gdal.WarpOptions( outputBounds=(120.123, 30.456, 120.456, 30.789), dstSRS="EPSG:4326", resampleAlg="bilinear" ) ds_crop = gdal.Warp(r"E:/data/tiff/crop_test.tif", src_ds, options=warp_opt)

outputBounds是 (minX, minY, maxX, maxY)，很多人按经纬度的直觉写成 (西, 东, 南, 北)，结果裁出来的范围完全不对。dstSRS指定输出坐标系，resampleAlg="bilinear"对多数分割场景够用，如果是分类前的原始数据，用near更保险。

波段取舍也很关键。灰度影像做 RGB 分割效果通常比单波段好，因为颜色信息在分割时是主要线索；如果有近红外波段，可以合成 CIR 或 NDVI 作为额外特征，但波段越多内存压力越大。我一般先只用 RGB 三个波段，跑通后再逐渐加特征。

2.3 可复用的读写模板：保住投影和仿射参数

读入影像时，投影字符串proj和仿射参数geot是矢量化最关键的资产。这两个值一旦丢掉，后面 Polygonize 出来的坐标就是像素行列号，而不是地理坐标。建议一开始就把它俩单独存成变量，和数组一起传下去：

import numpy as np from osgeo import gdal img_ds = gdal.Open(r"E:/data/tiff/crop_test.tif") width = img_ds.RasterXSize height = img_ds.RasterYSize proj = img_ds.GetProjection() geot = img_ds.GetGeoTransform() # 读取前三个波段，shape = (3, H, W) arr = img_ds.ReadAsArray().astype("float32") # 转成 skimage 习惯的 (H, W, C) img = np.transpose(arr, (1, 2, 0))

astype("float32")这一步很必要：原始整型数据在归一化时直接做除法会导致截断误差，而且 skimage 的很多算法内部会转 float，提前转换避免后续隐式转换。但要注意内存翻倍，2.2 里先裁剪的原因就在这里。

归一化直接关系到分割参数是否可迁移。Landsat 反射率数据通常除以 10000，8 位 DOM 除以 255。我习惯统一处理到 [0, 1]：

img_norm = np.clip(img / 10000.0, 0, 1)

归一化不是可有可无的步骤。SLIC 的 compactness 和 RAG 合并的 thresh 都是基于距离的，如果输入值是 0 到 10000 的原始反射率，距离尺度会大好几个量级，参数一换数据源就得重调。统一到 [0, 1] 后，不同来源的数据至少站在同一个尺度上说话。

3. 分割：先做“过分割”，别指望一步到位认出地物

3.1 算法选型：SLIC、分水岭、阈值分割谁更省事

传统分割算法在 skimage 里可选不少，但真正适合遥感整景影像的不多。我列一张常用对比表，方便快速决策：

大津阈值在水体提取里很好用，但它的本质是全局灰度分割，目标一旦有多种地物就失效。分水岭需要把影像转成梯度图，再提供标记点，参数稍微激进一点就会产生大量碎片，整景影像的处理性能也很差。SLIC 是稳定性和效果最平衡的选择：它本身就是为“过度分割”设计的，后面再用区域邻接图合并，正好对应标题里的“分割加合并”流程。

选型有个经验：先想清楚你要的是“直接出地物”还是“出斑块再合并”。如果直接出地物，你会被参数折磨到怀疑人生；如果接受先过分割再合并，整个流程的容错率立刻高一个档次。参数这关多少带点玄学，但方向选对了，后面只是调整量级的问题。

3.2 用 SLIC 跑初始分割：n_segments、compactness、sigma 怎么设

SLIC 的核心思想是把图像划分成颜色一致的超像素，每个超像素内部颜色接近、边界贴合地物边缘。遥感影像上的一小片农田、一小段河流，在 SLIC 里就是一个或多个相邻超像素。

from skimage.segmentation import slic # img_norm 是 (H, W, C) 的 float32 数组，值域 [0, 1] segments = slic( img_norm, n_segments=3000, compactness=30, sigma=1, start_label=1, # 从 1 开始编号，0 留给背景 channel_axis=-1 # 明确通道在最后一维 )

n_segments是最直观的参数，它表示期望生成的超像素数量。数值越大分割越碎，数值越小越粗糙。我一般根据影像面积和目标地物大小估算：1000 x 1000 的影像，目标是 200 米左右的农田地块，n_segments=3000~5000起步。先取偏大的值，后面合并阶段有回旋余地。

compactness控制超像素的形态方正程度。数值越大，超像素越接近正圆形，边缘贴合度越差；数值越小，超像素越能顺着颜色边界延伸，但整体形状会很不规则。遥感影像上我习惯从 30 开始，再按边界贴合度上下试。需要提醒的是，compactness 不是越大越好的单向参数，它和 n_segments 有耦合，调参时要固定一个变量。

sigma是进入 SLIC 前的平滑强度，越大越能抑制椒盐噪声，但也会把细小的田埂和沟渠抹平。对 10 米级别的遥感数据，sigma=1通常够用，亚米级影像可以降到 0.5。记一个检查习惯：把segments叠加到原图上，如果小块边界润到物体内部，就先加 sigma，而不是继续加 n_segments。

3.3 边界贴合度检查：分割结果的三个观察指标

分割参数调没调对，别等矢量化之后后悔。我每次跑完 SLIC 都会第一时间做可视化叠加，这一步能省去后面 80% 的返工：

import matplotlib.pyplot as plt from skimage.segmentation import mark_boundaries mark = mark_boundaries(img_norm, segments, color=(1, 0, 0)) plt.imshow(mark) plt.axis("off") plt.show()

读图时只看三个东西。第一，边界是否贴合目标地物的轮廓，比如农田地块边界与 SLIC 边界偏离是否在半像元以内；第二，目标地物是否被拆成了太多超像素，一个 300 米见方的地块如果被切出几十个超像素，说明 n_segments 偏大；第三，有没有超像素同时跨了两个明显不同的地物，比如一半在水里一半在岸上，说明 compactness 偏大或 sigma 偏小。

这三个观察指标是后续调整的基础。如果处处贴合、每个地物基本有独立超像素，那就直接进合并阶段；如果边界明显穿模，先调 compactness 和 sigma，不要动 n_segments。这个检查环节只需要一分钟，但它决定了整个流程能否稳定复现。

4. 合并：用区域邻接图把碎斑聚成地块

4.1 为什么合并是刚需：SLIC 输出离地物还差一步

SLIC 的输出是过度分割的超像素，它的底层的逻辑是颜色聚类，没考虑任何地物语义。一个完整农田地块在影像上因为喷灌不均、作物长势差异、阴影遮挡，颜色会出现渐进变化，SLIC 会把每个色差明显的子区域都划成独立超像素。这时候如果直接矢量化，交付图会碎得没法看。

合并阶段本质是把“颜色相似的相邻超像素”聚成更大的连通区域。skimage 的实现方式是区域邻接图 RAG：把每个超像素看成图上的一个节点，两个节点相邻就生成一条边，边的权重代表两个区域的合并代价。合并是迭代的，每次选权重最小的边合并，更新相邻关系，直到所有边权重都大于设定的阈值。

这个过程对新手像个黑匣子，但其实核心就一个数字：合并阈值thresh。它决定了两个超像素的颜色差异小到什么程度就值得变成同一个地块。阈值太小，合并效果不明显；阈值太大，河道、田埂这些细长地物会被背景吞掉。调 thresh 没有标准答案，只能靠多次试错，这也是合并且环节最容易被归为“玄学”的部分。

4.2 merge_hierarchical 用法：thresh 与权重函数

skimage 里完成 RAG 合并的标准函数是merge_hierarchical，配合rag_mean_color构建初始图：

from skimage.future.graph import rag_mean_color, merge_hierarchical, merge_mean_color, weight_mean_color # 构建区域邻接图，sigma 表示颜色距离计算时的高斯加权 g = rag_mean_color(img_norm, segments, sigma=5.0) # 迭代合并 merged = merge_hierarchical( img_norm, segments, g, thresh=15, # 合并阈值，值越大合并越激进 rag_copy=True, # 不修改原图 in_place_merge=True, # 允许原地合并，省内存 merge_func=merge_mean_color, weight_func=weight_mean_color )

thresh=15是归一化影像上的经验起点。注意这个值依赖输入值域，如果前面忘了归一化，这里的阈值就会完全失效。实际调参时我按翻倍和减半的节奏找区间：先thresh=15跑一次看合并效果，如果地块还是碎的，直接跳到 30；如果合并过头吞了河道，就退回 10。大概试三轮就能定位合适区间，然后在这个区间里细调到整数。

merge_mean_color是默认的合并函数，它把两个节点合并后的颜色取均值；weight_mean_color计算的是两个相邻区域的颜色距离。如果想加入阈值之外的信息，比如近红外波段的差异，可以自定义 weight_func，返回值越小代表越该合并。下面是个简化的思路：

def weight_func(g, src, dst, n): # g 是图，src/dst 是节点 id，这里以中心颜色距离作为代价 diff = np.abs(g.nodes[src]["mean color"] - g.nodes[dst]["mean color"]).sum() # 额外加上标准差惩罚，纹理差异大的区域不容易合并 texture_penalty = np.abs(g.nodes[src]["std color"] - g.nodes[dst]["std color"]).sum() return diff + 0.5 * texture_penalty

标准差的引入让合并同时考虑颜色和纹理，能明显减少“颜色相近但地物不同”的误合并。不过自定义函数会在每次迭代中被反复调用，超像素数量上万时性能会下降，建议先小图验证，再上大图。

4.3 合并后清理：重编号、去小斑、填洞的顺序不能乱

合并完的 label 数组并不干净。merge_hierarchical 和 remove_small_objects 都会产生标签缺失，甚至把某些像素置为 0，直接影响后续矢量化。这里有一个顺序问题，顺序错了会产生空洞和错位。

先重新编号，再去小斑，最后填洞：

from skimage.segmentation import relabel_sequential from skimage.morphology import remove_small_objects from scipy.ndimage import binary_fill_holes # 第一步：重新编号，确保标签连续 merged_clean, _, _ = relabel_sequential(merged) # 第二步：去掉面积小于 min_size 的碎斑 merged_clean = remove_small_objects( merged_clean.astype("int32"), min_size=200, # 小于 200 像素的斑块直接删除 connectivity=8 # 8 连通判断相邻 ) # 第三步：再次重编号，并把空洞填上 tmp = merged_clean.copy() merged_final = tmp.copy() for lab in np.unique(tmp): if lab == 0: continue mask = tmp == lab merged_final = merged_final.copy() merged_final[binary_fill_holes(mask)] = lab

min_size的单位是像素，取值取决于影像分辨率。10 米分辨率下 200 像素就是 2 公顷，对小地块提取要相应调低，否则地块会被整个删掉。为什么去小斑前必须重编号？因为合并过程中标签会乱跳，不重新编号直接去小斑，remove_small_objects认为的“小”可能统计到多个不相邻区域，误删大片地物。

填洞放在最后的原因也很直观：小斑被删掉后会在原地留下背景洞，如果先填洞再删小斑，洞就被当作真实地物填死了。这个顺序我踩过两次，都是矢量化完成后发现地块中间有镂空，最后回头改流程才解决。

5. 栅格转矢量避坑：5 个最常见的翻车现场

5.1 坐标错乱：投影变换参数在半路丢了

现象：输出 Shapefile 后用 GIS 打开，矢量完全对不上影像底图，坐标显示像是像素行列号。

原因：读取影像时投影字符串和仿射参数没有保存下来，或者用内存栅格保存分割结果时忘了写SetProjection和SetGeoTransform。GDAL 的 Polygonize 是按像素坐标生成几何的，栅格数据里没有地理信息，输出自然就是行列坐标。

解决：从第一步就把proj和geot存成变量，写内存栅格时回填：

from osgeo import gdal, ogr mem_ds = gdal.GetDriverByName("GTiff").Create("", cols, rows, 1, gdal.GDT_Int32) mem_ds.SetProjection(proj) mem_ds.SetGeoTransform(geot) mem_ds.GetRasterBand(1).WriteArray(final_labels)

检查是否写对，可以在 Polygonize 后打印图层范围，和原影像范围做比对。偏差超过一个像元尺寸，说明仿射参数没对上。

5.2 内存爆掉：整景影像直接塞进数组

现象：一个分区的大影像一跑就报 MemoryError，或者 Python 进程直接被系统杀掉。

原因：20000 x 20000 x 3 个 float32 就是 4.8 GB，SLIC 本身还要构建图结构，内存需求通常是数组的好几倍。我见过有人硬跑 2 米分辨率整县影像，最后把服务器拖挂了。

解决：先裁剪到可计算的尺寸，一般 3000 x 3000 以内比较稳妥；参数调试阶段先做 1/2 或 1/4 降采样，参数确定后再用原始分辨率算。降采样跑出来的参数不能直接照搬分辨率，但趋势是对的，先把流程走通，再细化。

5.3 合并过度：河道田埂被背景吞掉

现象：thresh 稍微从 15 调到 18，细长的河道、田埂就从结果里消失了，变成大片单一地块。

原因：纯颜色距离的合并天然对细长地物不友好。河道在影像上颜色和两边耕地可能很接近，但面积小、边界周长长，合并权重只要略低于阈值就会被相邻地块带走。

解决：给权重函数加边界惩罚。预先用 Sobel 或 Canny 算梯度图，在两个相邻区域共同边界上的平均梯度大就增大合并代价，让高边缘区域不容易被吞。另一个更稳的做法是在合并后做形状约束：按面积和周长比识别长条形斑块，如果它周围环境颜色太接近，说明它是线性地物，单独保留。

5.4 地块内空洞：去小斑和填洞的顺序错了

现象：一个完整地块矢量化后中间有一个或多个空白洞，拓扑破碎。

原因：remove_small_objects会把地块内部的小面积色差区直接置为 0，如果这个 0 和周围的背景连成一片，矢量化时就形成孔洞。更麻烦的是，merge 后的标签乱跳，填洞时对错对象操作，越填越乱。

解决：严格按 4.3 的顺序执行，先relabel_sequential再删小斑再填洞。填洞时用 label 掩膜循环逐类处理，不一次对整个数组做binary_fill_holes，否则把不同地物之间的间隙也填了。

5.5 边界锯齿：交付前少做了这两步

现象：矢量边界像长城，在 GIS 里放大后锯齿明显，甲方一眼就看出是机器自动提的，要求重做。

原因：栅格转矢量天然带像素棱角，分割结果不光滑直接转 Shapefile 就必然锯齿。解决手段不是转完再磨皮，而是转矢量前先做一次形态学闭运算，把细小凹槽和凸起磨平：

from skimage.morphology import closing, disk final_labels = closing(final_labels, disk(3))

闭运算之后再 Polygonize，边界的锯齿会少一多半。如果还嫌不够，再用 SimplifyPreserveTopology 做拓扑保持的简化：

for feat in dst_layer: geom = feat.GetGeometryRef() simple = geom.SimplifyPreserveTopology(tol) feat.SetGeometry(simple) dst_layer.SetFeature(feat)

tol就是简化容差，单位是坐标单位（米制投影下就是米），一般取 1~2 倍像元尺寸即可，太大会丢失角点。

6. 矢量化输出与边界优化：最后一步决定交付质量

合并后的final_labels是干净的标签栅格，下一步就是转成真正的矢量成果。用 GDAL 的 Polygonize 把每个标签值转成一个多边形：

from osgeo import gdal, ogr # 1. 把 final_labels 写入带地理信息的内存栅格 rows, cols = final_labels.shape mem_ds = gdal.GetDriverByName("GTiff").Create("mem", cols, rows, 1, gdal.GDT_Int32) mem_ds.SetProjection(proj) mem_ds.SetGeoTransform(geot) mem_ds.GetRasterBand(1).WriteArray(final_labels) # 2. 创建矢量图层 vec_drv = ogr.GetDriverByName("ESRI Shapefile") shp_path = r"E:/data/shp/segments.shp" dst_ds = vec_drv.CreateDataSource(shp_path) dst_layer = dst_ds.CreateLayer("segments", srs=None, geom_type=ogr.wkbPolygon) # 3. 加一个属性字段记录标签编号 field_defn = ogr.FieldDefn("class_id", ogr.OFTInteger) dst_layer.CreateField(field_defn) # 4. 矢量化 gdal.Polygonize(mem_ds.GetRasterBand(1), None, dst_layer, 0)

gdal.Polygonize会把每个独立连通区域转成一个多边形，属性值就是标签编号。这里不需要手动遍历连接关系，GDAL 自己处理，但要注意 SRS 参数如果传了投影字符串，也可以直接用mem_ds.GetSpatialRef()传给 CreateLayer，保持和栅格一致。

输出后还有一个容易忽略的点：面积过滤。RAG 合并之后通常还有少量碎多边形，尤其图斑边缘的小尖角，按面积阈值删一遍再交付会专业很多，而且不要在 WGS84 经纬度下直接算面积，先把矢量投影到米制坐标系再计算。用ogr.Geometry.GetArea()之前先确认坐标单位，否则面积值毫无意义。

我做模拟项目 X 时第一次输出矢量的教训是：线不简化、面不过滤、属性不重命名，结果在交付阶段被要求返工—所有能避免的问题都挤在了同一个环节暴露。从那以后我养成了一个习惯，输出前过一遍检查清单：投影有没有回填、边界有没有叠到影像上看过、小于最小面积的地块有没有清干净、属性字段名是否符合规范。这套流程跑熟之后，半分钟内就能完成一次交付级检查。希望帮到你。

本文还有配套的精品资源，点击获取