问题描述

假如目前有一个（多个）tif文件和一个shp文件，想要把tif中每个像素的值集成到shp文件的新字段中。如果栅格和像素是一一对应的，问题将会变得非常简单：直接把每个像素的值映射到每个栅格中的新字段里即可。如果栅格和像素不是一一对应的，这里我们需要用到zonal_stats，把所有涉及到的像素值求平均值再映射到栅格里去。接下来是一个简单的示例，假设我们要把不同日期的NDVI值映射到小栅格里面去。

代码

import geopandas as gpd
import rasterio
from rasterstats import zonal_stats
import osshapefile_path = 'your.shp'tif_directory = 'your_path'shapefile = gpd.read_file(shapefile_path)# 初始化一个字典来存储每个网格的 NDVI 值
ndvi_dict = {}for tif_file in os.listdir(tif_directory):if tif_file.endswith('.tif'):tif_path = os.path.join(tif_directory, tif_file)with rasterio.open(tif_path) as src:nodata = src.nodata  # 从 TIFF 文件中获取 nodata 值# 通过更精确的对齐计算区域统计数据stats = zonal_stats(shapefile, tif_path, stats='mean', nodata=nodata, all_touched=True)# 提取平均 NDVI 值mean_ndvi_values = [stat['mean'] if stat['mean'] is not None else 0 for stat in stats]ndvi_dict[tif_file] = mean_ndvi_valuesfor tif_file in ndvi_dict:shapefile[f'NDVI_{tif_file}'] = ndvi_dict[tif_file]shapefile.to_file('updated_shapefile.shp')

结果示例