shape文件自相交问题（Warning 1: Ring Self-intersection at or near point 112.48666420300003 34.8308993570000）

问题

今天在使用GDAL模块时，出现了下面的报错信息，就是新做的矢量文件没有办法在GDAL中使用。
Warning 1: Ring Self-intersection at or near point 112.48666420300003 34.830899357000078

多边形边界自相交是什么意思？

矢量文件自相交意味着该矢量文件中的多边形边界存在重叠或交叉的情况。在地理信息系统（GIS）中，矢量数据通常由点、线、面等几何元素组成，用于表示地理空间中的实体，如行政边界、土地利用类型等。

当一个矢量文件中的多边形自相交时，这通常意味着多边形的边界线与其他部分重叠或交叉，导致多边形内部出现不规则的分割或重叠区域。这种情况在GIS处理中是不合法的，因为它可能导致数据处理和分析时的错误。

例如，在剪裁栅格影像时，如果使用的矢量多边形文件存在自相交，可能会导致剪裁结果不准确，因为软件无法确定应该如何处理这些自相交的区域。解决这个问题通常需要对矢量文件进行编辑，消除自相交的部分，确保多边形的边界是封闭且不自相交的。这可以通过GIS软件中的编辑工具来完成，如使用拓扑编辑功能来修复这些错误。

为什么会出现自相交问题？

这个矢量文件我是在Arcgis中栅格转面要素再转shape做的，可能由以下原因造成：

• 复杂图形的简化：在对矢量数据进行简化或一般化（减少图形细节以减小文件大小或显示简化图形）时可能会引入自相交，尤其是如果简化算法不够鲁棒或未正确处理图形边缘。

• 导入或转换数据时的错误：在从一个格式转换到另一个格式，或者从一个GIS系统导入到另一个GIS系统时，坐标精度的变化或者处理的不当可能会引入自相交。

• 自动化处理的副作用：使用某些自动化工具或脚本来生成或修改矢量数据时，如果工具的算法没有正确处理几何约束条件，可能会创建出自相交的几何形状。

解决方法

参考了‘KKK开讲啦’博主提供的python+ogr解决shp文件自相交的问题
解决方案

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)

核心原理

在地理信息系统（GIS）和数据处理中，“Union” 通常指的是一种空间操作，它将两个或多个空间数据集（如多边形）合并在一起，生成一个新的数据集，其中包含了所有输入数据集的空间范围和特征。这种操作在GIS软件中非常常见，用于分析和处理地理空间数据。

在GDAL/OGR库中，ogr模块提供了一种方式来执行这种空间操作。例如，你可以使用ogr.Geometry对象的Union方法来合并两个几何对象。

以下是一个使用ogr模块执行Union操作的简单示例：

from osgeo import ogr

# 创建两个多边形几何对象
poly1 = ogr.CreateGeometryFromWkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))")
poly2 = ogr.CreateGeometryFromWkt("POLYGON((0.5 0.5, 0.5 1.5, 1.5 1.5, 1.5 0.5, 0.5 0.5))")

# 执行Union操作
union_geom = poly1.Union(poly2)

# 打印Union后的几何对象的WKT表示
print(union_geom.ExportToWkt())

在这个例子中，我们创建了两个多边形几何对象，然后使用Union方法将它们合并。结果是一个新的几何对象，它包含了两个原始多边形的所有空间范围。

在更复杂的GIS分析中，Union操作可以用于合并多个图层或数据集，以便进行进一步的空间分析或数据处理。例如，你可以使用Union来合并多个行政边界图层，以创建一个包含所有行政区域的新图层。

请注意，Union操作可能会导致数据集中的冗余或重叠区域，因此在处理和分析结果时需要考虑这些因素。

结果

最终在GDAL模块中能使用新的处理过的矢量文件了，就是有个问题发现回到Arcgis中不能再用了。