python+ogr解决shp文件自相交的问题

Warning1: RingSelf-intersectionatornearpoint112 .4866642030000334 .830899357000078 ERROR1: Cutlinepolygonisinvalid.

使用矢量shp文件剪裁栅格影像时会出现上述问题。提示shp文件自相交。搜集了大量资料，多数回答都是基于PostGIS进行处理。在下载安装过程中很复杂，本人试过多次也没有成功。这里给出一种使用python+ogr函数库的解决方法。

不废话，直接上代码。

from osgeo import ogr, gdal
import os

gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
gdal.SetConfigOption("SHAPE_ENCODING", "GBK")


# 思路:https://github.com/whosonfirst-data/whosonfirst-data/issues/1071
# Solve shapefile self-intersection
def s_shp_si(shpFile):
    out_shp = shpFile[:-4] + '_ssi' + '.shp'
    print(out_shp)
    # 打开数据
    ds = ogr.Open(shpFile, 0)
    if ds is None:
        print("打开文件 %s 失败！" % shpFile)
        return
    print("打开文件%s成功！" % shpFile)
    # 获取该数据源中的图层个数，一般shp数据图层只有一个，如果是mdb、dxf等图层就会有多个
    m_layer_count = ds.GetLayerCount()
    m_layer = ds.GetLayerByIndex(0)
    if m_layer is None:
        print("获取第%d个图层失败！\n", 0)
        return
    # 创建输出文件
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(out_shp):
        driver.DeleteDataSource(out_shp)
    outds = driver.CreateDataSource(out_shp)
    outlayer = outds.CreateLayer(out_shp[:-4], geom_type=ogr.wkbPolygon, options=["ENCODING=GBK"])
    # 获取输出层的要素定义
    outLayerDefn = outlayer.GetLayerDefn()
    # 对图层进行初始化，如果对图层进行了过滤操作，执行这句后，之前的过滤全部清空
    m_layer.ResetReading()
    # 获取投影
    prosrs = m_layer.GetSpatialRef()
    # 添加字段
    inLayerDefn = m_layer.GetLayerDefn()
    for i in range(0, inLayerDefn.GetFieldCount()):
        fieldDefn = inLayerDefn.GetFieldDefn(i)
        outlayer.CreateField(fieldDefn)

    # loop through the input features
    m_feature = m_layer.GetNextFeature()
    while m_feature:
        # print(m_feature)
        o_geometry = m_feature.GetGeometryRef()
        # 关键，合并几何
        o_geometry = o_geometry.Union(o_geometry)
        outfeature = ogr.Feature(outLayerDefn)
        outfeature.SetGeometry(o_geometry)
        # 遍历每个要素的字段，并设置字段属性
        for i in range(0, outLayerDefn.GetFieldCount()):
            print(outLayerDefn.GetFieldDefn(i).GetNameRef())
            outfeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), m_feature.GetField(i))
        outlayer.CreateFeature(outfeature)
        # dereference the features and get the next input feature
        outfeature = None
        m_feature = m_layer.GetNextFeature()

    ds.Destroy()
    outds.Destroy()
    # 保存投影文件

    prjfile = (out_shp[:-4] + '.prj')
    fn = open(prjfile, 'w')
    fn.write(prosrs.ExportToWkt())
    fn.close()


shpFile = r'E:/hhly/hhly.shp'  # 裁剪矩形
# # # 注册所有的驱动
ogr.RegisterAll()
s_shp_si(shpFile)

参考链接：
https://blog.csdn.net/gisuuser/article/details/103293038
https://github.com/whosonfirst-data/whosonfirst-data/issues/1071
https://blog.csdn.net/weixin_42924891/article/details/85469881