問題描述

假如目前有一個（多個）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')

結果示例