zl程序教程

您现在的位置是:首页 >  其它

当前栏目

geopandas拓扑检查(任意两个几何体不相交)

两个 检查 任意 拓扑 相交
2023-09-14 09:16:18 时间

之前使用rtree来进行拓扑检查,主要是使用GeoDataFramesindex来实现的,但是由于某种未知的原因,rtree总是会莫名奇妙的运行错误,且没有任何错误提示,因此博主自己实现了一版拓扑检查,思想和之前其实差不多:先得到所有几何体的外接矩形框,然后根据外接矩形框判断是否有可能相交(矩形框相交,几何体才有相交的可能性),通过外接矩形框的筛选之后,再对可能相交的几何体进行一个精确的相交判断,代码实现如下:

import geopandas
import copy
import matplotlib.pyplot as plt
import time

gdf=geopandas.read_file("D:\\problem.shp")
gdf=gdf[gdf.boundary.is_simple]
gdf.reset_index(drop=True,inplace=True)

import shapely.geometry as geo

def RectangleIntersection(r1,r2):      #判断两矩形是否相交 r[xmin,ymin,xmax,ymax]
    x1min,y1min,x1max,y1max=r1
    x2min,y2min,x2max,y2max=r2
    w1=x1max-x1min
    h1=y1max-y1min
    w2=x2max-x2min
    h2=y2max-y2min
    xc1=(x1min+x1max)/2
    yc1=(y1min+y1max)/2
    xc2=(x2max+x2min)/2
    yc2=(y2max+y2min)/2
    if abs(xc1-xc2) <= (w1+w2)/2 and abs(yc1-yc2) <= (h1+h2)/2:
        return True
    else:
        return False

def topology_check1(gdf):
    drop_list=[]
    n=len(gdf)
    for i in range(len(gdf)):
        polygon=gdf.loc[i,'geometry']
        polygons=gdf.loc[i+1:n,'geometry']
        intersect_result=polygons.intersects(polygon)
        if intersect_result.any():
            drop_list.append(i)
    return drop_list

def topology_check2(gdf):
    bounds=gdf.bounds
    rectangles=[]
    for i in range(len(gdf)):
        rectangles.append(list(bounds.loc[i,:]))
    #
    drop_list=[]
    for i in range(len(gdf)):
        polygon=gdf.loc[i,'geometry']
        polygon_bounds=rectangles[i]
        #根据矩形框查找可能相交的几何体
        possible_intersects_idx=[]
        for j in range(i+1,len(gdf)):
            if RectangleIntersection(polygon_bounds,rectangles[j]):
                possible_intersects_idx.append(j)
        print(i,possible_intersects_idx)
        #从可能相交的几何体里精确查询
        find_flag=0
        for index in possible_intersects_idx:
            if gdf.loc[index,'geometry'].intersects(polygon):
                find_flag=1
                break
        if find_flag==1:
            drop_list.append(i)
    return drop_list


start=time.time()
drop_list=topology_check2(gdf)
end=time.time()
print(drop_list)
print(end-start)
gdf.drop(axis=0,index=drop_list,inplace=True)
gdf.reset_index(drop=True,inplace=True)
gdf.to_file("D:\\solved.shp")


这里对矩形相交的判断可参考博客链接

更新

import geopandas
import shapely.geometry as geo
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

# polygon0=geo.Polygon([(-1,0),(0,0),(0,1),(-1,1)])
# polygon1=geo.Polygon([(0,0),(8,0),(8,2),(0,2)])
# polygon2=geo.Polygon([(3,-6),(5,-6),(5,6),(3,6)])
# polygon3=geo.Polygon([(2,-2),(6,-2),(4,10)])
#
# geopandas.GeoSeries([polygon1,polygon0,polygon2,polygon3]).boundary.plot()
# plt.show()
# df=geopandas.overlay(geopandas.GeoDataFrame(geometry=[polygon1]),geopandas.GeoDataFrame(geometry=[polygon0,polygon2,polygon3]),how='difference')
# df=df.explode()
# df.reset_index(drop=True,inplace=True)
# df.plot()
# plt.show()
# p1=df.loc[0,'geometry']
# p2=polygon3.buffer(distance=-0.001,resolution=1)
# geopandas.GeoSeries([p2,p1]).boundary.plot()
# plt.show()
# print(p1.intersects(p2))
# print(p1.touches(p2))
# print(p1.overlaps(p2))
# print(p1)
# print(p2)
# print(p1.intersection(p2))
# print(p1.overlaps(p2))
# print(p1.touches(p2))



gdf=geopandas.read_file("F:\\shapefile\\full\\Detroit.shp")
gdf=gdf[gdf.boundary.is_simple]
gdf.set_crs(crs=4326,inplace=True)
gdf.to_crs(crs=3395,inplace=True)
gdf.reset_index(drop=True,inplace=True)


import copy
def check_overlap(gdf):
    gdf_copy=copy.deepcopy(gdf)
    spatial_index = gdf_copy.sindex
    #建立R树索引.
    drop_list=[]
    for i in range(len(gdf_copy)):
        polygon=gdf_copy.loc[i,'geometry']
        possible_matches_index = list(spatial_index.intersection(polygon.bounds))
        possible_matches = gdf_copy.iloc[possible_matches_index]
        #通过R树进行简单包含查询
        precise_matches_intersects = possible_matches[possible_matches.intersects(polygon)]
        precise_matches_touches = possible_matches[possible_matches.touches(polygon)]
        #将简单查询结果进行进行匹配,得到精确查询结果
        spatial_index.delete(i, polygon.bounds)
        if (len(precise_matches_intersects)-len(precise_matches_touches))>1:
            drop_list.append(i)
    return drop_list

def RectangleIntersection(r1,r2):      #判断两矩形是否相交 r[xmin,ymin,xmax,ymax]
    x1min,y1min,x1max,y1max=r1
    x2min,y2min,x2max,y2max=r2
    w1=x1max-x1min
    h1=y1max-y1min
    w2=x2max-x2min
    h2=y2max-y2min
    xc1=(x1min+x1max)/2
    yc1=(y1min+y1max)/2
    xc2=(x2max+x2min)/2
    yc2=(y2max+y2min)/2
    if abs(xc1-xc2) <= (w1+w2)/2 and abs(yc1-yc2) <= (h1+h2)/2:
        return True
    else:
        return False

def topology_check3(gdf):
    bounds=gdf.bounds
    rectangles=[]
    for i in range(len(gdf)):
        rectangles.append(list(bounds.loc[i,:]))            #得到所有几何体的MBR
    #
    drop_list=[]
    for i in range(len(gdf)):
        polygon=gdf.loc[i,'geometry']
        polygon_bounds=rectangles[i]
        #根据矩形框查找可能相交的几何体
        possible_intersects_idx=[]
        for j in range(i+1,len(gdf)):
            if RectangleIntersection(polygon_bounds,rectangles[j]):
                possible_intersects_idx.append(j)

        #从可能相交的几何体里精确查询
        find_flag=0
        for index in possible_intersects_idx:
            if gdf.loc[index,'geometry'].intersects(polygon):
                find_flag=1
                print(i,possible_intersects_idx)
                break
        #difference
        if find_flag==1:
            df1=geopandas.GeoDataFrame(geometry=[polygon])
            if len(possible_intersects_idx)>1:
                df2=geopandas.GeoDataFrame(geometry=gdf.loc[possible_intersects_idx,'geometry'])
            else:
                df2=geopandas.GeoDataFrame(geometry=[gdf.loc[possible_intersects_idx[0],'geometry'] for _ in range(2)])
            new_polygon=geopandas.overlay(df1,df2,how='difference')
            if len(new_polygon)==0:
                drop_list.append(i)
            elif len(new_polygon)==1:
                gdf.geometry[i]=new_polygon.loc[0,'geometry']
                gdf.geometry[i]=gdf.geometry[i].buffer(distance=-0.01)
                for index in possible_intersects_idx:
                    if gdf.loc[index, 'geometry'].intersects(gdf.geometry[i]):
                        print("find again!")
                        break
            else:
                print("error")


    return drop_list


print(len(gdf))
drop=topology_check3(gdf)
print(len(gdf))
print(len(drop))
gdf.drop(axis=0,index=drop,inplace=True)
gdf=gdf.explode()
print(len(gdf))
gdf=gdf[gdf.boundary.is_simple]
gdf=gdf[~(gdf.is_empty)]
gdf.reset_index(drop=True,inplace=True)
print(len(gdf))
print(len(check_overlap(gdf)))