1 实现效果

矢量数据：

栅格数据：只有一个value值（像素值或DN值），为1，计算统计时nodata作为0值处理。

输出结果：

2 实现功能

基于单波段的栅格数据（一般常为分类数据）和矢量面要素数据，计算矢量数据内栅格数据的统计值（如最大值、平均值、总和、最小值等）。

3 实现代码

# -*- coding: utf-8 -*-
# -*- coding: utf-8 -*-
"""
@time: 2023-10-24 9:04
@author: RSer_gis
@description: 分区统计
计算栅格数据与矢量数据中相应多边形区域相关的统计信息。
"""import time
import rasterio
from rasterstats import zonal_stats
from osgeo import ogrdef create_field(shp_file, stats_list=['majority']):'''功能：为矢量数据创建属性字段:param shp_file:   输入矢量文件shapfile文件:param stats_list: 待添加字段列表:return:'''driver = ogr.GetDriverByName('ESRI Shapefile')layer_source = driver.Open(shp_file, 1)lyr = layer_source.GetLayer()defn = lyr.GetLayerDefn()# 获取图层字段数量featureCount = defn.GetFieldCount()exists_fields = []  # 创建一个空列表，用于存储已存在的字段名称for i in range(featureCount):field = defn.GetFieldDefn(i)  # 获取图层字段的定义（描述字段的信息）field_name = field.GetNameRef()  # 获取字段的名称exists_fields.append(field_name)  # 将当前字段的名称添加到exists_fields列表，以便后续检查字段是否存在# 判断待添加字段是否已存在当前矢量数据属性表（字段）中for ele in stats_list:if ele in exists_fields:print(f"{ele}字段已存在当前矢量数据属性表中")else:cls_name = ogr.FieldDefn(ele, ogr.OFTReal)  # 创建字段描述信息（字段结构），制定名称，数据类型# cls_name.SetWidth(64)lyr.CreateField(cls_name)  # 在图层中创建新字段driver = None  # 在创建字段后，将驱动程序设为“None”以释放资源def set_fieldvalue(raster_file, shp_file, stats_list=['majority']):''':功能   基于矢量要素，为新添加的字段赋值统计信息:param raster_file:  输入栅格数据tif文件:param shp_file:  输入矢量数据shapfile文件:param stats_list: 为添加的字段赋值:return:'''ras_driver = rasterio.open(raster_file, driver="GTiff", nodata=0)array = ras_driver.read(1)  # 如果是多个波段，只读取第1个波段；波段索引从1开始affine = ras_driver.transformzs = zonal_stats(shp_file, array, affine=affine, all_touched=True,stats=stats_list,nodata=0)  # 每个要素的统计信息以字典（Dictionary）的形式存储在列表中 zs = [{'min': 1.0, 'max': 5.0, 'mean': 2.5},{'min': 0.0, 'max': 3.0, 'mean': 1.5},...]driver = ogr.GetDriverByName('ESRI Shapefile')layer_source = driver.Open(shp_file, 1)lyr = layer_source.GetLayer()defn = lyr.GetLayerDefn()featureCount = defn.GetFieldCount()count = 0feature = lyr.GetNextFeature()  # 获取第一个要素while feature is not None:for i in range(featureCount):field = defn.GetFieldDefn(i)field_name = field.GetNameRef()# 判断图层要素中的字段是否存在于stats_list，如果存在，将相应的统计信息设置到要素的相应字段中if field_name in stats_list:feature.SetField(field_name, zs[count][field_name])lyr.SetFeature(feature)else:passcount += 1feature = lyr.GetNextFeature()  # 获取下一个要素# 清理资源if layer_source is not None:layer_source = None  # 关闭数据源并释放相关资源driver = None  # 释放驱动程序引用ras_driver.close()  # 关闭数据集if __name__ == "__main__":start = time.time()rasterPath = r"D:\Desktop\test\data\潍坊市1702361169438\WF2024.tif"shpPath = r"D:\Desktop\test\data\gd_AL.shp"# stats_list = ['min', 'max', 'mean', 'median', 'majority'] # 统计信息字段名称stats_list = ['majority', 'max']print("创建字段...".center(25, "-"))create_field(shpPath, stats_list)print("字段赋值...".center(25, "-"))set_fieldvalue(rasterPath, shpPath, stats_list)print("分区统计完成！".center(25, "-"))end = time.time()print((end - start) / 60.0)

说明：
zonal_stats 函数的参数解释如下：

• shp_path: 矢量数据的路径（如shapefile）。
• array: 栅格数据的NumPy数组。
• affine: 仿射变换矩阵，用于将栅格坐标映射到地理坐标。这通常是一个包含六个元素的元组或列表，如 (xoff, a, b, yoff, d, e)，其中 xoff 和 yoff 是左上角的x和y坐标，a 是x方向的像素分辨率（通常为负），d是y方向的像素分辨率（通常为正），b 和 e 通常是0（表示没有旋转或倾斜）。
• all_touched:
  如果为True，则即使多边形只与栅格的一个像素相交（而不仅仅是完全覆盖它），也会返回该像素的统计信息。如果为False，则只返回完全在多边形内的像素的统计信息。
• categorical: 如果为True，并且array是一个分类栅格（即每个值代表一个类别而不是一个数值），则返回每个类别的计数，而不是数值统计（如平均值、标准差等）。
• nodata: 用于标识array中无数据或缺失值的值。

zonal_stats 函数返回一个列表，其中每个元素都是一个字典，表示与shp_path中相应多边形区域相关的统计信息。