如何将聚类后的农业用地TIF文件按需求多边形化？

这是一张已完成聚类操作的农业用地影像，我希望将图中每个矩形区域转为多边形，不受像素值影响，但用Python代码和QGIS尝试后没达到预期效果。

我使用的多边形化代码如下：

input_path = "/home/saeed/Desktop/SatPlat_Goats/39SUR/mask/mask.tif"
poliginize_shape_path = "Polgonize"
polygonize_cmd = f"gdal_polygonize.py -8 {input_path} -b 1 {poliginize_shape_path}"
os.system(polygonize_cmd)

# Open the input shapefile
in_ds = ogr.Open('/home/saeed/Desktop/SatPlat_Goats/39SUR/polygon/out.shp')
print(in_ds)
# Get the input layer
in_layer = in_ds.GetLayer()

# Create a new shapefile to hold the filtered polygons
driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = driver.CreateDataSource("39SUR/polygon/poly100000")
out_layer = out_ds.CreateLayer('polygons', geom_type=ogr.wkbPolygon)
# Copy the input layer's fields to the output layer
in_layer_defn = in_layer.GetLayerDefn()
for i in range(in_layer_defn.GetFieldCount()):
    field_defn = in_layer_defn.GetFieldDefn(i)
    out_layer.CreateField(field_defn)
# Loop through the input layer's features
for in_feat in in_layer:
    # Get the feature's geometry and calculate its area
    geom = in_feat.GetGeometryRef()
    area = geom.GetArea()
    area = area * 10000000000
# print(area)

    # Check if the feature's area is greater than or equal to 10 square meters
    if area >= 100000:

        # Create a new feature in the output layer and set its geometry and attributes
        out_feat = ogr.Feature(out_layer.GetLayerDefn())
        out_feat.SetGeometry(geom.Clone())
        for i in range(out_feat.GetFieldCount()):
            out_feat.SetField(out_layer.GetLayerDefn().GetFieldDefn(i).GetNameRef(), in_feat.GetField(i))

        # Add the feature to the output layer
        out_layer.CreateFeature(out_feat)
print("3333333333333333333333")
# Close the shapefiles
in_ds = None
out_ds = None

当前代码能实现多边形化，但存在问题：部分农业用地因NDVI值偏低，会和主体农地区域被分割开。

核心思路是先消除影像内部的小异值区域，再进行多边形化，避免像素值差异导致的区域分割，以下是三种可行方法：

方法1：GDAL形态学预处理影像 #

先对原始mask影像做闭运算（先膨胀后腐蚀），填充小的孔洞或异值区域，再执行多边形化：

import os

# 1. 对原始影像做闭运算填充异值区域
input_tif = "/home/saeed/Desktop/SatPlat_Goats/39SUR/mask/mask.tif"
processed_tif = "/home/saeed/Desktop/SatPlat_Goats/39SUR/mask/mask_closed.tif"

# 用GDAL填充工具处理小范围异值（md参数为处理的最大距离，可根据实际调整）
fill_cmd = f"gdal_fillnodata.py -md 3 {input_tif} {processed_tif}"
os.system(fill_cmd)

# 2. 用处理后的影像执行多边形化
poliginize_shape_path = "Polgonize"
polygonize_cmd = f"gdal_polygonize.py -8 {processed_tif} -b 1 {poliginize_shape_path}"
os.system(polygonize_cmd)

# 后续面积筛选代码保持不变

方法2：QGIS可视化操作 #

1. 加载原始mask影像到QGIS
2. 打开工具箱 → 形态学处理 → 闭运算，选择合适的结构元素大小（如3x3或5x5，根据异值区域尺寸调整）
3. 对处理后的影像执行矢量 → 栅格转矢量（多边形）
4. 最后用按属性选择或提取矢量工具筛选面积符合要求的多边形

方法3：多边形化后合并相邻区域 #

如果不想预处理影像，可在多边形化后，将相邻的同地块小多边形合并：

# 加载多边形化后的图层
in_ds = ogr.Open('/home/saeed/Desktop/SatPlat_Goats/39SUR/polygon/out.shp')
in_layer = in_ds.GetLayer()

# 创建输出图层
driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = driver.CreateDataSource("39SUR/polygon/merged_poly.shp")
out_layer = out_ds.CreateLayer('merged_polygons', geom_type=ogr.wkbPolygon)

# 合并所有相邻多边形
merged_geom = None
for feat in in_layer:
    geom = feat.GetGeometryRef()
    if merged_geom is None:
        merged_geom = geom.Clone()
    else:
        merged_geom = merged_geom.Union(geom)

# 将合并后的几何写入输出图层
out_feat = ogr.Feature(out_layer.GetLayerDefn())
out_feat.SetGeometry(merged_geom)
out_layer.CreateFeature(out_feat)

# 关闭数据源
in_ds = None
out_ds = None

注：此方法适合相邻且属性相近的多边形，可根据实际需求添加合并条件（如属性匹配）。

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