1. 从一个实际需求说起：为什么需要 Polygonize
做遥感影像处理或者栅格数据分析的朋友，大概率都遇到过这样一个场景：手里有一张分类结果图，或者一张二值化的掩膜图，领导或者下游同事要的不是栅格，而是矢量面文件。比如把耕地分布图转成shp，把水体提取结果转成GeoJSON，把变化检测的斑块转成可编辑的多边形。这时候，gdal.Polygonize()就是绕不开的一个函数。
我自己第一次接触这个函数，是在做一个土地覆盖分类后处理的活儿。当时分类结果是一张GeoTIFF，每个像素值代表一种地类，我需要把相同地类的连通区域提取成独立的多边形，并且保留地类属性。一开始想着用ArcGIS手动转，但数据量一大就卡得不行，而且没法批量。后来翻GDAL的Python绑定文档，找到了osgeo.gdal.Polygonize()，写了几十行代码就搞定了，效率比手动操作高了不知道多少倍。
这篇文章就围绕osgeo.gdal.Polygonize()这个函数，把它的参数含义、调用逻辑、完整源码实例、常见坑点全部拆开讲清楚。不管你是刚接触GDAL的新手，还是已经用过但总觉得参数没搞明白的老手，应该都能从里面找到有用的东西。文章会给出可直接运行的完整代码，也会解释每一步为什么这么写，还会分享一些文档里不会写的实操经验。
2. Polygonize 函数的核心机制与参数拆解
2.1 函数签名与参数逐个说明
先看函数签名。在Python绑定里，gdal.Polygonize的调用形式是这样的：
gdal.Polygonize(srcBand, maskBand, outLayer, dstFieldIdx, options=[], callback=None)参数不多，但每一个都有讲究，我逐个拆解。
srcBand：源栅格波段。注意，这里传的是一个Band对象，不是Dataset。很多人第一次用的时候直接把dataset传进去，结果报类型错误。正确做法是先dataset.GetRasterBand(1)拿到波段。这个波段里的像素值就是后续多边形的属性来源。
maskBand：掩膜波段。这个参数可以传None，也可以传一个Band对象。它的作用是控制哪些像素参与多边形化。如果传None，那么源波段里所有非零像素都会被处理，零值像素被当作背景忽略。如果传一个掩膜波段，那么掩膜中非零的位置才会被处理。这个设计很实用，比如你只想处理某个行政区范围内的像素，就可以用行政区栅格做掩膜。
outLayer：输出矢量图层。这是一个Layer对象，需要提前创建好，并且定义好字段结构。Polygonize 不会自动帮你创建图层，也不会自动添加字段，这些都要你自己准备好。字段的类型和数量要和dstFieldIdx对应上。
dstFieldIdx：目标字段索引。这是一个整数，表示多边形化之后，源波段的像素值要写入输出图层的第几个字段。注意，这个索引是从0开始的。如果你传-1，表示不写入任何字段，只生成几何。实际使用中一般都会传一个有效索引，把像素值保留下来作为属性。
options：可选参数列表。这个参数用得比较少，一般传空列表就行。它对应的是GDAL命令行工具里的选项，比如"8CONNECTED=8"可以指定8连通，默认是4连通。这个后面会详细讲。
callback：回调函数。用于进度提示，传None就行，大数据量的时候可以传一个函数来显示进度。
2.2 4连通与8连通的本质区别
连通性是Polygonize里一个非常关键的概念，但很多人容易忽略。默认情况下，Polygonize 使用4连通，也就是说一个像素只和上下左右四个方向的像素算作连通。如果两个同值像素只是斜对角相邻，它们会被当作两个独立的多边形。
而8连通则把斜对角也算作连通。这个区别在实际数据上影响很大。举个例子，一条斜着走的道路，在栅格上表现为阶梯状的像素排列。用4连通，这条道路会被切成很多个小多边形；用8连通，它会被合并成一个完整的多边形。
设置8连通的方式是在options里传"8CONNECTED=8"：
gdal.Polygonize(srcBand, maskBand, outLayer, 0, options=["8CONNECTED=8"])我个人的经验是，做地块提取、水体提取这类需要保持区域完整性的任务，用8连通更合适；做像素级精确统计的任务，用4连通更严谨。这个选择没有绝对的对错，取决于你的业务需求。
2.3 输出图层的字段设计逻辑
输出图层的字段设计是很多人容易踩坑的地方。Polygonize 本身不会创建字段，它只是把像素值写入你指定的字段索引。所以你必须提前用CreateField创建好字段。
字段类型的选择也有讲究。源波段的像素值如果是整数类型（比如Byte、Int16、Int32），字段用OFTInteger就行。如果是浮点类型，字段要用OFTReal。如果类型不匹配，写入的时候可能会报错或者丢失精度。
还有一个细节：如果你需要保留多个信息，比如像素值加上多边形面积，那面积字段需要在Polygonize完成之后再计算，因为Polygonize执行的时候几何还没生成完。这个顺序不能搞反。
3. 完整可运行的源码实例
3.1 环境准备与依赖确认
在写代码之前，先确认环境。需要安装GDAL的Python绑定，也就是osgeo模块。安装方式取决于你的系统，用conda的话一般是conda install gdal，用pip的话是pip install GDAL，但pip安装有时候会因为编译问题失败，建议优先用conda。
验证安装是否成功：
from osgeo import gdal, ogr, osr print(gdal.__version__)能打印出版本号就说明环境没问题。我实测下来，GDAL 3.x 版本的API和2.x有一些差异，下面代码基于3.x编写，2.x基本兼容。
3.2 完整源码逐段解析
下面是一个完整的实例，功能是：读取一张分类栅格，把相同值的连通区域提取成多边形，保留像素值作为属性，并计算每个多边形的面积。
from osgeo import gdal, ogr, osr def polygonize_raster(src_path, dst_path, field_name="value", connectedness=4, mask_path=None): """ 将栅格数据多边形化，输出为矢量面文件 参数: src_path: 输入栅格路径 dst_path: 输出矢量路径 field_name: 属性字段名称 connectedness: 连通性，4或8 mask_path: 可选的掩膜栅格路径 """ # 1. 打开源栅格 src_ds = gdal.Open(src_path, gdal.GA_ReadOnly) if src_ds is None: raise RuntimeError(f"无法打开栅格文件: {src_path}") src_band = src_ds.GetRasterBand(1) src_proj = src_ds.GetProjection() src_geotrans = src_ds.GetGeoTransform() # 2. 处理掩膜波段 mask_band = None mask_ds = None if mask_path: mask_ds = gdal.Open(mask_path, gdal.GA_ReadOnly) if mask_ds is None: raise RuntimeError(f"无法打开掩膜文件: {mask_path}") mask_band = mask_ds.GetRasterBand(1) # 3. 创建输出矢量数据源 driver = ogr.GetDriverByName("ESRI Shapefile") if ogr.GetDriverByName("ESRI Shapefile") is None: raise RuntimeError("Shapefile驱动不可用") # 如果文件已存在，先删除 import os if os.path.exists(dst_path): driver.DeleteDataSource(dst_path) dst_ds = driver.CreateDataSource(dst_path) if dst_ds is None: raise RuntimeError(f"无法创建输出文件: {dst_path}") # 4. 创建空间参考 srs = osr.SpatialReference() srs.ImportFromWkt(src_proj) # 5. 创建图层 dst_layer = dst_ds.CreateLayer("polygons", srs, ogr.wkbPolygon) # 6. 创建属性字段 field_defn = ogr.FieldDefn(field_name, ogr.OFTInteger) dst_layer.CreateField(field_defn) # 7. 执行多边形化 options = [] if connectedness == 8: options.append("8CONNECTED=8") result = gdal.Polygonize(src_band, mask_band, dst_layer, 0, options=options, callback=None) if result != 0: raise RuntimeError("Polygonize执行失败") # 8. 计算面积并添加字段 dst_layer.CreateField(ogr.FieldDefn("area", ogr.OFTReal)) area_field_idx = dst_layer.GetLayerDefn().GetFieldIndex("area") dst_layer.ResetReading() for feature in dst_layer: geom = feature.GetGeometryRef() if geom is not None: area = geom.GetArea() feature.SetField(area_field_idx, area) dst_layer.SetFeature(feature) # 9. 清理资源 dst_ds = None src_ds = None if mask_ds: mask_ds = None print(f"多边形化完成，输出文件: {dst_path}") # 调用示例 if __name__ == "__main__": polygonize_raster( src_path="classification.tif", dst_path="output_polygons.shp", field_name="class_id", connectedness=8 )这段代码看起来不长，但每一段都有它的必要性。第1步打开栅格，拿到波段、投影和地理变换信息，这些后面都要用。第2步处理掩膜，如果不需要掩膜就传None。第3步创建输出数据源，这里用的是Shapefile驱动，你也可以换成GeoJSON或者GPKG，只要把驱动名改掉就行。
第4步和第5步创建空间参考和图层，这一步很关键。空间参考必须从源栅格继承过来，否则输出的矢量没有坐标系，后续没法做空间分析。图层的几何类型指定为wkbPolygon，因为Polygonize输出的是面。
第6步创建属性字段，这里只创建了一个整型字段用来存像素值。第7步执行Polygonize，把连通性选项传进去。第8步遍历所有要素，计算面积并写入新字段。最后清理资源，释放文件句柄。
3.3 关键步骤的参数计算与选择依据
关于连通性的选择，我再展开说一下。假设你有一张二值图，前景是1，背景是0。如果前景是一个规则的矩形，4连通和8连通的结果是一样的。但如果前景是一个菱形或者斜线，两者差别就出来了。
举个具体的例子：一个2x2的像素块，四个像素都是1。4连通下，这四个像素互相上下左右相邻，会被合并成一个多边形。但如果是一个斜对角排列的两个像素，比如(0,0)和(1,1)都是1，4连通下它们是两个独立的多边形，8连通下它们是一个多边形。
面积字段的计算也有讲究。geom.GetArea()返回的面积单位取决于你的坐标系。如果坐标系是地理坐标系（经纬度），返回的是平方度，这个值没有实际意义。如果坐标系是投影坐标系（比如UTM），返回的是平方米。所以做面积统计之前，一定要确认坐标系类型。如果是地理坐标系，建议先做投影转换再计算面积。
4. 常见问题与排查技巧实录
4.1 报错类型与对应解决方案
在实际使用中，我遇到过不少报错，整理成表格方便大家速查。
| 报错信息 | 可能原因 | 解决方案 |
|---|---|---|
| TypeError: in method 'Polygonize', argument 1 of type 'GDALRasterBandShadow *' | 传入了Dataset而不是Band | 用GetRasterBand(1)获取波段 |
| RuntimeError: Failed to create field | 字段类型不支持或图层未创建 | 检查图层创建顺序，确认字段类型 |
| 输出文件为空 | 源波段全为0或掩膜全为0 | 检查数据范围，确认非零像素存在 |
| 多边形数量异常多 | 使用了4连通但数据是斜向分布 | 改用8连通 |
| 属性值全部为0 | dstFieldIdx传了-1或字段未创建 | 确认字段索引和字段创建 |
| 内存溢出 | 栅格太大，多边形数量过多 | 分块处理或降低分辨率 |
这个表格里的每一条都是我或者身边同事实际踩过的坑。特别是第一条，新手特别容易犯。GDAL的Python绑定里，很多函数要的是Band，不是Dataset，这个要养成习惯。
4.2 性能优化的几个实操心得
Polygonize 的性能主要受两个因素影响：栅格大小和多边形数量。一张10000x10000的栅格，如果里面有几万个连通区域，处理时间可能会到几分钟甚至更久。
第一个优化思路是分块处理。把大栅格切成小块，逐块Polygonize，最后合并。但这样有个问题：跨块的多边形会被切断。解决办法是块之间留一定的重叠区域，合并的时候再去重。这个方案实现起来比较复杂，适合对性能要求极高的场景。
第二个思路是先用gdal.SieveFilter做小斑块过滤，把面积小于阈值的连通区域先去掉，减少多边形数量。这个操作能显著提升后续Polygonize的速度，而且能去掉很多噪点。
# 先做小斑块过滤 gdal.SieveFilter(src_band, None, src_band, threshold=10, connectedness=8)第三个思路是降低分辨率。如果业务上不需要像素级的精度，可以先用gdal.Warp做重采样，把栅格缩小，再Polygonize。这样多边形数量会大幅减少，速度提升明显。
4.3 输出格式选择的经验之谈
输出格式的选择也影响使用体验。Shapefile是最常见的格式，但它有几个限制：字段名不能超过10个字符，单个文件不能超过2GB，不支持UTF-8编码的中文字段。如果你的属性字段名比较长，或者需要存中文，建议用GeoJSON或者GPKG。
GeoJSON的好处是文本格式，方便查看和传输，但文件体积比较大。GPKG是SQLite数据库格式，支持大数据量，字段名和编码都没有限制，我个人最推荐这个格式。切换格式只需要改驱动名：
driver = ogr.GetDriverByName("GPKG") dst_ds = driver.CreateDataSource("output.gpkg")注意GPKG的图层创建方式略有不同，需要指定图层名和选项，但整体逻辑是一样的。
5. 进阶用法与场景延展
5.1 多波段栅格的批量处理
实际项目中，经常遇到多波段栅格，每个波段都需要单独Polygonize。这时候可以写一个循环，遍历所有波段，每个波段输出一个独立的矢量文件，或者输出到同一个数据源的不同图层。
def polygonize_all_bands(src_path, dst_dir): src_ds = gdal.Open(src_path) band_count = src_ds.RasterCount for i in range(1, band_count + 1): band = src_ds.GetRasterBand(i) dst_path = f"{dst_dir}/band_{i}.shp" # 这里复用前面的polygonize逻辑 # 注意每次都要重新创建输出数据源这个用法在时间序列分析里很常见，比如每个波段代表一个月份的NDVI分类结果，需要逐月提取多边形做变化分析。
5.2 结合属性过滤做条件提取
有时候我们不需要提取所有像素，只想提取特定值的区域。这时候有两种做法：一种是用掩膜波段，把不需要的像素设为0；另一种是先Polygonize全部，再用属性查询过滤。
第一种做法效率更高，因为Polygonize处理的像素少了。构造掩膜的方法很简单：
import numpy as np from osgeo import gdal # 读取源数据 ds = gdal.Open("input.tif") band = ds.GetRasterBand(1) arr = band.ReadAsArray() # 构造掩膜：只保留值为1的像素 mask_arr = np.where(arr == 1, 1, 0).astype(np.uint8) # 写入掩膜文件 driver = gdal.GetDriverByName("GTiff") mask_ds = driver.Create("mask.tif", ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte) mask_ds.SetGeoTransform(ds.GetGeoTransform()) mask_ds.SetProjection(ds.GetProjection()) mask_ds.GetRasterBand(1).WriteArray(mask_arr) mask_ds = None然后用这个掩膜文件去调用Polygonize，就只有值为1的区域会被提取出来。
5.3 多边形后处理：简化与平滑
Polygonize输出的多边形是像素级的，边缘呈锯齿状。如果用于制图或者展示，通常需要做简化和平滑。GDAL提供了gdal.Simplify相关的功能，但更常用的是OGR的几何处理方法。
# 对每个多边形做简化 for feature in dst_layer: geom = feature.GetGeometryRef() simplified = geom.Simplify(tolerance=0.5) feature.SetGeometry(simplified) dst_layer.SetFeature(feature)Simplify的tolerance参数控制简化程度，值越大越简化。这个操作会改变多边形的形状和面积，所以如果面积精度要求高，要谨慎使用。我一般会在简化之后重新计算面积，保证属性数据的准确性。
6. 写在最后的一些个人体会
Polygonize这个函数，看起来简单，但真正用好需要理解栅格和矢量的转换逻辑。我刚开始用的时候，总觉得参数太多记不住，后来用多了发现，核心就是三个东西：源波段、输出图层、字段索引。把这三个搞明白，其他都是细节。
还有一个体会是，GDAL的文档虽然全，但很多细节藏在C++ API的说明里，Python绑定的文档相对简略。遇到不确定的参数，直接去看GDAL源码或者用help(gdal.Polygonize)看签名，比在网上搜半天管用。
最后分享一个小技巧：如果你的数据量特别大，Polygonize跑一次要很久，建议先用小范围的数据测试参数，确认连通性和字段设置没问题之后，再跑全量数据。这样能避免跑了一半发现参数错了要重来的尴尬。我自己就吃过这个亏，一张几万乘几万的栅格跑了二十分钟，结果发现字段索引传错了，只能重跑。
