统计
  • 文章总数:710 篇
  • 评论总数:0 条
  • 分类总数:14 个
  • 最后更新:7月9日

shp文件自相交处理的方法

本文阅读 3 分钟
首页 地理信息 正文

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

运行结果

 

本文来自投稿,不代表本站立场,如若转载,请注明出处:
用ArcGIS desktop几种制作“热力图”的方法
« 上一篇 03-13
CAD制图的几个实用技巧!
下一篇 » 03-13

作者信息

作者有点忙,还没写简介
TA的最新作品
    请设置要调用的作者ID

动态快讯

    请配置好页面缩略名选项

热门文章

标签TAG