### 使用Python实现基于五天滑动窗口法和伽玛分布拟合计算逐日降水数据中每个格点每日的SPI值 为了实现这一目标，主要分为几个部分处理： - **读取并预处理数据**：确保有合适的结构存储多年逐日降水量的数据。 - **应用五天滑动窗口方法**：对于每一天，考虑其前后两天共五个连续日子内的总降水量作为输入给定日期的新特征向量。 - **利用Scipy库进行伽玛分布参数估计**：针对每一个新创建的时间序列片段执行最大似然估计得到形状(shape)、位置(location)及尺度(scale)三个参数[^1]。 - **标准化转换获得SPI分数**：依据上述求得的最佳匹配模型，将原始观测转化为标准正态变量。 下面给出具体代码框架用于指导实践操作： ```python import pandas as pd from scipy.stats import gamma import numpy as np def moving_window_sum(data_series, window_size=5): """Calculate sum over a sliding window.""" return data_series.rolling(window=window_size).sum().shift(-(window_size//2)) def calculate_spi(precipitation_data): """ Calculate SPI values for each day using Gamma distribution. Parameters: precipitation_data (DataFrame): DataFrame containing daily precipitation with columns ['date', 'precip']. Returns: DataFrame: Original dataframe extended by an additional column named 'spi'. """ # Ensure date is datetime type and sort index accordingly precipitation_data['date'] = pd.to_datetime(precipitation_data['date']) precipitation_data.set_index('date', inplace=True) # Apply five-day rolling sum on precipitations precipitation_data['five_day_precip'] = moving_window_sum( precipitation_data['precip'], 5) spi_values = [] dates = [] unique_dates = sorted(set(precipitation_data.index.date)) for current_date in unique_dates: mask = precipitation_data.index.date == current_date if not any(mask): continue subset = precipitation_data.loc[mask] try: shape, loc, scale = gamma.fit(subset['five_day_precip'].dropna(), floc=0) cdf_value = gamma.cdf(subset['five_day_precip'], shape, loc, scale)[0] if cdf_value >= 0.5: spi = norm.ppf(cdf_value) else: spi = -norm.ppf(1-cdf_value) spi_values.append(spi) dates.append(current_date) except Exception as e: print(f"Error processing {current_date}: ", str(e)) pass result_df = pd.DataFrame({'date':dates,'spi':spi_values}) result_df.set_index('date', inplace=True) final_result = precipitation_data.join(result_df, how='left') return final_result.reset_index() if __name__ == '__main__': df = pd.read_csv('daily_rainfall_1980_to_2022.csv') # Replace path/to/yourfile.csv with actual file location results_with_spi = calculate_spi(df) results_with_spi.to_csv('results_with_spi.csv', index=False) ``` 这段脚本首先导入必要的包，并定义辅助函数`moving_window_sum()`来进行移动窗口累加运算；接着核心逻辑封装于`calculate_spi()`内完成整个流程——从准备到最终输出带有SPI指标的结果表单。注意这里假设输入CSV文件至少包含两列：“date”代表记录的具体时间，“precip”表示对应时刻测得的降雨量数值。