欢迎光临草根哥
www.cgg6.com

shp文件自相交处理的方法

今天基于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()

 

运行结果

 

赞(0)
版权声明:
文章名称:《shp文件自相交处理的方法》
文章链接:http://www.cgg6.com/967.html
声明:文章版权归本站自创作者所有,未经允许不得转载

评论 抢沙发

评论前必须登录!