Geopandas对地图文件按区域进行分割并统计

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这个函数的使用，这个函数非常强大，详细使用方法见链接
上述代码实现了:
对地图文件按照矩形进行分割，得到一个网状结构，然后统计每个网格内建筑物的平均高度。
以下是运行结果：

这是对网络的分割，下面是获取每个网格内的建筑物：

这是每个网格内建筑物的平均高度（每个网格以元祖存储）：