GDAL/Python：如何修复多边形自相交错误以解决栅格裁剪失败问题

问题描述 #

使用GDAL的gdal.Warp函数通过Shapefile裁剪栅格时执行失败，报错显示裁剪多边形无效，同时提示存在环自相交问题。

原代码 #

from osgeo import ogr, osr, gdal
shpfile="a.shp"
input_file="area.tiff"
output_file = "subarea.tiff"
outTile = gdal.Warp(output_file, input_file, cutlineDSName=shpfile, format="GTiff")
OutTile = None

报错信息（翻译后） #

警告1：在或靠近点114.56538630983039 30.263472880169331处存在环自相交
outTile = gdal.Warp(output_file, input_file, cutlineDSName=shpfile, format="GTiff")
文件 "C:\ProgramData\Anaconda3\envs\base1\lib\site-packages\osgeo\gdal.py", 第625行, 在Warp函数中
返回 wrapper_GDALWarpDestName(destNameOrDestDS, srcDSTab, opts, callback, callback_data)
文件 "C:\ProgramData\Anaconda3\envs\base1\lib\site-packages\osgeo\gdal.py", 第3410行, 在wrapper_GDALWarpDestName函数中
返回 _gdal.wrapper_GDALWarpDestName(*args)

错误1：裁剪多边形无效。

原因分析 #

报错核心是Shapefile中的多边形存在自相交问题，GDAL无法识别无效几何图形完成裁剪操作。

解决方案 #

方案1：修复Shapefile的自相交几何（推荐） #

确保裁剪边界有效是最稳妥的解决方式，可通过以下途径修复：

• QGIS可视化修复：打开Shapefile，选择处理工具→矢量几何→修复几何图形工具，一键修复自相交问题。
• OGR命令行修复：

  ogr2ogr -f "ESRI Shapefile" fixed_a.shp a.shp -dialect sqlite -sql "SELECT ST_MakeValid(geometry) AS geometry FROM a"
  

• Python代码批量修复：

  from osgeo import ogr
  import os
  
  # 打开原始Shapefile
  src_ds = ogr.Open("a.shp", 0)
  src_layer = src_ds.GetLayer()
  
  # 创建修复后的Shapefile
  driver = ogr.GetDriverByName("ESRI Shapefile")
  if os.path.exists("fixed_a.shp"):
      driver.DeleteDataSource("fixed_a.shp")
  dst_ds = driver.CreateDataSource("fixed_a.shp")
  dst_layer = dst_ds.CreateLayer("fixed_a", src_layer.GetSpatialRef(), ogr.wkbPolygon)
  
  # 复制属性字段
  src_layer_def = src_layer.GetLayerDefn()
  for i in range(src_layer_def.GetFieldCount()):
      field_def = src_layer_def.GetFieldDefn(i)
      dst_layer.CreateField(field_def)
  
  # 逐个修复几何并写入新文件
  for feature in src_layer:
      geom = feature.GetGeometryRef()
      if not geom.IsValid():
          geom = geom.MakeValid()
          # 处理修复后可能出现的多几何类型
          if geom.GetGeometryType() == ogr.wkbMultiPolygon:
              for poly in geom:
                  new_feature = ogr.Feature(dst_layer.GetLayerDefn())
                  new_feature.SetGeometry(poly)
                  for i in range(src_layer_def.GetFieldCount()):
                      new_feature.SetField(i, feature.GetField(i))
                  dst_layer.CreateFeature(new_feature)
          else:
              feature.SetGeometry(geom)
              dst_layer.CreateFeature(feature)
      else:
          dst_layer.CreateFeature(feature)
  
  # 释放资源
  src_ds = None
  dst_ds = None
  

  修复完成后，将代码中的shpfile路径改为fixed_a.shp重新执行即可。

方案2：修改GDAL Warp参数忽略无效几何（临时应急） #

若暂时无法修复Shapefile，可添加参数让GDAL跳过有效性检查，但可能导致裁剪结果不准确：

from osgeo import gdal
shpfile="a.shp"
input_file="area.tiff"
output_file = "subarea.tiff"
# 关闭裁剪多边形有效性检查
outTile = gdal.Warp(output_file, input_file, 
                    cutlineDSName=shpfile, 
                    format="GTiff",
                    options=["CUTLINE_VALIDITY_CHECK=NO"])
outTile = None

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