原標(biāo)題：shp文件自相交處理的方法

今天基于GDAL使用shp文件對(duì)柵格影像進(jìn)行裁剪時(shí)出現(xiàn)了下面的問(wèn)題，提示多邊形自相交了

Warning1: RingSelf-intersectionatornearpoint112 .4866642030000334 .830899357000078

ERROR1: Cutlinepolygonisinvalid.

很多人的第一反應(yīng)是使用ArcGIS進(jìn)行拓?fù)錂z查，或使用ArcToolBox里的修復(fù)幾何。確實(shí)我的第一反應(yīng)也是去做這些東西。但是結(jié)果卻沒(méi)有檢查出任務(wù)拓?fù)溴e(cuò)誤，幾何修復(fù)也沒(méi)有檢查出問(wèn)題。

使用ArcGIS檢查不到的原因

強(qiáng)大的ArcGIS居然檢查不到，最終找到了這個(gè)原因。

在ArcGIS 中無(wú)論是拓?fù)洹hapefile文件、還是個(gè)人地理數(shù)據(jù)庫(kù)都是設(shè)置有容差的，小于這個(gè)容差的自相交，都是無(wú)法檢測(cè)到的。

解決方案

查閱了很多資料，最終整理了如下的解決方案。

1.使用PostGIS將shape文件導(dǎo)入Postgresql數(shù)據(jù)庫(kù)，記得導(dǎo)入的時(shí)候要勾選下面的選項(xiàng)。

2.從表里提取出自相交的多邊形

CREATETABLEtemp1 asselect* fromm2 whereST_IsValid(geom) = false

3.刪除原表中的自相交圖形

delete fromm2 whereST_IsValid( geom)= false

4.修復(fù)多邊形

updatetemp1 setgeom =ST_Buffer(geom, 0.0)

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

5.修復(fù)完的數(shù)據(jù)恢復(fù)到原來(lái)的表

insertintom2 select* fromtemp1

6.最后通過(guò)PostGIS插件導(dǎo)出shp文件即可

結(jié)果檢測(cè)

使用gdal對(duì)結(jié)果進(jìn)行檢測(cè)

fromosgeo importogr

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

# # # 注冊(cè)所有的驅(qū)動(dòng)ogr.RegisterAll

defcheck_shp:# 打開(kāi)數(shù)據(jù)ds = ogr.Open(shpFile, 0)ifds isNone:print( "打開(kāi)文件【%s】失敗！", shpFile)returnprint( "打開(kāi)文件【%s】成功！", shpFile)# 獲取該數(shù)據(jù)源中的圖層個(gè)數(shù)，一般shp數(shù)據(jù)圖層只有一個(gè)，如果是mdb、dxf等圖層就會(huì)有多個(gè)m_layer_count = ds.GetLayerCountm_layer = ds.GetLayerByIndex( 0)ifm_layer isNone:print( "獲取第%d個(gè)圖層失敗！n", 0)return# 對(duì)圖層進(jìn)行初始化，如果對(duì)圖層進(jìn)行了過(guò)濾操作，執(zhí)行這句后，之前的過(guò)濾全部清空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( "無(wú)效多邊形共"+ str(count) + "個(gè)")

check_shp

運(yùn)行結(jié)果

本文作者：GIS無(wú)情老博士

總結(jié)

以上是生活随笔為你收集整理的java shp求相交面积_shp文件自相交处理的方法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。