import numpy as np import matplotlib.pyplot as plt from shapely.geometry import Polygon import geopandas as gpd import math def make_mesh(box,w,h): #按照w,h对box进行网格划分 [xmin,ymin,xmax,ymax]=box
list_x=np.arange(xmin,xmax,w) list_y=np.arange(ymin,ymax,h) polygon_list=[] for i in range(len(list_x)): for j in range(len(list_y)): xleft=list_x[i] ydown=list_y[j] if i==len(list_x)-1: xright=xmax else: xright=list_x[i+1] if j==len(list_y)-1: yup=ymax else: yup=list_y[j+1] rectangle=Polygon([(xleft, ydown), (xright, ydown), (xright, yup), (xleft, yup)]) polygon_list.append(rectangle) return gpd.GeoSeries(polygon_list) path="f:\\shapefile\\OUTPUT\\test.shp" gdf=gpd.read_file(path) print(gdf.total_bounds) [xmin,ymin,xmax,ymax]=gdf.total_bounds
w=(xmax-xmin)/10 h=(ymax-ymin)/10 #[xmin,ymin,xmax,ymax]=eval(input("输入box:")) #[w,h]=eval(input("输入网格宽和高:")) rect=make_mesh([xmin,ymin,xmax,ymax],w,h) ax=gdf.plot() rect.boundary.plot(ax=ax,color='yellow') plt.show() rectangles=gpd.GeoDataFrame(geometry=rect) rectangles['RectangleID']=range(len(rectangles)) res=gpd.sjoin(left_df=rectangles, right_df=gdf, how='right', op='intersects') g=res.groupby('RectangleID') g=dict(list(g)) count=0 ax=gdf.plot() color=['yellow','green','blue','black','purple','red'] mean_heights=[] for k,v in g.items(): row=k//10 col=k%10 mean_heights.append((row,col,v['height'].mean())) print("row={},col={},mean height={}".format(row,col,v["height"].mean())) gpd.GeoSeries(rectangles['geometry'][k]).boundary.plot(ax=ax,color=color[(row%6+col)%6]) v.plot(ax=ax,color=color[(row%6+col)%6]) plt.show() print(res.columns.tolist()) print(res.loc[:,['index_left','RectangleID']]) for mean_height in mean_heights: print(mean_height)
这里需要关注的是geopandas.sjoin这个函数的使用，这个函数非常强大，详细使用方法见链接
上述代码实现了:
对地图文件按照矩形进行分割，得到一个网状结构，然后统计每个网格内建筑物的平均高度。
以下是运行结果：
这是对网络的分割，下面是获取每个网格内的建筑物：
这是每个网格内建筑物的平均高度（每个网格以元祖存储）：
