简介ml_drought 是一套面向气候科学机器学习实践的端到端管道资源适合具备Python基础、关注干旱预测与气候数据分析的学习者和研究者。它围绕统一训练/测试数据格式、方法互比与评估展开将数据准备、模型实验等任务拆分到src目录的不同类中并借助Copernicus 2019相关场景与Jupyter Notebook工作流帮助读者理解如何把机器学习方法落到干旱预测任务。压缩包约49.31MB文件总数显示为0、类型明细暂无统计从说明可知核心为Python管道代码与environment.yml环境配置实际文件清单以压缩包为准可通过创建esowc-drought环境后运行。资源已有143人学习适合用来梳理气候机器学习项目的目录组织、实验流程和评估思路也可作为课程设计或研究复现的参考起点。1. 干旱预测项目最难的环节先把「干旱」翻译成模型能吃的标签拿到 ml_drought 这套机器学习干旱预测 Notebook 时我第一反应不是去看它用了什么模型而是翻到构造标签那一段——因为干旱预测这件事模型从来不是瓶颈把「干旱」这个模糊的气候概念翻译成一个可计算、可对齐时间、可回测的数值标签才是真正卡人的地方。这套项目把 ERA5 再分析气象网格数据、干旱指数计算、特征工程和回归/分类建模串成了一条完整的 Jupyter Notebook 流水线适合两类人一类是想把气候数据分析和机器学习接起来的同学另一类是想复现一个「从原始网格到事件预警」完整链路的开发者。它不追求 SOTA 精度价值在于把每个环节都摊开给你看方便你换数据、换区域、换模型。2. 项目定位与技术栈一条跑在 Jupyter Notebook 里的端到端流水线这套资源的形态是 Notebook 集合而非工程化的 Python 包这一点直接决定了它的使用方式和二次开发的起点。它做的事情可以概括成一句话给定一片区域的月度气象网格算出历史干旱指数再用过去若干月的降水、温度、土壤湿度、蒸散等变量去预测未来一到三个月的干旱状态。整条链路是「读数据 → 算标签 → 造特征 → 训练 → 评估 → 出图」每一步都留了可视化单元格跑一遍就能看到中间结果长什么样。为什么值得拆因为大部分气候类机器学习教程要么只教你调 sklearn要么只教你画等值线图真正把「网格数据怎么变成样本表」讲清楚的很少。而这套 Notebook 的中间层——也就是特征表那一层——恰恰是能直接抄走复用到其他气象任务上的部分。2.1 仓库目录与 Notebook 分工这类项目的仓库通常按流水线顺序拆成几个 Notebook命名上会带编号方便按序执行。常见的组织方式是这样ml_drought/ ├── data/ │ ├── raw/ # 原始网格数据NetCDF按月或按年存放 │ ├── interim/ # 裁剪、重采样后的中间结果 │ └── processed/ # 对齐好的样本表parquet / csv ├── notebooks/ │ ├── 01_load_and_clip.ipynb # 读取网格、裁剪研究区、单位换算 │ ├── 02_drought_index.ipynb # 计算 SPI / SPEI 干旱指数 │ ├── 03_feature_build.ipynb # 滞后项、滚动统计、空间邻域特征 │ ├── 04_train_baseline.ipynb # 线性/树模型基线训练 │ ├── 05_evaluate.ipynb # 指标计算、事件命中率、绘图 │ └── 06_predict_future.ipynb # 用最新数据做推理出图 ├── src/ │ ├── indices.py # 干旱指数函数Notebook 里 import 复用 │ ├── features.py # 特征构造函数 │ └── utils.py # 路径、时间对齐、绘图辅助 ├── requirements.txt └── README.md几个要点值得单独说。data/raw和data/processed分开是必须的因为网格数据动辄几个 G一旦在下游 Notebook 里反复重算裁剪一个单元格能跑十分钟。src/目录的存在说明作者没有把所有逻辑都塞在 Notebook 单元格里——这一点很重要凡是需要被两个以上 Notebook 调用的函数比如 SPI 计算都应该下沉到src/否则改一个公式要改三个地方血泪经验。另外要注意 Notebook 的编号顺序代表执行依赖。02_drought_index.ipynb依赖01产出的interim文件如果你跳过01直接跑02报的错大概率是文件找不到而不是公式错别急着怀疑代码。2.2 技术栈选型为什么是 xarray scikit-learn 而不是深度学习技术栈基本是气候数据分析的标准组合库用途关键点xarray读写 NetCDF 网格、按经纬度/时间切片维度命名time, latitude, longitude要统一pandas样本表、时间索引、滚动窗口时间索引必须是DatetimeIndex且按月频率numpy数值计算、滞后拼接注意 NaN 传播scipy.statsgamma 拟合、正态分位数SPI 计算的核心scikit-learn切分、Pipeline、指标、线性/树模型时间序列不能用随机切分xgboost / lightgbm非线性基线网格小样本时容易过拟合matplotlib / cartopy空间分布图投影和纬度顺序容易翻车为什么不用 LSTM、Transformer我的判断是样本量撑不住。一片中等大小的研究区月度数据三十年也就三百六十个时间步再乘以格点数看似很多但格点之间高度相关有效样本远小于表面数字。这种情况下树模型加手工滞后特征通常比深度模型稳而且可解释性更好——你能直接看出「前三个月降水」这个特征贡献最大这在干旱预警场景里比多零点零几的 R² 值钱得多。另一个常被忽略的选型点是标签用回归还是分类。回归预测连续的 SPI 值信息量大但评估麻烦分类预测「是否干旱事件」业务上好用但阈值一改标签全变。常见做法是两条都做主任务回归、副任务分类用回归的输出再卡阈值做事件预警这样阈值可以随时调而不用重训模型。3. 数据管道拆解从 ERA5 网格到可训练特征表这一章是整个项目里最值得细看的部分。建模那一步其实几行 sklearn 就完事了真正决定结果好坏的是数据管道——标签算错、特征泄漏、时间没对齐模型再花哨也是白搭。3.1 网格读取、区域裁剪与单位换算ERA5 这类再分析数据以 NetCDF 存储变量名是缩写tp 表示总降水、t2m 表示两米气温、pev 表示潜在蒸散、swvl1 表示表层土壤湿度单位也和直觉不符。读取时最常见的三个坑降水单位是「米」不是「毫米」、纬度是递减的、月度数据里可能带一个expver维度需要合并。import xarray as xr import numpy as np import pandas as pd def load_region(path, lon_range, lat_range): ds xr.open_dataset(path) # 纬度递减slice 必须从大到小否则裁出来是空的 ds ds.sel(longitudeslice(lon_range[0], lon_range[1]), latitudeslice(lat_range[1], lat_range[0])) # 部分产品带 expver 维度取第一个非全空切片合并 if expver in ds.dims: ds ds.isel(expver0) # 降水累积米 - 月总毫米 ds[precip_mm] ds[tp] * 1000.0 # 气温开尔文 - 摄氏度 ds[t2m_c] ds[t2m] - 273.15 return ds[[precip_mm, t2m_c, pev, swvl1]] ds load_region(data/raw/era5_monthly.nc, lon_range(100, 125), lat_range(30, 45)) print(ds.sizes) # {time: 360, latitude: 16, longitude: 26}load_region里的纬度切片顺序是最容易翻车的一行。ERA5 的 latitude 从 90 递减到 -90如果你按slice(30, 45)写xarray 会返回空数组而且不报错后面所有计算都在空数据上跑直到最后画图才发现一张白纸。单位换算也别省tp是累积米直接拿去算 SPI 会得到全零的荒谬结果。expver维度是 ERA5 特有的历史数据拼接标记不合并的话后续to_dataframe会多出一层索引特征表直接错位。裁剪完之后一般要重采样到月尺度并对齐时间标签。这里有个细节月度数据的时间戳通常标在月初或月中如果你把两个不同来源的数据比如再分析数据和观测站点直接 merge时间戳差半个月就会全部对不上。稳妥做法是统一resample(MS).mean()归到月初。3.2 干旱标签SPI 与 SPEI 的滑窗计算SPI标准化降水指数的逻辑是先算滑动累积降水比如 3 个月累计再把这个累积序列拟合成一个概率分布最后转成标准正态分位数。SPI 小于 -1 通常算中度干旱小于 -1.5 算重度。SPEI 把降水换成「降水减潜在蒸散」在升温背景下更能反映实际水分亏缺。工程上有两种实现路子参数法拟合 gamma 分布和经验法用秩序统计量。参数法在样本量小时拟合不稳定我一般用经验法因为它不依赖分布假设结果更可复现。from scipy.stats import norm import pandas as pd def spi_empirical(precip: pd.Series, scale: int 3) - pd.Series: 基于经验频率的 SPI 计算。 precip: 月度降水序列DatetimeIndex单位毫米 scale: 滑动窗口长度月常用 3 / 6 / 12 acc precip.rolling(scale, min_periodsscale).sum() valid acc.dropna() n valid.size if n 30: # 样本太少别硬算 raise ValueError(序列不足 30 个月SPI 不可靠) ranks valid.rank(methodaverage) # Gringorten 绘点位置公式避免概率取到 0 或 1 prob (ranks - 0.44) / (n 0.12) spi pd.Series(norm.ppf(prob), indexvalid.index) return spi.reindex(acc.index) # 补回前 scale-1 个 NaN spi3 spi_empirical(df[precip_mm], scale3) label (spi3.shift(-1) -1.0).astype(Int64) # 预测下个月是否干旱spi_empirical有几个参数必须盯着。scale3是业务上最常用的「季度干旱」定义改成 6 或 12 表示看更长期的干旱但窗口越长开头丢的 NaN 越多。min_periodsscale保证窗口不满时不输出半截值避免前期出现虚高的累积量。(ranks - 0.44) / (n 0.12)是 Gringorten 公式用 0.44 和 0.12 这两个常数是为了让极端值的概率不落在 0 和 1 上——norm.ppf(0)是负无穷会让整列数据被污染成 inf。最关键的一行是spi3.shift(-1)。这是「用当月信息预测下月状态」的标准做法shift 的方向反了就变成用未来数据预测现在属于典型的数据泄漏模型指标会好看得离谱上线就废。这个坑我在项目里踩过不止一次后面避坑章节还会再提。3.3 特征工程滞后项、滚动统计与空间邻域特征构造是这套 Notebook 里最能体现作者经验的部分。气象预测的特征大体分三类时间滞后、窗口统计、空间邻域。def build_features(df: pd.DataFrame, lags(1, 2, 3, 6, 12)) - pd.DataFrame: df 需含列precip_mm, t2m_c, swvl1索引为月度 DatetimeIndex feats pd.DataFrame(indexdf.index) for col in [precip_mm, t2m_c, swvl1]: for lag in lags: feats[f{col}_lag{lag}] df[col].shift(lag) # 滚动统计3/6/12 月累计降水和均值气温 for w in (3, 6, 12): feats[fprecip_roll{w}_sum] df[precip_mm].rolling(w).sum() feats[ft2m_roll{w}_mean] df[t2m_c].rolling(w).mean() # 季节项干旱有强季节性必须显式给模型 feats[month_sin] np.sin(2 * np.pi * df.index.month / 12) feats[month_cos] np.cos(2 * np.pi * df.index.month / 12) # 水分亏缺的物理量降水减蒸散 feats[water_deficit] df[precip_mm] - df[pev] * 1000 return feats逐段说明。滞后项lag1到lag12覆盖了短期记忆和年际记忆但滞后阶数不是越多越好十二阶对三百多个月的样本来说已经接近过拟合边缘树模型会自动筛线性模型就需要自己做相关性过滤。滚动统计里的roll12_sum实际上和年降水高度共线但它对识别「长期累积亏缺」有用值得保留。季节项用正余弦编码而不是 one-hot原因是正余弦保留了月份之间的连续性——十二月和一月在 one-hot 里完全无关在正余弦里距离很近这对干旱这种强季节现象更合理。water_deficit是物理驱动的特征不做标准化直接用会有量纲问题所以下一步必然要接一个标准化器。我一般把它包进 sklearn 的 Pipeline 里而不是在特征表阶段就标准化因为标准化参数必须只用训练集拟合否则又是一次泄漏。空间邻域特征在 Notebook 里通常是可选的把目标格点周围一圈格点的 SPI 均值也作为输入因为干旱有空间连片性。做法是用xarray.rolling(latitude3, longitude3).mean()做一次滑动平均再展平这一步在网格维度上操作内存占用不小建议先裁剪到研究区再做。4. 建模与评估把基线跑成可复现的数字数据备好之后建模部分反而简单了。这套 Notebook 的做法是先建立几个基线模型再逐步加复杂度每一步都记指标方便判断新特征到底有没有用。4.1 时间切分别用 train_test_split 随机打乱这是最容易犯也最致命的错误。气象时间序列有强自相关随机切分会让相邻月份同时出现在训练集和测试集里模型只要记住「上个月什么状态」就能拿到很高的分数但换成真实的未来预测立刻崩掉。import pandas as pd from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler from sklearn.linear_model import Ridge def time_split(feats, label, train_end2015-12-31, val_end2018-12-31): 按时间顺序切分保证训练集永远在验证/测试集之前 X_tr feats.loc[:train_end] y_tr label.loc[:train_end] X_va feats.loc[train_end:val_end].iloc[1:] y_va label.loc[train_end:val_end].iloc[1:] X_te feats.loc[val_end:].iloc[1:] y_te label.loc[val_end:].iloc[1:] return (X_tr, y_tr), (X_va, y_va), (X_te, y_te) (X_tr, y_tr), (X_va, y_va), (X_te, y_te) time_split(feats, label) print(X_tr.shape, X_va.shape, X_te.shape) # 例如 (336, 40) (36, 40) (24, 40)切分点选在 2015 年底和 2018 年底不是随便定的。留出三年验证集够跑一轮超参搜索留出两年测试集能覆盖至少一次完整的干旱事件周期。如果你的研究区在这几年里没有发生过显著干旱测试集的召回率会不可信这时候要么换切分点要么补数据。注意.iloc[1:]这个细节用loc切片时两端都是闭区间切分点那个月会同时出现在训练集和验证集里虽然只差一个月但在强自相关的数据上足以让验证指标虚高一点点。这种小泄漏最难查因为指标看起来「只是有点好」不会好到让你警觉。4.2 模型对比从 Ridge 到梯度提升树基线顺序建议是 Ridge → RandomForest → XGBoost/LightGBM。Ridge 用来验证特征和目标之间有没有线性关系如果 Ridge 的 R² 已经是 0.6说明任务本身不难后面的树模型提升有限如果 Ridge 接近 0 而树模型突然跳到 0.7那多半是某个特征偷偷泄漏了信息。from sklearn.ensemble import RandomForestRegressor from xgboost import XGBRegressor from sklearn.metrics import mean_squared_error, r2_score import numpy as np models { ridge: Pipeline([(sc, StandardScaler()), (m, Ridge(alpha1.0))]), rf: RandomForestRegressor(n_estimators400, max_depth8, min_samples_leaf5, random_state42), xgb: XGBRegressor(n_estimators500, learning_rate0.03, max_depth4, subsample0.8, colsample_bytree0.8, random_state42), } for name, model in models.items(): model.fit(X_tr, y_tr) pred model.predict(X_va) print(f{name:6s} RMSE{np.sqrt(mean_squared_error(y_va, pred)):.3f} fR2{r2_score(y_va, pred):.3f})参数选择上max_depth压到 4 到 8 是有意的。气候特征之间高度共线树太深会沿着噪声分裂训练集拟合得漂亮、验证集一塌糊涂。subsample0.8和colsample_bytree0.8是给样本和特征都加随机性在小样本场景下比调学习率有效。random_state固定住不然每次跑出来的指标都不一样没法判断改动是好是坏。Pipeline里包StandardScaler的意义在于标准化参数只从训练集学。如果你在特征表阶段就全局标准化验证集和测试集的均值方差信息会渗进训练过程这属于隐蔽的泄漏指标通常只虚高百分之几很难肉眼发现。4.3 评估指标RMSE 之外必须看事件命中率回归任务只看 RMSE 和 R² 会错过真正关心的东西。干旱预警的业务价值在于「该报的时候报出来」所以至少要补三个指标事件命中率召回、误报率、以及提前量。指标计算方式关注点RMSE预测 SPI 与真实 SPI 的均方根误差整体拟合程度R²解释方差比例相对气候态均值的提升事件召回真实 SPI -1 中被预测出的比例漏报代价高时优先看误报率预测为干旱但实际未发生的比例预警可信度提前量首次命中事件时距事件发生还有几个月预警可用性def event_scores(y_true, y_pred, thr-1.0): hit ((y_true thr) (y_pred thr)).sum() n_event (y_true thr).sum() n_pred (y_pred thr).sum() recall hit / max(n_event, 1) precision hit / max(n_pred, 1) return {recall: round(recall, 3), precision: round(precision, 3), miss: int(n_event - hit), false_alarm: int(n_pred - hit)} print(event_scores(y_te.values, models[xgb].predict(X_te)))注意阈值thr-1.0是业务的定义不是模型的一部分。你可以把模型输出的连续 SPI 直接发布让下游自己卡阈值这样更灵活。很多项目的做法是两者都出回归值给科研用二分类标签给预警系统用。还有个容易被忽略的验证角度分季节看指标。干旱在夏季和冬季的可预测性差别很大因为冬季主要受大尺度环流控制、夏季更多受局地土壤湿度记忆影响。如果整体 R² 是 0.55拆开看可能夏季 0.7、冬季 0.3那模型的实际适用窗口就很清楚了。5. 避坑与排查四个让结果悄悄作废的细节跑通和跑对是两件事。下面这几条都是那种「不报错、指标还挺好看」的问题属于最难发现的一类。5.1 数据与标签侧的坑现象模型 R² 高达 0.95换成下一年数据立刻掉到 0.1。原因特征里混进了未来信息。最常见的来源是滚动窗口没加shiftrolling(3).sum()默认包含当前月如果你用它预测当月标签等于把答案喂给了模型。另一个来源是全局标准化。解决所有滚动统计在入模前统一shift(1)标准化一律放进 Pipeline用fit_transform只作用于训练集验证测试集只transform。现象SPI 序列大量为 NaN或者开头几十个月全是空值。原因滑动窗口min_periods设成了窗口长度而序列开头本来就不满窗口或者经验频率计算里样本量不足 30norm.ppf在极端分位上返回了 inf。解决先检查序列长度和缺失情况再算指数SPI 至少需要 30 年月度数据才稳定如果数据不够改用更短的窗口或者换成百分位阈值法别硬算标准正态分位数。现象不同格点算出来的 SPI 量级差异很大同一场干旱在图上东边红西边蓝。原因SPI 是相对当地的标准化指标每个格点用自己历史序列拟合分布本来就该如此。但如果某个格点历史数据缺失多、有效样本少拟合出来的分位数会失真。解决对每个格点单独检查有效样本数低于阈值的格点标记为无效而不是强行输出出图时把这些格点遮掉别让噪声冒充信号。5.2 建模与评估侧的坑现象训练集和验证集指标都很好但预测出的干旱月份在时间上整体偏移了一两个月。原因标签的shift方向和特征的时间对齐错了。用当月特征预测下月标签应该是label.shift(-1)写反成shift(1)就变成用未来预测过去模型在训练集上依然能学得很好因为特征和标签的相对关系是一致的只是整个时间轴反了。解决拿一个已知的干旱年份做单点验证把预测值和真实 SPI 按月份并排打印出来看峰值是不是对齐。这一步比看任何指标都直接。现象换了个研究区同一套代码跑出来的召回率是 0。原因新区域的气候态和原区域差别太大SPI 阈值 -1 在新区域可能几十年都不触发或者某个特征在新区域根本不存在比如土壤湿度层深度不同。解决换区域后第一件事是画标签的分布直方图和事件年份时间线确认这个区域确实有干旱事件再谈模型特征表按列检查缺失率和方差方差接近零的列直接删掉。提示每次改动数据管道后固定跑一遍「标签分布 事件年份 单格点时间序列对照」这三张图比盯着指标数字可靠得多。6. 二次开发换数据源、加模型、把 Notebook 变成能定期跑的任务这套 Notebook 的扩展点很清晰按投入产出比排序我建议从这三处下手。第一处是换数据源。所有读取逻辑集中在01_load_and_clip.ipynb和src/utils.py只要新数据能转成带time / latitude / longitude三个维度的 xarray Dataset后面的流程一行都不用改。替换时唯一要小心的是变量名和单位映射——我一般会在src/utils.py里加一张映射表把不同来源的变量名统一成precip_mm / t2m_c / pev / swvl1四个标准名这样下游 Notebook 完全无感。VAR_MAP { era5: {tp: (precip_mm, 1000.0), t2m: (t2m_c, 1.0)}, custom: {precipitation: (precip_mm, 1.0), temp: (t2m_c, 1.0)}, }加模型则更简单04_train_baseline.ipynb里的models字典加一项就完事但要注意新模型必须走同一个time_split和同一个Pipeline别自己另起一套切分逻辑否则指标之间没法横向比较。如果想试序列模型建议先用sklearn的HistGradientBoostingRegressor做过渡它对缺失值友好、训练快效果通常不输调过参的 LSTM能用它解决的问题就不必上深度模型。第三处是把 Notebook 变成脚本。做法是把每个 Notebook 里的核心单元格抽到src/下的函数里用papermill或nbconvert定时执行或者干脆写一个run_pipeline.py按顺序调用。切换成脚本后会暴露一批 Notebook 里被隐藏的问题——最典型的是路径硬编码和中间变量依赖上一步的全局状态这时候把路径统一走pathlib加配置文件就基本能跑通了。验证改造是否成功我有个固定习惯改完之后用同一份历史数据从头跑一遍把关键指标RMSE、召回、事件年份记进一个 csv和改造前的记录逐行对比。指标差超过 5% 就必须回去查是哪一步变了哪怕只是换了个随机种子。从那以后我每次动数据管道哪怕只改一个shift的方向都强制走一遍「标签时间线对照 指标回归对比」这两步吃过一次亏就知道气候类项目里最贵的不是算力是你以为跑对了的那两周。希望帮到你。本文还有配套的精品资源点击获取
