Python 分区统计 rasterstats zonal s

需求：给定一个栅格数据和一个矢量面数据，计算面数据内栅格数据的统计值（如最大值、平均值、总和、最小值等）。在QGIS和ArcGIS中，有现成的工具可以直接使用，如分区统计（zonal statistic）工具箱

解决方法：如果需要在Python中批量操作，可以采用rasterstats包，源代码参见 https://github.com/perrygeo/python-rasterstats ，使用手册参见https://pythonhosted.org/rasterstats/manual.html#design-goals

安装如下，

conda install rasterio
pip install rasterstats

使用代码，

import pandas as pd
from rasterstats import zonal_stats

def zonal_statistic():
    # 矢量的边界文件 polygon
    shp_poly = r"E:\Data\poly.shp"
    input_path = "E:\\Data\\"
    file = raster.tif

    df_out = pd.DataFrame()
    
    stats = zonal_stats(shp_poly, input_path + file, stats=['min', 'max', 'count', 'sum', 'std'])
    df_out = pd.DataFrame.from_dict(data=stats)
   
    df_out.to_csv(r"E:\Data\prov_stat_34.csv", header=True, index_label='fid', encoding='gbk')

更多使用方法参见：

http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#clip-a-geotiff-with-shapefile

以及一个课程学习资料：
https://www.gis.usu.edu/~chrisg/python/2009/