geopandas拓扑检查(任意两个几何体不相交)
两个 检查 任意 拓扑 相交
2023-09-14 09:16:18 时间
之前使用rtree来进行拓扑检查,主要是使用GeoDataFrame的sindex来实现的,但是由于某种未知的原因,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)))
相关文章
- [转]java 计算两个时间差得出:时分秒
- Java实现 LeetCode 350 两个数组的交集 II(二)
- (剑指Offer)面试题50:树中两个结点的最低公共祖先
- LeetCode-1662. 检查两个字符串数组是否相等【数组,双指针,字符串】
- 合并两个有序数组
- 83. 一静一动,一张一弛 - 通过具体的两个例子,学习 ABAP 动态断点的使用诀窍
- Python语言学习:利用python获取当前/上级/上上级目录路径(获取路径下的最后叶目录的文件名、合并两个不同路径下图片文件名等目录/路径案例、正确加载图片路径)之详细攻略
- 数据库的两个好帮手:pagehack和pg_xlogdump
- 合并两个有序的链表Python和Go代码
- Linux下的有名管道---使用两个管道实现两个进程之间的通信(手机模式)
- Android两个注意事项.深入了解Intent和IntentFilter(两)
- udp内网穿透 两个内网互联
- LabVIEW中比较两个VI