### 计算 SPEI 空间分布的 Python 代码 为了计算标准化降水蒸发指数 (SPEI)，可以使用 `spei` 库来处理时间序列数据，并利用 `rasterio` 和 `numpy` 来读取和写入 GeoTIFF 文件。以下是完整的 Python 实现： #### 安装所需库 首先安装所需的 Python 包： ```bash pip install rasterio numpy pandas scipy statsmodels spei ``` #### 导入库并定义函数 接下来，导入必要的库并编写主要功能函数。 ```python import os import numpy as np import pandas as pd from datetime import datetime, timedelta import rasterio from rasterio.enums import Resampling from spei import spei def read_tiff(file_path): """Read a single band TIFF file into an array.""" with rasterio.open(file_path) as src: data = src.read(1).astype(float) profile = src.profile.copy() return data, profile def write_tiff(output_file, image_data, meta_profile): """Write the output to a new TIFF file using provided metadata.""" with rasterio.open( output_file, 'w', driver='GTiff', height=image_data.shape[0], width=image_data.shape[1], count=1, dtype=str(image_data.dtype), crs=meta_profile['crs'], transform=meta_profile['transform'] ) as dst: dst.write(image_data, indexes=1) def compute_spei(precipitation_array, pet_array, scale=3): """ Compute SPEI based on precipitation and PET arrays. Parameters: - precipitation_array: Array of monthly precipitation values over time at each pixel location. - pet_array: Array of corresponding potential evapotranspiration values. - scale: Time window length used for calculating moving averages; default is 3 months. Returns: An array containing computed SPEI indices per month. """ # Calculate water balance deficit/surplus series by subtracting PET from PPT wb_series = precipitation_array - pet_array # Convert this difference matrix back into DataFrame format suitable for SPEI calculation df_wb = pd.DataFrame(wb_series.T) # Apply SPEI function row-wise across all pixels within dataframe result_df = df_wb.apply(lambda col: spei(col.values, scales=[scale]).values[-1], axis=0) return result_df.to_numpy().reshape((precipitation_array.shape[1], precipitation_array.shape[2])) # Example usage: input_precipitation_tif = "path/to/your/monthly_precipitation.tif" input_pet_tif = "path/to/your/monthly_potential_evapo_transpiration.tif" output_dir = "./results/" os.makedirs(output_dir, exist_ok=True) # Read input rasters ppt_arr, ppt_meta = read_tiff(input_precipitation_tif)[^2] pet_arr = read_tiff(input_pet_tif)[0] # Ensure both inputs have same dimensions & resolution before proceeding further... assert ppt_arr.shape == pet_arr.shape, f"Mismatched shapes between {input_precipitation_tif} ({ppt_arr.shape}) vs. {input_pet_tif} ({pet_arr.shape}). Please resample one or more datasets." # Perform computation spei_result = compute_spei(ppt_arr, pet_arr) # Save results as geotiffs write_tiff(os.path.join(output_dir, "spei_output.tif"), spei_result, ppt_meta) print(f"Spatial SPEI successfully calculated! Results saved under '{output_dir}'") ``` 此脚本假设输入的时间序列为逐月累积的数据集。如果实际数据不是按月存储，则需调整日期范围以匹配具体需求[^3]。