文章目录
- 1 实现效果
- 2 实现功能
- 3 实现代码
1 实现效果
矢量数据:
栅格数据:只有一个value值(像素值或DN值),为1,计算统计时nodata作为0值处理。
输出结果:
2 实现功能
基于单波段的栅格数据(一般常为分类数据)和矢量面要素数据,计算矢量数据内栅格数据的统计值(如最大值、平均值、总和、最小值等)。
3 实现代码
# -*- 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 gdal, 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)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 # 释放驱动程序引用if __name__ == "__main__":start = time.time()rasterPath = r"D:\Desktop\test\data\test.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中相应多边形区域相关的统计信息。