【GDAL应用】基于rasterstats的矢量数据分区统计栅格值信息

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 ogr


def 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.transform

    zs = 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 = 0
    feature = 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:
                pass
        count += 1
        feature = 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中相应多边形区域相关的统计信息。