![[草根哥]学习资料大杂烩-自学逆袭风水宝地](http://www.cgg6.com/wp-content/uploads/2023/01/zhengchang-logo.png){kind=logo}
草根哥今天基于GDAL使用shp文件对栅格影像进行裁剪时出现了下面的问题，提示多边形自相交了
Warning 1: Ring Self-intersection at or near point 112.48666420300003 34.830899357000078ERROR 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()
运行结果
评论前必须登录！
注册