做干旱监测的人对标准化降水指数（SPI）应该都不陌生。这个指数用一段降水序列拟合概率分布，再映射到标准正态空间，最后得到一个正负浮点值，正值偏湿、负值偏旱，用来判断某个区域某个时段的水分异常非常直观。但真正动手算过 SPI 的人都知道，麻烦不在于公式，而在于“降水数据从哪来、怎么整理成可计算的序列”。
我这次要聊的工作流是：通过 GEE 的 Python 接口取降水栅格数据，在本地用 climate_indices 包批量计算 SPI，最后做可视化。很多教程把这两步混在一起讲，实际上 GEE 只负责数据供给，SPI 的计算全在本地 Python 环境里完成。这样拆开有一个好处：换数据源、换区域、换时间尺度都只改一小段代码，不用动核心计算逻辑。
下面按我实际操作时的顺序展开，从数据选型到导出、从 SPI 参数设置到画图，最后复盘几个踩过的坑。这篇文章适合正在做干旱监测、水文气象分析，或者刚接触遥感降水产品的研究生和从业者，照着自己的需求替换参数就能跑通。
1. 降水数据源的选择：GEE 里三款主流产品的取舍
1.1 为什么把取数这一步放在 GEE
以前做站点级别的 SPI 分析，最头疼的是降水数据整理：要么从气象站要逐月观测值，要么下载再分析资料再转格式，数据质量参差不齐，时间长度也不统一。用 GEE 取数之后，这些问题被大幅简化——所有降水栅格都在同一个平台里，时间范围、空间范围可以靠代码直接控制，导出结果又是一份规整的表格或栅格，不用手工拼接。
尤其做区域尺度（比如一个流域、一个省、一块几万平方公里的样区）的 SPI 时，GEE 的优势更明显。你不需要在本地维护一堆历史栅格文件，只需要在代码里写清楚“从哪一天到哪一天、哪一块区域、哪个数据集”，剩下的交给人家的服务器去算。算完把结果拉到本地，进入 Python 分析环节。
需要说明的是，GEE 取到的降水是栅格产品的估值，不是站点实测。这意味着它适合做空间上的连续分析，但如果你的研究区有密集的站点数据，还是以站点资料为准，栅格产品只做参考或填补空白。
1.2 CHIRPS、ERA5-Land、PERSIANN-CDR 怎么选
我在 GEE 里常用的日尺度降水产品主要有三个，各有各的脾气。
| 产品 | GEE 集合 ID | 空间分辨率 | 起始年份 | 原始单位 | 适合场景 |
|---|---|---|---|---|---|
| CHIRPS | UCSB-CHG/CHIRPS/DAILY | 约 0.05 度（5.5 公里） | 1981 年 | mm/day | 站点稀疏区、干旱监测、季风区 |
| ERA5-Land | ECMWF/ERA5_LAND/DAILY_AGGR | 约 0.1 度（11 公里） | 1950 年 | m（日累计总量） | 长时序分析、全球覆盖、再分析研究 |
| PERSIANN-CDR | PERSIANN-CDR | 约 0.25 度（25 公里） | 1983 年 | mm/day | 长时序粗尺度趋势分析 |
CHIRPS 是“站点观测 + 卫星反演 + 地形修正”融合出来的产品，时空连续性在干旱监测里表现稳定，也是很多国际干旱监测系统的默认输入之一。它的分辨率是三款里最高的，做中小区域（比如几万平方公里的流域）很适合。代价是站点较少的地区，比如某些内陆高山区，CHIRPS 可能会被低估或高估，最好和当地站做一期对比再定。
ERA5-Land 是再分析产品，它的优势是时间跨度长，从 1950 年就有逐日数据，做 40 年以上的长时序 SPI 非常合适。但它的原始降水单位不是毫米，而是米，导出时如果忘了换算，后面算出来的 SPI 数值会整体漂移，这一点我放在最后的踩坑清单里专门说。
PERSIANN-CDR 分辨率最低，一般我不做精细区域分析用，只有在需要 1983 年至今的全球尺度长序评估时才会考虑它。实际上它和 CHIRPS 在月尺度上的趋势很接近，但空间细节差不少。
1.3 我的选型建议
如果你只需要一套数据，优先 CHIRPS。原因是 SPI 本身是一种“标准化相对指数”，它看重的是降水序列的统计分布，而不是绝对降水量，CHIRPS 的几十年连续序列恰好能满足这个需求。
如果你的研究区在高纬度或者需要 1981 年之前的历史时段，ERA5-Land 是更稳的选择，但务必先把单位从米换算成毫米。如果论文或报告要求多源交叉验证，那就把 CHIRPS 和 ERA5-Land 都跑一遍，看看 SPI 的时空格局是否一致；如果两者在关键旱季上的结论差异很大，那说明该区域降水产品不确定性较高，需要用站网数据进一步仲裁。
2. 从影像到一维时间序列：GEE 导出这一步别卡壳
2.1 定义研究区和影像集
我建议在本地用 Python 初始化 GEE 环境，直接在 Notebook 里操作。首先定义研究区，我这里用一个虚构的矩形区域做演示，坐标大致是东经 100 到 104 度、北纬 30 到 34 度，大约 4 度见方：
import ee # 如果是从零开始，先初始化 # ee.Authenticate() ee.Initialize() roi = ee.Geometry.Polygon([ [100.0, 30.0], [104.0, 30.0], [104.0, 34.0], [100.0, 34.0], ]) chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY").select("precip")这里有个细节：我在加载影像集时就用了 select 只保留降水波段。单看这一步无所谓，但如果你同时加载多个波段，后面做 reduceRegion 的时候会多出一堆用不上的数据，白白增加运行时间和导出体积。尽早瘦身是 GEE 取数的一个好习惯。
2.2 按自然月聚合的代码骨架
SPI 一般基于月尺度降水。CHIRPS 的原始数据是逐日 mm/day，所以一个自然月的月总量，就是把该月各天的降水值累加。这里必须按日历月来聚合，不能简单地用滚动 30 天，否则 SP 序列的相位会偏。
我常用的做法是生成一个“年 × 月”的编号序列，比如 198101、198102……然后在 GEE 端逐月做 reduceRegion。一个可用于中小区域面均值的示例代码如下：
def monthly_mean(ym): year = ym.divide(100).floor() month = ym.mod(100) start = ee.Date.fromYMD(year, month, 1) end = start.advance(1, "month") img = chirps.filterDate(start, end).sum() mean = img.reduceRegion( reducer=ee.Reducer.mean(), geometry=roi, scale=5500, ).get("precip") return ee.Feature(None, {"date": start, "precip": mean}) ym_list = ee.List.sequence(198101, 202312, 1) fc = ee.FeatureCollection(ym_list.map(monthly_mean))注意这里的 scale 参数。CHIRPS 分辨率约为 0.05 度，对应赤道附近约 5500 米。GEE 的 reduceRegion 允许你指定一个比原始分辨率略大的尺度来做区域平均，但不要给一个明显不合理的值，比如 100 米，那会触发大量不必要的计算。5500 米对 CHIRPS 来说比较合适。
2.3 导出到本地后的数据形态
上面的 fc 是一份 FeatureCollection，每个要素包含 date 和 precip 两个字段。导出到云端硬盘再下载，用标准的 export 函数就行：
task = ee.batch.Export.table.toDrive( collection=fc, description="chirps_monthly_roi", folder="GEE_EXPORT", fileNamePrefix="chirps_monthly", fileFormat="CSV", ) task.start()导出的 CSV 大致长这样：
system:index,date,precip 0,1981-01-01T00:00:00,7.83 1,1981-02-01T00:00:00,4.12 ...在本地用 pandas 读进来时，最好先把时间字段统一成“年-月”的格式，避免后续处理时区的问题：
import pandas as pd df = pd.read_csv("chirps_monthly.csv") df["month"] = pd.to_datetime(df["date"], utc=True).dt.strftime("%Y-%m") df = df.set_index("month") df.index = pd.PeriodIndex(df.index, freq="M") monthly = df["precip"].astype(float)这里把日期统一成 PeriodIndex 之后，后面计算 SPI 只需要取.values和一个起始年份，逻辑会非常干净。
2.4 单点直接取 vs 面均值导出
如果你研究区是一个点，或者只有几个站点位置，不必导出 CSV，直接用 sample 就能把像元值拉回本地。比如：
point = ee.Geometry.Point([102.0, 32.0]) def monthly_mean_at_point(ym): year = ym.divide(100).floor() month = ym.mod(100) start = ee.Date.fromYMD(year, month, 1) end = start.advance(1, "month") img = chirps.filterDate(start, end).sum() value = img.reduceRegion( reducer=ee.Reducer.mean(), geometry=point, scale=5500, ).get("precip") return ee.Feature(None, {"date": start, "precip": value})不过说实话，除非你只需要一两个点，否则还是做成面均值的 FeatureCollection 再一次性导出更省事。单个点一条条拉回本地的写法在 GEE 里会有大量往返请求，点一多就容易踩到计算配额限制。
3. Climate Indices 包的 SPI 计算：参数、原理与版本差异
3.1 安装与新旧版 API 的坑
这一步的主角是 climate_indices 包，pip 直接安装：
pip install climate-indices这个包的版本升级有点让人头疼。老版本（0.2.x 时期）的写法是直接用 indices 模块里的函数：
from climate_indices import indices spi_values = indices.spi( precip, timescale=3, distribution=indices.Distribution.gamma, data_start_year=1981, calibration_year_start=1981, calibration_year_end=2010, )新版本（0.3 以上）改成了类接口：
from climate_indices import ClimateIndices ci = ClimateIndices() spi_values = ci.spi( precip=precip, frequency="monthly", timescale=3, distribution="gamma", data_start_year=1981, calibration_year_start=1981, calibration_year_end=2010, )如果你在网上搜到的是旧写法，而本地装的是新版本，import 或者调用时会直接报错。我的建议是装好后先print(climate_indices.__version__)，再决定用哪种 API。下文统一用新版类接口。
3.2 核心参数逐个说明
把月度降水序列读成 numpy 数组之后，调用 SPI 的主代码是这样：
import numpy as np from climate_indices import ClimateIndices precip = monthly.values.astype(float) start_year = int(monthly.index[0].year) end_year = int(monthly.index[-1].year) ci = ClimateIndices() spi3 = ci.spi( precip=precip, frequency="monthly", timescale=3, distribution="gamma", data_start_year=start_year, calibration_year_start=start_year, calibration_year_end=2010, )逐项拆开看：
precip必须是一维数组，按时间升序排列，长度要足够覆盖校准期。如果序列里有缺失值，最好先补全再传入，否则结果里会出现意想不到的空洞。frequency填"monthly"或"daily"。注意这个参数描述的是输入序列的采样频率，不是你要算的时间尺度。你输入的是月序列，就填月，输入日序列就填日。timescale是累积窗口。timescale=3得到 SPI-3，表示最近 3 个月的累积降水在历史同期分布中的位置。1、3、6、12 是最常用的四档，分别对应短期土壤湿度异常、季节性水分异常和长期水文干旱。distribution填"gamma"或"pearson3"。SPI 的经典定义用的是伽马分布，这也是多数国际用户的默认选项。某些国家干旱监测系统偏好 Pearson III 型分布，两种在绝大多数情形下结果差异很小。data_start_year是输入序列的起始年份，包内部会用它来定位每个值的时间位置，所以必须填准确，尤其是当你的序列不是从 1 月开始的时候。calibration_year_start和calibration_year_end是分布拟合的校准期。常规做法是取气候标准期，比如 30 年。校准期越长，拟合出的分布参数越稳定，后面的 SPI 值也就越可信。
返回的数组长度和输入一样，但前几个值会是 NaN。具体来说，SPI-3 的前 3 个值算不出来，因为累积窗口还没填满。这不是 bug，后面画图的时候可以用同样的 x 轴索引直接画，让 NaN 自然断掉即可。
3.3 SPI 到底是怎么算出来的
为什么一个看起来只是“降水标准化”的指数，能成为干旱监测的标配？因为它的计算过程把气候背景差异抹掉了：某个区域年均降水量 300 毫米，另一个区域年均 1200 毫米，直接用绝对降水量没法比较，但转换成 SPI 之后，两地的干湿程度就处于同一尺度上。
计算过程大致分四步：
- 先把最近
timescale个月的降水累加。比如 SPI-3 在 5 月，就用 3 月、4 月、5 月三个月的降水之和。 - 对累计降水序列拟合伽马分布，得到形状参数和尺度参数。
- 把每个累计值转成累积概率。这里有个关键细节：当月降水存在大量零值时，伽马分布在零处没有定义。通常的做法是构造一个混合累积概率，把“零降水出现的经验频率”单独拿出来，非零部分再用伽马分布拟合。
- 最后把累积概率映射到标准正态分布的分位数上。比如累计概率是 0.1，对应的标准正态分位数约是 -1.28，这个值就是 SPI。
所以 SPI 正值表示降水偏多，负值表示偏少，数值离 0 越远说明异常程度越高。常用分级如下：
| SPI 区间 | 分类 |
|---|---|
| 2.0 及以上 | 极度湿润 |
| 1.5 到 1.99 | 非常湿润 |
| 1.0 到 1.49 | 中等湿润 |
| -0.99 到 0.99 | 接近正常 |
| -1.0 到 -1.49 | 中等干旱 |
| -1.5 到 -1.99 | 严重干旱 |
| -2.0 及以下 | 极端干旱 |
这也是为什么很多人喜欢在图上画 -1、-1.5、-2 三条虚线：它们正好对应干旱从启动到加强的几道门槛。
4. 可视化落地的两条路线：时间曲线与空间分布图
4.1 时间序列图：参考线、填充和双轴
算完 SPI 之后，最常用的可视化就是时间序列图。我一般用 matplotlib 画，三步走：画参考线、填充正负区域、叠加降水量做双轴。
import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(12, 5)) x = monthly.index.to_timestamp() # SPI 曲线 ax.plot(x, spi3, color="#2166ac", lw=1.2) ax.axhline(0, color="#333333", lw=0.8) for level, color in [(-1, "#8c8c8c"), (-1.5, "#f39c12"), (-2, "#c0392b")]: ax.axhline(level, ls="--", lw=0.8, color=color) ax.fill_between( x, spi3, 0, where=(spi3 >= 0), color="#2166ac", alpha=0.25 ) ax.fill_between( x, spi3, 0, where=(spi3 < 0), color="#b2182b", alpha=0.25 ) ax.set_ylabel("SPI-3")填充颜色的意义在于，一眼扫过去就能看到哪些时段持续偏旱、哪些时段偏湿，比单看折线直观很多。如果你想同时呈现降水量的变化，可以在同一个图上加一个次级坐标轴，把月降水画成柱状：
ax2 = ax.twinx() ax2.bar(x, monthly.values, color="#d9d9d9", alpha=0.6, label="Precip") ax2.set_ylabel("Precipitation (mm)")这时 SPI 曲线叠加在灰色柱子上，降水异常和 SPI 波动的对应关系非常清楚。类似的图在干旱监测报告里很常见，做法不复杂，但信息密度很高。
4.2 空间分布图：逐像素算 SPI 的正确姿势
时间序列图只适合呈现一个点或一个区域均值。如果想看某年某月的 SPI 在空间上怎么分布，就得扩展到每个像素。这也是“GEE python 算 SPI”最容易让人误入歧途的地方：GEE 服务器端并不适合直接做伽马分布拟合和分位数映射，因为你没法方便地在服务器端代码里调用 climate_indices 包。我试过几种方案，最终的落地做法是：
- 在 GEE 端把整个区域按月聚合，导出月降水栅格序列，格式可以选 GeoTIFF，也可以整理成 NetCDF。
- 本地用 xarray 打开，得到一个三维数组（time, latitude, longitude）。
- 自定义一个包装函数，调用 climate_indices 里的 SPI 函数。
- 用
xr.apply_ufunc沿着时间维逐像素点应用这个函数。
核心思路如下：
import xarray as xr ds = xr.open_dataset("monthly_precip.nc") precip_da = ds["precip"] def calc_spi(ts): arr = ci.spi( precip=ts, frequency="monthly", timescale=6, distribution="gamma", data_start_year=1981, calibration_year_start=1981, calibration_year_end=2010, ) return arr spi_da = xr.apply_ufunc( calc_spi, precip_da, input_core_dims=[["time"]], output_core_dims=[["time"]], vectorize=True, )vectorize=True会把这个函数逐像素跑一遍，代价是速度比较慢。如果区域很大，先在小范围测试，再放开到全区域。画图的时候可以用 cartopy 叠加海岸线、河流等地理要素，这是目前最稳妥的组合。
这里要提醒一句：导出的月降水栅格序列在时间维上必须是完整的，不能缺月份。缺一个月份，累积窗口就会错位，SPI 结果会在该时段附近出现畸形跳变。检查时间维的长度是最基本的一步。
5. 实操复盘：六个容易翻车的细节
5.1 数据预处理相关的三个坑
第一个坑是单位。CHIRPS 的原始值是 mm/day，月度累加之后得到 mm/month，这个量级是对的。但 ERA5-Land 的 total_precipitation 单位是米，换算成毫米需要乘 1000。如果忘记这一步，SPI 的符号可能不会变，因为标准化处理对整体缩放不敏感，但校准期分布参数会受影响，最终 SPI 数值和正确结果有明显偏差。这种错误很隐蔽，因为它不会直接报错。
第二个坑是时间聚合边界。CHIRPS 的逐日影像按 UTC 日期归档，聚合月份时一定要用日历月的边界，也就是从当月 1 日开始到下月 1 日之前。如果你在中途用了滚动 30 天窗口，SPI 序列会和真实月尺度错位，尤其月末几天降水多的时候，差异会很明显。
第三个坑是缺失月份。某些极端干旱区域或高云覆盖区，月聚合可能会出现全 NaN 的像元。最好在导出之前先检查每个月份影像的数量，比如chirps.filterDate(start, end).size()，如果数量不足，需要决定是补 0 还是插值。直接把 NaN 丢给 climate_indices 包，大概率会在拟合阶段翻车。
5.2 包使用相关的两个坑
第四个坑是新旧版本 API 混用。很多人会在 Stack Overflow 上找到 0.2.x 时期的代码，粘贴后发现indices.Distribution不存在。其实解法很简单，升级到新版本后统一用ClimateIndices类接口；如果项目环境必须锁定老版本，那就要用老接口。最省心的办法是在项目 requirements 里固定版本号。
第五个坑是校准期参数选择。校准期太短，比如只有 5 年，SPI 的标准差会被样本量拉偏，导致后期所有值都偏大或偏小。至少要有 30 年连续降水序列，这是 WMO 推荐的气候标准期长度。如果序列本身只有十几年，那算出来的 SPI 只能算“相对某段时期的异常指数”，跨区域比较时要非常谨慎。
5.3 一个藏得最深的单位坑
第六个坑其实还是单位，但这次藏在 PERSIANN-CDR 这类产品里。这些产品的日降水量字段有时是 mm/day，有时是 mm/hr，取决于 GEE 集合里的波段说明。直接用之前一定要看image.getInfo()或者数据集的波段属性，确认单位换算系数。我就遇到过把 mm/hr 当天尺度直接用，结果月总量差了一个数量级的案例，当时整条 SPI 曲线看起来都在极旱区徘徊，排查了很久才发现根因在原始数据单位。
另外，SPI 计算结果的数组前几个值是 NaN，这是正常行为。有人会在这一步用dropna()把开头空值删掉，结果后面画图时 Excel 或 pandas 对齐错位，导致不同尺度的 SPI 曲线错开几个月。处理方式是保留原始索引，画图时让 NaN 自然断裂，不要轻易重排数组。
整个流程跑通之后，我的体会是：SPI 这个指数本身不算难，难的是把数据源、GEE 导出、本地包调用、可视化四段链路衔接好。每段链路单独看都很简单，但放在一起，任何一处的细节出错都会让你在下一步面对一个看起来“完全合理实则错误”的结果。建议拿到一个新数据集时，先找一个有站点降水记录的区域做一次验证，把 SPI-1 的曲线和当地实际旱情报道对照一下，确认整条链路没有系统偏差，再扩展到更大区域。这套流程稳定下来之后，换数据源、换时间尺度都只是改参数的事情。
