从多光谱到高光谱:如何用 Python 进行植被指数(NDVI/EVI)时间序列特征提取

NDVI = (ρNIR - ρRed) / (ρNIR + ρRed)

EVI = G * (ρNIR - ρRed) / (ρNIR + C1 * ρRed - C2 * ρBlue + L)

import rasterioimport numpy as npimport osimport globimport matplotlib.pyplot as pltdef read_band_from_stack(stack_path, band_index):    """    从多波段栅格文件中读取指定索引的波段数据。    参数:        stack_path (str): 输入的多波段 TIF 文件路径。        band_index (int): 要读取的波段索引 (1-based)。    返回:        numpy.ndarray: 波段数据数组。    """    with rasterio.open(stack_path) as src:        band_data = src.read(band_index)        # 获取元数据用于后续写入或参考        profile = src.profile        transform = src.transform        crs = src.crs    return band_data, profile, transform, crs

def calculate_ndvi(nir_band, red_band):    """    计算 NDVI。    """    # 避免除以零    with np.errstate(divide='ignore', invalid='ignore'):        ndvi = (nir_band.astype(float) - red_band.astype(float)) / (nir_band.astype(float) + red_band.astype(float))        # 将无效值（如 NaN）设置为 -1 或其他特定填充值        ndvi[np.isnan(ndvi)] = -1.0    return ndvidef calculate_evi(nir_band, red_band, blue_band, G=2.5, C1=6, C2=7.5, L=1):    """    计算 EVI。    """    with np.errstate(divide='ignore', invalid='ignore'):        evi = G * (nir_band.astype(float) - red_band.astype(float)) / \              (nir_band.astype(float) + C1 * red_band.astype(float) - C2 * blue_band.astype(float) + L)        evi[np.isnan(evi)] = -1.0    return evi

def process_annual_timeseries(monthly_stack_paths, output_dir, pixel_coords=None):    """    处理一年的多时相影像，计算每月 NDVI 和 EVI。    参数:        monthly_stack_paths (list): 包含12个月影像路径的列表，需按时间排序。        output_dir (str): 输出目录。        pixel_coords (tuple, optional): 特定像元的 (row, col) 坐标，用于提取物候曲线。    """    os.makedirs(output_dir, exist_ok=True)    monthly_ndvi = []    monthly_evi = []    time_axis = []    for i, path in enumerate(monthly_stack_paths):        month = i + 1        print(f"正在处理第 {month} 月数据: {os.path.basename(path)}")        # 读取三个波段 (假设索引：1-Blue, 2-Red, 3-NIR)        blue, _, _ = read_band_from_stack(path, 1)  # 蓝波段        red, _, _ = read_band_from_stack(path, 2)   # 红波段        nir, profile, transform, _ = read_band_from_stack(path, 3)  # 近红外波段        # 计算指数        ndvi_month = calculate_ndvi(nir, red)        evi_month = calculate_evi(nir, red, blue)        # 保存结果为新的 GeoTIFF        ndvi_out_path = os.path.join(output_dir, f'NDVI_2023_{month:02d}.tif')        evi_out_path = os.path.join(output_dir, f'EVI_2023_{month:02d}.tif')        # 更新 profile 为单波段        profile.update(count=1, dtype='float32')        with rasterio.open(ndvi_out_path, 'w', **profile) as dst:            dst.write(ndvi_month.astype(np.float32), 1)        with rasterio.open(evi_out_path, 'w', **profile) as dst:            dst.write(evi_month.astype(np.float32), 1)        # 收集时间序列数据        monthly_ndvi.append(ndvi_month)        monthly_evi.append(evi_month)        time_axis.append(month)        # 如果指定了特定像元，提取该点的时间序列值        if pixel_coords is not None:            row, col = pixel_coords            # 检查边界            if 0 <= row < ndvi_month.shape[0] and 0 <= col < ndvi_month.shape[1]:                ndvi_val = ndvi_month[row, col]                evi_val = evi_month[row, col]                print(f"  像元 ({row},{col}) - NDVI: {ndvi_val:.3f}, EVI: {evi_val:.3f}")    return monthly_ndvi, monthly_evi, time_axis# --- 主程序入口 ---if __name__ == '__main__':    # 假设数据目录结构    DATA_DIR = r'E:\Data\Sentinel2_Monthly'    OUTPUT_DIR = r'E:\Data\VegetationIndex_TimeSeries'    # 按月份排序获取文件列表    monthly_files = sorted(glob.glob(os.path.join(DATA_DIR, '*_stack.tif')))    if len(monthly_files) != 12:        print(f"警告：预期12个月的数据，实际找到 {len(monthly_files)} 个文件。")    else:        # 假设我们想追踪一个位于影像中心的像元        # 首先获取影像尺寸        with rasterio.open(monthly_files[0]) as src:            height, width = src.shape        center_pixel = (height // 2, width // 2)        print(f"将追踪像元坐标: {center_pixel}")        # 执行处理        ndvi_stack, evi_stack, months = process_annual_timeseries(            monthly_files, OUTPUT_DIR, pixel_coords=center_pixel        )        print("年度时间序列处理完成！")        print(f"NDVI 影像已保存至: {OUTPUT_DIR}")

def plot_phenology_curve(months, ndvi_stack, evi_stack, pixel_coords):    """    绘制指定像元的 NDVI/EVI 物候曲线。    """    row, col = pixel_coords    ndvi_pixel = [ndvi_stack[m][row, col] for m in range(len(months))]    evi_pixel = [evi_stack[m][row, col] for m in range(len(months))]    plt.figure(figsize=(12, 6))    plt.plot(months, ndvi_pixel, 'g-o', linewidth=2, markersize=6, label='NDVI')    plt.plot(months, evi_pixel, 'r-s', linewidth=2, markersize=6, label='EVI')    plt.xlabel('月份', fontsize=12)    plt.ylabel('指数值', fontsize=12)    plt.title(f'像元 ({row}, {col}) 的年度植被指数物候曲线', fontsize=14)    plt.legend(fontsize=11)    plt.grid(True, alpha=0.3)    plt.xticks(months)    plt.ylim(-0.1, 1.0)    plt.tight_layout()    # 保存图表    plt.savefig(os.path.join(OUTPUT_DIR, f'phenology_curve_{row}_{col}.png'), dpi=150)    plt.show()    print(f"物候曲线图已保存。")# 在主程序末尾调用# plot_phenology_curve(months, ndvi_stack, evi_stack, center_pixel)