【技术教程】使用Python进行区域统计

1. 公众号回复0615即可获得下载地址。

2. 为了支持我得学城，欢迎转发，感谢！‍

全文约3500字，预计阅读时间9分钟。

conda activate pygeo38

from pathlib import Path # 处理文件路径
import rasterio as rio
from rasterio.plot import show
from shapely.geometry import box
from rasterstats import zonal_stats

ax = gdf.plot(facecolor='none', edgecolor='red')
with rio.open(dem_file) as f:
    show(f, cmap='terrain', ax=ax)

gdf.crs 

# 确定是投影坐标系

# 获取做500m缓冲区后的外轮廓的外包矩形的四角坐标
cell_width = 5000  # 网格边长
minx, miny, maxx, maxy = gdf.buffer(cell_width  2).bounds.values[0]
# 获取每个小方格的左下角坐标
xmin_lst = np.arange(minx // cell_width, maxx // cell_width, 1) * cell_width
ymin_lst = np.arange(miny // cell_width, maxy // cell_width, 1) * cell_width
# 生成矩阵
x1, y1 = np.meshgrid(xmin_lst, ymin_lst)
# 获取右上角坐标
x2 = x1 + cell_width
y2 = y1 + cell_width
x1, y1, x2, y2 = [i.flatten() for iin [x1, y1, x2, y2]]
# 生成网格
boxes = [box(*i) for i in zip(x1, y1, x2, y2)]
gdf_fishnet = gpd.GeoDataFrame(geometry=boxes, crs=gdf.crs).reset_index()
# 存为shapefile
gdf_fishnet.to_file(data_path / 'fishnet.shp')
gdf_fishnet.shape

(8554, 2)

f, ax = plt.subplots(1,1, figsize=(10,10))
gdf.plot(facecolor='none', edgecolor='red', ax=ax)
with rio.open(dem_file) as f:
    show(f, cmap='terrain', ax=ax)

gdf_fishnet.plot(facecolor='none', edgecolor='blue', ax=ax,lw=0.3)

stat_result = zonal_stats(
    (data_path /'fishnet.shp').__str__(), 

# 无法识别pathlib的类型 转为字符串 也可以直接使用相对路径和绝对路径
    (data_path /'dem.tif').__str__()
    stats="count min mean max"
)

gdf_result = gdf_fishnet.join(pd.DataFrame(stat_result))
gdf_result.plot(column='mean', figsize=(10,10))

1. 科研论文中的技术方法，比如论文中的关键技术；

2. 机器学习与深度学习算法原理，比如像之前介绍PCA与K-means原理的文章；

3. 机器学习与深度学习在科研与业务中的应用，这是后面关注的重点，会有很多文章，敬请期待；

4. 学术论文写作方法与技巧，后面会请学术大咖做客我得学城，介绍论文写作方法与技巧；

5. 业务咨询类的量化分析技术，比如国土空间规划中的相关量化方法与技术。