shp文件自相交处理 |
您所在的位置：网站首页 › arcgis拓扑自相交 › shp文件自相交处理 |
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 中无论是拓扑、shapefile文件、还是个人地理数据库都是设置有容差的，小于这个容差的自相交，都是无法检测到的。  解决方案
查阅了很多资料，最终整理了如下的解决方案。 使用PostGIS将shape文件导入Postgresql数据库，记得导入的时候要勾选下面的选项。 从表里提取出自相交的多边形 CREATE TABLE temp1 as select * from m2 where ST_IsValid(geom) = false 删除原表中的自相交图形 delete from m2 where ST_IsValid(geom) = false
修复多边形 update temp1 set geom =ST_Buffer(geom, 0.0) -- update temp1 set geom =ST_MakeValid(geom) 也可以
修复完的数据恢复到原来的表 insert into m2 select * from temp1
最后通过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()运行结果
|
【本文地址】
今日新闻 |
推荐新闻 |