原标题：shp文件自相交处理的方法

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

Warning1: RingSelf-intersectionatornearpoint112 .4866642030000334 .830899357000078

ERROR1: Cutlinepolygonisinvalid.

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

使用ArcGIS检查不到的原因

强大的ArcGIS居然检查不到，最终找到了这个原因。

在ArcGIS 中无论是拓扑、shapefile文件、还是个人地理数据库都是设置有容差的，小于这个容差的自相交，都是无法检测到的。

解决方案

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

1.使用PostGIS将shape文件导入Postgresql数据库，记得导入的时候要勾选下面的选项。

2.从表里提取出自相交的多边形

CREATETABLEtemp1 asselect* fromm2 whereST_IsValid(geom) = false

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

delete fromm2 whereST_IsValid( geom)= false

4.修复多边形

updatetemp1 setgeom =ST_Buffer(geom, 0.0)

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

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

insertintom2 select* fromtemp1

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

结果检测

使用gdal对结果进行检测

fromosgeo importogr

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

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

defcheck_shp:# 打开数据ds = ogr.Open(shpFile, 0)ifds isNone:print( "打开文件【%s】失败！", shpFile)returnprint( "打开文件【%s】成功！", shpFile)# 获取该数据源中的图层个数，一般shp数据图层只有一个，如果是mdb、dxf等图层就会有多个m_layer_count = ds.GetLayerCountm_layer = ds.GetLayerByIndex( 0)ifm_layer isNone:print( "获取第%d个图层失败！n", 0)return# 对图层进行初始化，如果对图层进行了过滤操作，执行这句后，之前的过滤全部清空m_layer.ResetReadingcount = 0m_feature = m_layer.GetNextFeaturewhilem_feature isnotNone:o_geometry = m_feature.GetGeometryRefifnotogr.Geometry.IsValid(o_geometry):print(m_feature.GetFID)count = count + 1

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

check_shp

运行结果

本文作者：GIS无情老博士