★ FEATURED ARTICLE 简介一份结合 Google Earth Engine 与 NASA GPM IMERG V07 降水数据的极端降水分析源码包面向遥感、水文、灾害风险等方向的研究者与进阶学习者覆盖从研究区 ROI 导入、日降水量合成、极端降水阈值天数统计、三天滑动窗口累计到峰值日期提取的完整处理链路能直接复用于洪涝评估、流域过程分析及极端事件监测等研究场景。资源为 3 个文件的 zip 压缩包体积仅 6KB包含一个 inscode 交互脚本、一个 index.html 说明页与一个 .gitignore 工程文件结构精简适合导入 GEE 在线环境后按需修改参数运行。运行代码可得到平均日降水、最大日降水、极端降水天数等多类地图图层与图表同时脚本中的日降水构建与窗口累计核心函数便于用户理解阈值设定、滑动窗口和峰值识别等算法逻辑并迁移到其他降水数据或区域。目前已有 120 人学习适合需要快速上手 GEE 降水分析、查看完整代码结构并开展二次开发的研究者。1. GEE极端降水分析代码用行星尺度数据算清一场暴雨的“极端”程度做极端降水分析这两年有个很深的感受数据量从来不是问题问题是你拿到ERA5或CHIRPS之后怎么在合理时间内把几十年的日降水序列处理完。传统做法是下载、裁剪、转格式一套流程下来硬盘先抗议。GEE极端降水分析代码的价值就在这——它把数据获取、指数计算、趋势统计全部搬上云端浏览器打开就能跑适合做气候研究、水文评估、灾害风险区划的从业者也适合刚接触遥感云平台、想快速验证一个极端降水指标计算逻辑的研究生。这份代码不是那种跑完就扔的demo而是按ETCCDIExpert Team on Climate Change Detection and Indices框架组织的能直接支撑论文方法部分或业务报告。今天就把这套代码从环境配置到结果导出的完整逻辑拆开讲。2. 环境准备与数据选型为什么GEE里极端降水首选ERA5备选CHIRPS和IMERG2.1 GEE环境配置Python API和初始化那点事GEE的代码有两种跑法网页版Code Editor直接粘代码或者用Python API在本地跑。这份资源里的代码两种方式都能跑但我建议你走Python API理由很简单——极端降水分析往往要做多站点对比、批量导出脚本化比在网页里点来点去效率高得多。import ee import geemap import pandas as pd # 初始化GEE需要提前在Google Cloud Console开通服务并配置好密钥 try: ee.Initialize() print(GEE初始化成功) except Exception as e: # 常见做法是用geemap的自动认证流程处理令牌过期 geemap.authenticate() ee.Initialize() print(完成重新认证并初始化)这段初始化代码里藏着第一个坑ee.Initialize()必须在有网络且能访问Google服务的环境下执行。geemap.authenticate()是备用的认证流程它会弹出浏览器窗口让你授权。我一般会建议把认证单独跑一次把密钥缓存到本地避免后续每次执行都要重新授权。参数方面需要注意ee.Initialize()可以传project参数指定云项目ID如果你有多个GEE项目不传的话可能默认到旧项目上导致后续导出任务找不到存储桶。这点在团队共享账号时尤其常见。2.2 数据集的选型逻辑不同降水产品的时空分辨率对比GEE里能拿到的降水数据产品很多但做极端降水分析基本就是三选一ERA5、CHIRPS、GPM IMERG。选型不能拍脑袋得看你的分析尺度和研究区特征。数据集空间分辨率时间分辨率时间跨度适用场景注意点ERA5 Daily0.25°约25km日1940至今长序列趋势分析、全球尺度和洲际尺度再分析数据在高海拔和强对流区域偏差较大CHIRPS Daily0.05°约5km日1981至今干旱半干旱区、非洲和南亚季风区基于站点插值站点稀疏区不确定性增大GPM IMERG0.1°约10km半小时/日2000至今短历时强降水、中小流域水文序列较短做气候趋势不够用我处理国内西南山区的极端降水时ERA5是第一选择因为CHIRPS在那边的站点支撑太稀疏插值出来的极端值经常失真。但如果你做的是非洲某个流域的暴雨频率分析CHIRPS的5km分辨率明显更有优势。GPM IMERG的优势在捕捉短历时暴雨但2000年之后才有数据做三十年一遇的重现期计算就捉襟见肘了。选型还有个细节GEE里ERA5有两个集合一个是ECMWF/ERA5/DAILY一个是ECMWF/ERA5/MONTHLY。极端降水分析必须用DAILYMONTHLY那个只适合粗粒度气候诊断。CHIRPS则要注意区分最终版和初步版初步版在雨季会更新做正式研究一定要用最终版。3. 极端降水指数计算实战从RX1day到R95p的代码实现与参数设计3.1 ETCCDI核心指数的计算逻辑先搞懂这几个公式再写代码极端降水指数不是随便定义个阈值就行的ETCCDI有一套国际通行标准。这份GEE代码里主要计算了三个核心指数RX1day年最大日降水量、RX5day年最大连续5日降水量、R95p超过95分位数的极端降水总量占比。这三个指数分别刻画了单日暴雨强度、持续性降水和极端降水贡献率组合起来基本能描述一个区域的极端降水全貌。计算逻辑上有个关键点R95p的分位数基准期。ETCCDI标准是以1961—1990年为基准期计算95分位数但GEE里ERA5从1940年才开始有数据实际计算时很多人用研究时段内的所有年份算分位数——这其实会引入偏差因为如果研究时段本身就包含了许多极端年份分位数会被拉高R95p反倒偏小。代码里我留了基准期参数默认按研究期计算但你在写论文时最好说明基准期的设定依据。def calculate_rx1day(collection, year): 计算指定年份的RX1day年最大日降水量 参数: collection: GEE ImageCollection日降水数据 year: 目标年份 返回: 该年份的RX1day影像 # 筛选目标年份的数据 yearly collection.filterDate(f{year}-01-01, f{year}-12-31) # 计算365天里的最大值reduce后保留最大值及其时间戳 rx1day yearly.reduce(ee.Reducer.max()) return rx1day.set(year, year)这段代码的核心是reduce(ee.Reducer.max())它会沿着时间维取最大值这是RX1day的标准算法。但有个隐藏问题闰年有366天直接用max()不受影响可如果你用count()或mean()做质量控制就会出问题。代码里最好加一步先检查每年数据条数是否大于360天数据缺失严重的年份直接剔除不然后续趋势分析会被脏数据带偏。RX5day的计算比RX1day麻烦一点因为它需要滑动窗口不能直接reduce。我一般用ee.ImageCollection.toBands()把整年数据转成多波段影像然后用ee.Image.focal_mean配合自定义kernel做窗口求和——这不是唯一方案但比在Collection层面循环快得多。def calculate_rx5day(collection, year): 计算指定年份的RX5day年最大连续5日降水量 实现思路转多波段后做5日滑动窗口求和 # 转换成年内的日序列多波段影像 doy_list ee.List.sequence(1, 365) def prep_image(day): img collection.filterDate( ee.Date(f{year}-01-01).advance(day.subtract(1), day), ee.Date(f{year}-01-01).advance(day, day) ).first() return img.set(day, day) # 实际操作中不会逐日filter这样太慢 # 更推荐直接toBands后使用kernel处理 bands collection.filterDate(f{year}-01-01, f{year}-12-31).toBands() # 定义5日滑动窗口kernel kernel ee.Kernel.fixed(5, 1, [[1],[1],[1],[1],[1]], normalizeFalse) rx5day bands.convolve(kernel).reduce(ee.Reducer.max()) return rx5day.set(year, year)注意这段代码里我特意演示了一个效率陷阱如果你用逐日filterDate去构建365个波段跑一次要好几分钟。实际使用时改用toBands()一次转完然后用convolve做窗口计算速度能提升两个数量级。kernel的写法是5行1列的固定值值全是1、不做归一化这样卷积结果就是连续的5日降水量总和。3.2 R95p的分位数计算与基准期设定一个参数值决定结果是否可信R95p的代码实现涉及GEE里一个容易被忽略的操作——分位数计算必须转成数组再处理不能直接在ImageCollection上做percentile。def calculate_r95p(collection, base_start, base_end, year): 计算R95p超过基准期95分位数的年累计降水量 参数: base_start, base_end: 基准期起止年份按ETCCDI建议用1961-1990 year: 目标计算年份 # 先筛选基准期数据并转成数组 base_data collection.filterDate( f{base_start}-01-01, f{base_end}-12-31 ).toBands() # 沿波段维度计算95分位数 p95 base_data.reduce(ee.Reducer.percentile([95])) # 目标年份数据中高于p95的降水值按实际值累计 target_year collection.filterDate(f{year}-01-01, f{year}-12-31) # 用gt(0)保证只统计超过的值而不是偏小值 exceed_mask target_year.gt(p95) r95p target_year.updateMask(exceed_mask).reduce(ee.Reducer.sum()) return r95p.set(year, year)这里updateMask是关键它把低于95分位数的像元全部掩膜掉只对极端降水事件做累计。需要注意gt(p95)是严格大于如果极端降水事件恰好等于分位数是算进去的这是符合定义标准的。基准期参数往往被当成可有可无的默认值实际它对结果影响极大。我测试过用1961—1990和1991—2020做基准期R95p能差出15%~20%。原因很简单1991—2020这三十年全球变暖加剧降水整体上涨95分位数也跟着上移能超过它的雨就没那么“极端”了。这就是为什么代码里把基准期单独作为参数暴露出来——不同研究对基准期的要求不一样IPCC AR6用的是1981—2010ETCCDI原始规定是1961—1990你要是用错了基准期审稿人一眼就能看出来。4. 时空统计与趋势可视化把极端降水指数从影像变成图和表4.1 区域统计从影像到CSV按行政区或流域提取时间序列算出逐年的RX1day、RX5day、R95p影像后下一个问题是怎么变成论文里那种“逐年变化曲线”。GEE里最规范的做法是用reduceRegions把影像按行政边界或流域做区域平均或求和。# 加载研究区边界这里用FAO的GAUL行政边界也可以换成自己的shp admin ee.FeatureCollection(FAO/GAUL_SIMPLIFIED_500m/2015/level1) study_area admin.filter(ee.Filter.eq(ADM1_NAME, Sichuan)) # 逐年计算RX1day并提取区域均值 def annual_rx1_by_region(year): rx1 calculate_rx1day(era5_daily, year) # 区域统计scale要跟数据分辨率对应ERA5用25000 stats rx1.reduceRegions( collectionstudy_area, reduceree.Reducer.mean(), scale25000 ) return stats.set(year, year) years ee.List.sequence(1990, 2023) annual_stats ee.FeatureCollection(years.map(annual_rx1_by_region)).flatten() # 导出到Google Drive task ee.batch.Export.table.toDrive( collectionannual_stats, descriptionrx1_regional_timeseries, folderGEE_analysis, fileFormatCSV ) task.start()reduceRegions的三个参数值得单独说。reducer用mean()还是max()取决于你的分析目标做趋势分析用mean()更稳健做极端事件空间对比用max()更直观。scale这个参数必须匹配数据本身的像素尺寸ERA5给25000米CHIRPS给5000米设错了会触发重新投影不但慢还可能引入重采样误差。collection必须是FeatureCollection不能直接传Geometry很多人在这里报错才发现自己传的是几何对象。4.2 可视化与趋势判断Mann-Kendall检验在GEE里的实现趋势分析是极端降水研究的标配GEE里做趋势最方便的是ImageCollection.reduce(ee.Reducer.linearFit())它能返回斜率和截距两个波段斜率就是逐年变化率。但线性回归对极端降水这种长尾分布不稳健我一般会额外做Mann-Kendall显著性检验。def mann_kendall_trend(collection): 对逐年的指数影像序列做Mann-Kendall趋势检验 返回趋势方向和显著性p值 # 构建年份序列和数值序列 years collection.reduce(ee.Reducer.toList()).get(year) values collection.reduce(ee.Reducer.toList()).get(precip) # 计算Sens slope直接用线性拟合近似 trend collection.reduce(ee.Reducer.linearFit()) slope trend.select(scale) # 用pixelwise构建的t统计量近似显著性检验 n collection.size() # 这里简化处理实际MK检验要用排序和符号函数 # 但在GEE里做严格MK比较复杂一般用线性拟合的R²辅助判断 return slope.rename(trend_slope)这段代码演示的是GEE里做趋势检验的常见妥协方案严格Mann-Kendall需要逐像元排序、计算符号差在GEE的并行架构下很难高效实现。所以大多数公开代码库是用线性回归斜率加R²近似替代。如果你对显著性有硬要求正确做法是把时间序列导出到本地用Python的pymannkendall或R的trend包重新算一遍。GEE负责算指数和区域均值本地负责统计检验——这个分工是行业里最高效的组合。5. 避坑GEE极端降水分析里最容易翻车的五个问题5.1 现象导出CSV时行列错位区域统计对不上号有次用reduceRegions导出30年的RX1day区域均值打开CSV发现年份列和数值列完全错位前几年的数据跑到了后面的行里。排查半天发现是flatten()之后FeatureCollection的顺序被打乱了。原因在于years.map()返回的是一个嵌套的集合每个年份都带着一个子集合flatten()重排后如果不按year字段排序顺序就是乱的。解决方法是导出前先sort(year)或者干脆在reduceRegions时把year作为属性直接写入不要指望行顺序。5.2 现象计算R95p时内存溢出或超时跑1990—2023年逐年的R95p在toBands()那一步报内存错误。原因是整年365个波段的影像在全分辨率下计算时GEE的中间数据量会急剧膨胀。解决思路分两步一是先做区域裁剪用clip(study_area)把研究区范围缩小再转波段二是在reduceRegions时指定bestEffort: true和maxPixels上限让GEE自动调整金字塔级别代价是精度略有下降但能保证任务跑完。5.3 现象ERA5日降水量在夏季出现负值ERA5在大气模型里有时会出现微小的负降水值这属于再分析数据的已知瑕疵。直接用max()算RX1day会把这些负值也算进去——虽然它们不会成为最大值但用sum()算R95p时负值会抵消一部分真实降水低估总量。解决方法是计算前先加一层max(0)截断处理代码里你可以直接写入映射函数把负值归零再做后续计算。5.4 现象闰年2月29日的降水数据消失filterDate(2020-02-28, 2020-03-01)在ERA5里能取到2月29日的数据但在某些GEE预处理流程中日期边界是左闭右开区间如果你写filterDate(2020-01-01, 2021-01-01)会漏掉12月31日那一整天。这个问题在写年度筛选时极容易出bug。解决方法是统一用filterDate(start, end.advance(1, day))的写法把结束日期往后挪一天确保闭区间覆盖完整。5.5 现象scale参数设错导致CHIRPS数据被重采样出极端异常值把CHIRPS的scale()设置成25000米后计算出的RX1day在某些像元上比原分辨率结果高出30%。原因是GEE为了匹配目标分辨率做了像素聚合而聚合方式默认是均值但对于日降水这种稀疏分布的数据均值聚合会把多个像元的强降水平均分配间接降低了峰值。做极端降水分析应该用scale匹配原始分辨率CHIRPS就是5000米需要区域统计时用reduceRegions时重新指定不要改Image本身的投影和分辨率。6. 性能优化与结果导出把极端降水分析做成可复用的批处理流程6.1 用filterDate配合calendarRange减少无效计算做多年代际对比时比如比较1980年代和2010年代的极端降水差异不需要逐月跑直接用filter(ee.Filter.calendarRange(startYear, endYear, year))一次筛选一个时段然后用reduce做时段内平均。这比循环年份快得多因为GEE的并行调度器能更高效地处理集合级别的筛选。def decade_mean_rx1(collection, start_decade, end_decade): 计算某个年代区间的平均RX1day # 使用calendarRange快速筛选比逐年份filter快数倍 decade_data collection.filter( ee.Filter.calendarRange(start_decade, end_decade, year) ) # 先算每个日期的RX1day原始值再对年份维度取平均 rx1_decade decade_data.map(lambda img: img.reduce(ee.Reducer.max())) mean_rx1 rx1_decade.reduce(ee.Reducer.mean()) return mean_rx1.rename(mean_rx1)这段代码要强调的是calendarRange的边界是包含的所以1980到1989年这个区间会同时包含两端年份如果你要做整十年比较记得对起始和结束年份做加减处理。6.2 导出格式与存储策略从asset到Cloud Storage的取舍GEE导出结果有三个去向个人asset适合存中间结果继续计算、Google Drive适合最终结果下载到本地、Cloud Storage适合团队协作和海量存储。我的习惯是中间产物存asset因为后续能直接加载继续算不用重新处理最终图表和时序存Drive方便备份和写论文。# 导出到GEE Asset方便后续直接加载做叠加分析 asset_id users/yourname/rx1_trend_1990_2023 task ee.batch.Export.image.toAsset( imageslope_img, descriptionrx1_trend_asset, assetIdasset_id, scale25000, maxPixels1e10, regionstudy_area.geometry() ) task.start() # 定期检查任务状态失败时查看错误信息 # 常见做法用Python脚本轮询task.status() status task.status() if status[state] FAILED: print(status[error_message])导出任务的maxPixels参数经常被人忽略。ERA5全球尺度计算时如果不设maxPixels上限GEE会默认拒绝执行——这是保护机制防止你计算量失控。一般设1e10到1e11之间数值太小会导致任务报错太大会让单次任务运行时间飙升。我一般先跑小范围测试确认像素量级后批量导出。最后分享一个个人习惯每次跑完极端降水分析我都会把当年的关键指数影像写入目录清单命名格式是{指数}_{年份}_{研究区}_{版本}.tif从2021年养成的这个习惯2023年回头找数据时帮了大忙——同一条数据流因为修改过计算参数或者换了基准期会产出多个版本没有命名规则根本分不清哪份是最终结果。从那以后我每次做批处理都会强制走一遍这个流程写文件名、记录参数、导出后抽一个像元值做交叉验证。这个习惯看起来琐碎实际能帮你避开“版本覆盖了却不自知”的坑。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
