shp文件自相交处理

今天基于GDAL使用shp文件对栅格影像进行裁剪时出现了下面的问题，提示多边形自相交了

Warning 1: Ring Self-intersection at or near point 112.48666420300003 34.830899357000078
ERROR 1: Cutline polygon is invalid.

很多人的第一反应是使用ArcGIS进行**拓扑检查**，或使用ArcToolBox里的**修复几何**。确实我的第一反应也是去做这些东西。但是结果缺没有检查出任务拓扑错误，几何修复也没有检查出问题。

使用ArcGIS检查不到的原因

强大的ArcGIS居然检查不到，最终找到了这个原因。
在ArcGIS 中无论是拓扑、shapefile文件、还是个人地理数据库都是设置有容差的，小于这个容差的自相交，都是无法检测到的。

解决方案

查阅了很多资料，最终整理了如下的解决方案。

1. 使用PostGIS将shape文件导入Postgresql数据库，记得导入的时候要勾选下面的选项。
2. 从表里提取出自相交的多边形

   CREATE TABLE temp1 as select * from m2 where ST_IsValid(geom) = false

    

3. 删除原表中的自相交图形

   delete from  m2 where ST_IsValid(geom) = false

    

4. 修复多边形

   update temp1 set geom =ST_Buffer(geom, 0.0) 
   -- update temp1 set geom =ST_MakeValid(geom) 也可以

    

5. 修复完的数据恢复到原来的表

   insert into m2 select * from temp1

    

6. 最后通过PostGIS插件导出shp文件即可

 结果检测

使用gdal对结果进行检测
 

from osgeo import ogr

shpFile = 'F:/m2.shp'  # 裁剪矩形

# # # 注册所有的驱动
ogr.RegisterAll()


def check_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
    # 对图层进行初始化，如果对图层进行了过滤操作，执行这句后，之前的过滤全部清空
    m_layer.ResetReading()
    count = 0
    m_feature = m_layer.GetNextFeature()
    while m_feature is not None:
        o_geometry = m_feature.GetGeometryRef()
        if not ogr.Geometry.IsValid(o_geometry):
            print(m_feature.GetFID())
            count = count + 1

        m_feature = m_layer.GetNextFeature()
    print("无效多边形共" + str(count) + "个")


check_shp()

运行结果