使用python实现多光谱区域统计分析

基于geopandas和rasterio实现区域统计分析，批量获取我们感兴趣的ROI的统计值

数据准备

假设我们有一张多光谱影像和一个shp图，shp是一些取样点或者小区的面信息

这个shp中的属性表信息如下：

这里面其实只需要name信息即可

具体代码

# coding=utf-8
"""
Author: Shuaijie
Blog: https://blog.csdn.net/weixin_45452300
公众号：AgbioIT
date: 2023/12/8 9:29
desc: 根据shp图的名字来区域统计各个波段的值
"""
import geopandas as gpd
import rasterio
import pandas as pd
from rasterio.windows import from_bounds
from rasterio.mask import mask

# 读取Shapefile文件
shp_file = 'path/to/shapefile.shp'
shape_data = gpd.read_file(shp_file)

# 读取多光谱影像
image_file =  'path/to/original_image.tif'
image_data = rasterio.open(image_file)

# 获取影像的波段数量
num_bands = image_data.count
# 创建一个空的DataFrame用于保存结果
result_df = pd.DataFrame(columns=['name'] + [f'Band_{i}' for i in range(1, num_bands+1)])

# 遍历每个区域
for index, region in shape_data.iterrows():
    region_id = region['name']  

    # 提取区域的几何形状
    region_shape = region['geometry']
    if region_shape == None:
        continue
    # 在影像上裁剪出区域的部分
    mask_shape = region_shape.__geo_interface__
    region_image, masked_transform = mask(image_data, [mask_shape], all_touched=False, crop=True)
# 基于两点位置矩形裁剪
#     minx, miny, maxx, maxy = region['geometry'].bounds
#     region_image = image_data.read(window=from_bounds(minx, miny, maxx, maxy, image_data.transform))

    # 统计每个波段的区域像素值
    band_stats = []
    data_dic = {}
    data_dic['name'] = region_id
    for band in range(num_bands):
        band_values = region_image[band, :, :]
        # 将四周没有切进去的值替换为空值
        band_values = np.where(band_values, band_values, np.nan)
        # 计算除去空值后的均值
        data_dic['Band_'+str(band+1)] = np.nanmean(band_values)
        # 计算像元数
#         data_dic['Band_'+str(band+1)] = np.count_nonzero(~np.isnan(band_values))       
    band_stats.append(data_dic.copy())
    # 将统计结果添加到DataFrame中
    result_df = result_df.append(pd.DataFrame(band_stats))

# 保存结果为表格
output_file = 'path/to/your_excel.xlsx'
result_df.to_excel(output_file, index=False)

结果