企业官网建设流程全解析

简介：这份资源面向气象、水文、生态及农业领域的科研人员与数据分析学习者，提供多时间尺度、多分辨率、多PET计算方式的日/周/月干旱指标SPEI计算代码及配套测试文件。SPEI通过降雨与潜在蒸散发的水平衡累积反映不同时长干旱，3个月尺度对应农业与土壤干旱，6个月尺度对应水文干旱，而日分辨率SPEI能捕捉仅持续数周的短期干旱事件，弥补月尺度指标的不足。资源包共22个文件，包含10个txt说明与数据文件、4个csv测试数据、3个py计算与绘图脚本、3张jpg结果示意图及2个pyc缓存文件，压缩包约1.07MB，结构清晰，覆盖输入数据、计算流程与输出结果。目前已有1533人学习下载。读者可据此复现日尺度SPEI构建方法，掌握多PET方式下的指标计算，并基于游程理论提取干旱事件、分析区域干湿趋势，适合具备一定Python基础、希望深入干旱监测研究的学习者参考使用。

1. 干旱指数计算工具：从日尺度到月尺度的 SPEI 全流程拆解

做气象、水文或者农业干旱监测的开发者，大概率都绕不开 SPEI（标准化降水蒸散指数）这个指标。它比 SPI 多考虑了蒸散项，能同时反映降水和温度对干旱的贡献，所以在气候变暖背景下的适用性更好。但真正动手算的时候，问题就来了：不同业务场景要的时间尺度不一样，有的要日尺度预警，有的要周尺度评估，有的要月尺度做长期趋势；空间分辨率也五花八门，站点数据、格点数据、区域汇总数据混在一起；PET（潜在蒸散）的计算方式更是有好几种，Thornthwaite、Hargreaves、Penman-Monteith 各有各的适用条件。这套代码就是冲着这个痛点来的——把多时间尺度、多分辨率、多 PET 计算方式的 SPEI 计算流程整合成一套可配置、可测试的工程代码，附带测试文件，适合做干旱监测系统二次开发或者科研复现的开发者直接上手。

2. 项目架构与目录结构：模块怎么分、数据怎么流

2.1 整体模块划分

这套代码的核心思路是「配置驱动 + 模块解耦」。整个项目围绕三个维度展开：时间尺度（日/周/月）、空间分辨率（站点/格点/区域）、PET 计算方式（Thornthwaite/Hargreaves/Penman-Monteith）。如果按传统写法，三个维度交叉会产生大量重复代码，改一个参数要动好几个文件。所以项目在架构上做了分层：

• 数据层：负责读取不同分辨率的输入数据，统一转成内部标准格式（时间序列 + 空间索引）。
• 计算层：核心算法模块，包括 PET 计算、水分亏缺量计算、概率分布拟合、SPEI 标准化。
• 配置层：用配置文件或参数字典控制时间尺度、分辨率、PET 方式的选择。
• 测试层：提供模拟数据和预期结果，验证各组合下的计算正确性。

这种分层的好处是，新增一种 PET 计算方式只需要在计算层加一个函数，配置层加一个选项，不用动数据层和测试层的核心逻辑。

2.2 目录结构说明

典型的目录组织如下（不同版本可能略有差异，但核心结构一致）：

spei_calculator/ ├── config/ │ └── settings.yaml # 全局配置：时间尺度、PET方式、分辨率类型 ├── data/ │ ├── sample_station.csv # 站点示例数据 │ ├── sample_grid.nc # 格点示例数据 │ └── sample_region.csv # 区域汇总示例数据 ├── src/ │ ├── data_loader.py # 数据读取与格式统一 │ ├── pet_models.py # 三种PET计算方式实现 │ ├── spei_core.py # SPEI核心计算流程 │ ├── utils.py # 时间尺度转换、缺失值处理等工具 │ └── main.py # 主入口，串联全流程 ├── tests/ │ ├── test_pet.py # PET计算单元测试 │ ├── test_spei.py # SPEI计算单元测试 │ └── test_data/ # 测试用模拟数据 ├── requirements.txt └── README.md

这个结构里，src/是核心，tests/是验证，config/是控制入口。实际二次开发时，大部分改动集中在pet_models.py和spei_core.py，数据层和测试层相对稳定。

2.3 数据流与调用链路

一次完整的 SPEI 计算，数据流是这样的：

1. main.py读取settings.yaml，确定本次计算的时间尺度（比如月尺度）、分辨率类型（比如站点）、PET 方式（比如 Thornthwaite）。
2. data_loader.py根据分辨率类型加载对应数据，统一成「时间 × 空间」的二维结构，时间索引对齐到目标尺度。
3. pet_models.py根据配置调用对应的 PET 计算函数，输出与降水数据同维度的 PET 序列。
4. spei_core.py计算水分亏缺量 D = P - PET，然后对 D 序列做概率分布拟合（常用 Log-logistic 分布），再标准化得到 SPEI。
5. 结果写回文件，同时tests/下的测试用例验证关键步骤的数值正确性。

这个链路里，时间尺度转换和 PET 计算是两个最容易出问题的环节，后面会单独展开。

3. 核心模块拆解：PET 计算与 SPEI 标准化怎么实现

3.1 三种 PET 计算方式的代码实现

PET 是 SPEI 计算的前置步骤，选哪种方式直接影响最终结果。项目里把三种常用方式都实现了，通过配置切换。

Thornthwaite 方式：只依赖月平均气温和纬度，适合气温数据完整但其他气象要素缺失的场景。

def pet_thornthwaite(temp_monthly, latitude, heat_index=None): """ Thornthwaite 方法计算月 PET temp_monthly: 月平均气温序列（摄氏度） latitude: 站点纬度（度） 返回: 月 PET 序列（mm/month） """ # 计算热指数 I，如果未提供则从气温序列估算 if heat_index is None: # 只对气温大于0的月份累加 valid = temp_monthly[temp_monthly > 0] heat_index = np.sum((valid / 5.0) ** 1.514) # 经验系数 a a = (6.75e-7 * heat_index**3 - 7.71e-5 * heat_index**2 + 1.792e-2 * heat_index + 0.49239) # 基础 PET 计算 pet_base = 16 * ((10 * temp_monthly / heat_index) ** a) pet_base[temp_monthly <= 0] = 0 # 气温非正时 PET 为 0 # 根据纬度做日长修正 lat_rad = np.radians(latitude) day_length = 24 / np.pi * np.arccos( np.clip(-np.tan(lat_rad) * np.tan(0.409 * np.sin(2*np.pi*(np.arange(12)+1)/12 - 1.39)), -1, 1) ) correction = day_length / 12.0 return pet_base * correction

这段代码的关键参数是heat_index（热指数），它反映了一个地区全年的热量积累。如果输入数据跨年，热指数应该用多年平均来算，而不是逐年单独算，否则年际之间的 PET 会有系统性偏差。纬度修正部分用的是 FAO 推荐的日长公式，np.clip防止极地地区出现无效值。

Hargreaves 方式：需要气温和纬度，但增加了气温日较差项，在干旱区表现更好。

def pet_hargreaves(temp_mean, temp_max, temp_min, latitude, extraterrestrial_rad): """ Hargreaves 方法计算 PET temp_mean/temp_max/temp_min: 日均/最高/最低气温 extraterrestrial_rad: 天文辐射（MJ/m2/day） """ # 气温日较差 temp_range = temp_max - temp_min # 经验系数 0.0023，温度项 0.5 为偏移 pet = 0.0023 * extraterrestrial_rad * (temp_mean + 17.8) * np.sqrt(temp_range) return pet

Hargreaves 的核心优势是只需要气温数据，extraterrestrial_rad可以根据纬度和日序计算得到。注意temp_range如果为负（数据异常），需要做截断处理，否则开根号会出问题。

Penman-Monteith 方式：最完整但也最复杂，需要气温、湿度、风速、辐射四类数据。

def pet_penman_monteith(temp, humidity, wind_speed, radiation, pressure=101.3, albedo=0.23): """ FAO-56 Penman-Monteith 计算日 PET 所有单位需统一：temp(°C), humidity(%), wind(m/s), radiation(MJ/m2/day) """ # 饱和水汽压曲线斜率 delta = 4098 * (0.6108 * np.exp(17.27 * temp / (temp + 237.3))) / (temp + 237.3)**2 # 实际水汽压 ea = 0.6108 * np.exp(17.27 * temp / (temp + 237.3)) * humidity / 100 # 净辐射 rn = (1 - albedo) * radiation - 2.45e-9 * (temp + 273.16)**4 * (0.34 - 0.14 * np.sqrt(ea)) # 土壤热通量（日尺度近似为0） g = 0 # 干湿表常数 gamma = 0.665e-3 * pressure # 最终公式 pet = (0.408 * delta * (rn - g) + gamma * 900 / (temp + 273) * wind_speed * (0.6108 * np.exp(17.27*temp/(temp+237.3)) - ea)) / (delta + gamma * (1 + 0.34 * wind_speed)) return pet

Penman-Monteith 的参数最多，也最容易因为单位不统一而出错。pressure默认 101.3 kPa 是海平面标准大气压，高海拔地区需要按海拔修正。albedo默认 0.23 是参考作物表面反射率，一般不用改。

3.2 SPEI 标准化流程与分布拟合

PET 算完之后，SPEI 的核心步骤是：水分亏缺量 D = P - PET，然后对 D 序列拟合概率分布，再标准化。

def calculate_spei(precip, pet, scale=1, dist='loglogistic'): """ 计算指定时间尺度的 SPEI precip/pet: 同维度的降水与PET序列 scale: 时间尺度（1=月，3=季，12=年） dist: 概率分布类型 """ # 水分亏缺量 d = precip - pet # 按尺度做滑动累加 if scale > 1: d = np.convolve(d, np.ones(scale), mode='valid') # 拟合分布 if dist == 'loglogistic': # 三参数 Log-logistic 分布拟合 params = fit_loglogistic(d) cdf = cdf_loglogistic(d, *params) elif dist == 'gamma': params = fit_gamma(d) cdf = cdf_gamma(d, *params) # 标准化：标准正态分布的分位数 spei = norm.ppf(cdf) return spei

这里有几个关键点。第一，scale参数控制时间尺度，scale=1是月尺度，scale=3是季尺度，scale=12是年尺度。滑动累加用np.convolve实现，注意mode='valid'会让输出长度减少scale-1，时间索引要对应调整。第二，分布拟合用的是 Log-logistic，这是 SPEI 的推荐分布，但样本量小于 30 时拟合不稳定，需要做参数约束。第三，norm.ppf在 cdf 接近 0 或 1 时会返回无穷大，实际代码里要做截断，比如限制在 0.001 到 0.999 之间。

3.3 多时间尺度与多分辨率的配置切换

项目用配置文件控制不同组合，避免硬编码。settings.yaml的典型内容：

time_scale: monthly # daily / weekly / monthly resolution: station # station / grid / region pet_method: thornthwaite # thornthwaite / hargreaves / penman_monteith distribution: loglogistic output_dir: ./output

main.py里根据配置分发：

def run_spei(config): # 加载数据 data = load_data(config['resolution'], config['data_path']) # 时间尺度转换 data = resample_time(data, config['time_scale']) # PET 计算 pet_func = PET_METHODS[config['pet_method']] pet = pet_func(data['temp'], data['latitude'], **data.get('extra', {})) # SPEI 计算 spei = calculate_spei(data['precip'], pet, scale=config.get('scale', 1)) # 输出 save_result(spei, config['output_dir'])

这种配置驱动的写法，新增一种 PET 方式只需要在PET_METHODS字典里注册一个函数，新增一种分辨率只需要在load_data里加一个分支。实际二次开发时，建议把配置项做成命令行参数或者环境变量，方便批量跑不同组合。

4. 运行与测试：从零跑通一套 SPEI 计算

4.1 环境准备与依赖安装

项目依赖主要是科学计算栈：numpy、scipy、pandas、pyyaml，如果处理格点数据还需要xarray和netCDF4。requirements.txt里一般会固定版本，避免因为 scipy 版本差异导致分布拟合结果不一致。

# 创建虚拟环境 python -m venv venv source venv/bin/activate # Windows 用 venv\Scripts\activate # 安装依赖 pip install -r requirements.txt

如果只是跑测试，pytest也要装上。格点数据那块，netCDF4在 Windows 上偶尔会有编译问题，常见做法是用conda install netcdf4代替 pip。

4.2 用测试文件验证计算正确性

项目自带的tests/目录是验证环境是否搭好的最快方式。测试用例通常包含模拟的站点数据、已知的 PET 预期值、以及 SPEI 的参考结果。

# 跑全部测试 pytest tests/ -v # 只跑 PET 测试 pytest tests/test_pet.py -v # 跑特定测试函数 pytest tests/test_spei.py::test_monthly_spei -v

测试文件里一般会构造一个理想化的输入，比如 12 个月的恒定气温和降水，然后断言 PET 和 SPEI 的输出在合理范围内。如果测试全绿，说明环境没问题，可以开始跑真实数据。

4.3 跑通一次完整的月尺度 SPEI 计算

以站点数据为例，完整流程如下：

# 修改配置 vim config/settings.yaml # 设置 time_scale: monthly, resolution: station, pet_method: thornthwaite # 运行主程序 python src/main.py --config config/settings.yaml --data data/sample_station.csv --output ./output

运行成功后，output/下会生成 SPEI 结果文件，通常是 CSV 格式，包含时间列和 SPEI 值列。如果要做日尺度或周尺度，改time_scale配置即可，但要注意日尺度下 PET 计算需要日值气象数据，Thornthwaite 方法在日尺度下精度有限，建议换 Hargreaves 或 Penman-Monteith。

4.4 结果解读与常见数值范围

SPEI 的输出是标准化值，理论上服从标准正态分布。一般解读标准：

如果算出来的 SPEI 大量超出 ±3，大概率是分布拟合出了问题，检查样本量是否足够、是否有异常值。

5. 避坑与排查：SPEI 计算中最容易翻车的几个点

5.1 时间尺度转换后时间索引错位

现象：月尺度转季尺度后，结果的时间标签对不上，比如 3 月的数据被标成了 1 月。

原因：滑动累加用np.convolve时，mode='valid'会丢弃前scale-1个值，但时间索引没有同步调整。

解决：累加后手动截取时间索引，time_index[scale-1:]，或者用 pandas 的rolling方法，它会自动对齐索引。

5.2 PET 计算中单位不统一导致结果偏差

现象：Penman-Monteith 算出来的 PET 比 Thornthwaite 大一个数量级。

原因：辐射单位用了 MJ/m²/day 但公式里按 W/m² 处理，或者风速用了 km/h 但公式要求 m/s。

解决：在函数入口处加单位检查和转换，所有输入统一到 FAO-56 标准单位。常见做法是在data_loader.py里做一次全局单位归一化。

5.3 分布拟合样本量不足导致 SPEI 异常

现象：短序列（少于 30 个时间步）算出来的 SPEI 出现 ±5 以上的极端值。

原因：Log-logistic 分布的三参数拟合在小样本下不稳定，参数估计方差大。

解决：样本量不足时改用两参数的 Gamma 分布，或者对拟合参数加约束（比如形状参数限制在合理范围内）。项目测试文件里一般会包含小样本的边界测试。

5.4 格点数据分辨率不一致导致广播错误

现象：降水是 0.5° 格点，温度是 0.25° 格点，直接相减报维度不匹配。

原因：不同来源的格点数据空间分辨率不同，没有做重采样对齐。

解决：在data_loader.py里加一步空间重采样，统一到较粗的分辨率，或者用xarray的interp_like方法对齐。注意重采样会引入误差，尽量用原始分辨率一致的数据。

5.5 缺失值处理不当导致整段结果失效

现象：某个月降水数据缺失，导致该月及之后所有尺度的 SPEI 都是 NaN。

原因：滑动累加时缺失值会传播，np.convolve遇到 NaN 会污染整个窗口。

解决：在累加前做缺失值插补，或者用np.nansum配合掩码数组。如果缺失比例超过 10%，建议直接标记该时段为无效，不要强行插补。

6. 二次开发与改造：从能跑到好用

这套代码跑通之后，二次开发的空间主要在三个方向：新增 PET 计算方式、适配自定义数据格式、优化计算性能。

新增 PET 方式是最常见的需求。比如想加一个 Priestley-Taylor 方法，只需要在pet_models.py里写一个函数，然后在PET_METHODS字典里注册，配置里加一个选项就行。函数签名保持和其他 PET 函数一致，输入输出维度对齐，测试文件里补一个单元测试验证数值合理性。

适配自定义数据格式，重点改data_loader.py。如果输入是 Excel 或者自定义二进制格式，在load_data里加一个分支，输出统一成「时间 × 空间」的 DataFrame 或 ndarray。这里有个血泪经验：时间索引一定要用 pandas 的 DatetimeIndex，不要用字符串，否则后续重采样和对齐会非常痛苦。

性能优化方面，如果处理大区域格点数据，逐格点循环会非常慢。常见做法是用xarray的apply_ufunc做向量化，或者用numba加速核心循环。我一般会先用小区域测试向量化结果的正确性，确认无误后再上全量数据。

验证改造是否成功，最直接的方法是跑tests/下的全部用例，然后拿一小段真实数据对比改造前后的输出。如果差异在合理范围内（比如 PET 差异小于 5%），基本可以认为改造没问题。

从那以后我每次改核心计算模块，都强制走一遍「单元测试 + 小样本对比 + 全量跑通」三步，再急也不跳过。希望帮到你。

本文还有配套的精品资源，点击获取