把tif的值映射到shp中

问题描述

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

代码

import geopandas as gpd
import rasterio
from rasterstats import zonal_stats
import os

shapefile_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_values

for tif_file in ndvi_dict:
    shapefile[f'NDVI_{tif_file}'] = ndvi_dict[tif_file]

shapefile.to_file('updated_shapefile.shp')

结果示例