RClimDex极端气候指数计算：质量控制与bootstrap重采样

简介极端气候指数计算软件配套文档面向气候科研、气象业务及相关专业学习者聚焦极端气候事件的监测与指数化分析。压缩包仅含1个doc文件约600KB为RClimDex用户手册介绍该软件由加拿大环境部气象研究分支开发维护的背景、安装运行、数据质量控制与指数计算流程并附已知问题与缺陷报告指引。正文及附录覆盖ETCCDMI推荐的27项核心气候指数清单、输入数据格式、各指数定义以及阈值和基准期温度指数的计算方法与统计意义另含R在Windows环境下的常见问题解答便于读者理解软件功能边界、核对计算结果并排查数据格式问题。已有3839人学习下载适合需要查阅RClimDex使用说明、核对指数定义、排查数据格式问题或进行结果解释的研究者作为案头参考。1. 极端气候指数不是算不出来是算得不一致手里攥着一份六十年的逐日气温降水表想回答“这地方的霜冻日到底有没有在减少”最省事的做法是丢进表格软件拉一个条件计数。可只要换个人、换个起始年、换一种缺测处理方式同一个站能算出两条方向相反的趋势线。RClimDex 这款极端气候指数计算软件就是冲着这个不一致来的它把 ETCCDMI 推荐的 27 个核心极端气候指数固化成一套带强制质量控制的流程喂进一个纯文本日值文件吐出逐年指数序列和趋势图。它从表格环境迁到 R 上核心动机只有一个——百分位类温度指数需要 bootstrap 重采样表格软件做不了。手上有多站逐日数据、要出可横向对比指数序列的人值得把它拆开看一遍只想画条气温曲线的这套东西偏重了。2. 跑起来之前R 运行时、六列文本和三个输出目录这套工具的门槛不在算法在“喂进去的东西必须长得一模一样”。它不是一个装完就能双击的 exe而是一段 R 源码加上一层 TclTk 图形界面所有的交互都建立在“当前工作目录里有什么文件”这件事上。目录搞错、列数错一位、缺失值写成 -999界面会直接卡住或者弹一个语焉不详的报错。所以先把底座铺平再谈算指数。2.1 R 与 TclTk为什么二十年前的界面今天还能开RClimDex 依赖两样东西R 的基础环境以及 TclTk 图形库。手册成书年代的 R 版本号还停留在 1.8.x那个年代 TclTk 需要单独装从 R 1.9.0 之后Windows 版的默认安装包里已经带上了 TclTk所以你今天下载任意一个 4.x 的 Windows 安装包一路下一步走完图形库就齐了。Linux 下则建议用发行版仓库里的预编译包源码编译虽然也能过但 TclTk 的开发头文件经常缺会白折腾一轮。装完之后先做一件事确认library(tcltk)不报错。这一步比什么都重要因为 RClimDex 的主菜单、参数窗口、文件选择框全部由 TclTk 渲染缺了它source()执行到一半会抛出一个跟界面毫无关系的错误新手很容易误判成源码损坏。# 环境自检三行都通过再往下走 getRversion() # 建议 4.0 以上 capabilities(tcltk) # 必须为 TRUE library(tcltk) # 不报错即正常Windows 下还有一个老问题R 启动时会读一个叫HOME的环境变量用来定位用户目录。部分机器上这个变量没设source()用相对路径找文件时就会跑到系统目录里去报“找不到 rclimdex.r”。解决办法不复杂在系统环境变量里加一条HOME值指向你的工作目录即可。这是手册附录里专门拎出来讲的一条说明当年踩的人不少。2.2 输入文件的六列约定一行都不能多RClimDex 只认一种输入纯 ASCII 文本文件名形如站号.txt无表头空格或制表符分隔固定六列顺序是年、月、日、降水量、日最高温、日最低温。缺测统一写-99.9这是程序内部识别的哨兵值写-999、NA、空字符串都会被当成非法数据。列序字段单位约定1year四位年如 19612month1–12不补零3day1–31不补零4prcpmm缺测写 -99.95tmax℃缺测写 -99.96tmin℃缺测写 -99.9先用 shell 过一遍字段数比在界面里反复试错快得多。任何一行字段数不等于 6程序在读入阶段就会失败而报错信息只会告诉你“格式不对”不会告诉你第几行。# 逐行检查字段数把非 6 列的行连行号一起打出来 awk NF ! 6 {print NR: $0} 21946.txt | head -20 # 顺便看看有没有人用逗号做了分隔 grep -c , 21946.txt字段数对了还要看日期链是否连续、有没有重复日。RClimDex 内部按日索引缺一天不会报错但会让某些需要滑动窗口的指数比如连续干旱日数、生长期长度在跨年处算错位。用下面这段做一次全量体检raw - read.table(21946.txt, header FALSE, na.strings -99.9, col.names c(year,month,day,prcp,tmax,tmin)) # 1) 日期能否正常构造 d - as.Date(sprintf(%04d-%02d-%02d, raw$year, raw$month, raw$day)) cat(非法日期数:, sum(is.na(d)), \n) # 2) 重复日期 cat(重复日期数:, sum(duplicated(d)), \n) # 3) 缺失的日历日 full - seq(min(d, na.rm TRUE), max(d, na.rm TRUE), by day) miss - setdiff(as.character(full), as.character(d)) cat(日历缺口数:, length(miss), \n); head(miss, 10)na.strings -99.9这一句是关键它把哨兵值在读入阶段就转成 R 的NA后面的统计才不会把 -99.9 当成真实的低温参与平均。日期缺口如果集中在某几年八成是原始数据本身断档不用强行补齐让它在指数序列里体现为缺测年份反而更诚实。2.3 加载源码与工作目录source()之后发生了什么手册给的加载方式是source(rclimdex.r)不是library()。这点经常被搞混——早期分发形态就是一份单文件源码靠source()注入当前会话执行完主菜单窗口就弹出来了。后来也有打包成 R 包的版本用install.packages()加library()加载两种方式对应不同分发渠道拿到手先看对方给的是.r文件还是包名不要混着用。# 方式一本地源码最原始也最稳 setwd(D:/climdex) # 必须先切到 rclimdex.r 所在目录 source(rclimdex.r) # 成功后弹出主菜单 # 方式二已打包版本 # install.packages(RClimDex) # library(RClimDex)setwd()这一步不是可选项。RClimDex 的所有输出都写在“当前工作目录”下它不会问你存哪儿而是自动建三个子目录log/质量控制阶段的四个 CSV 和四张 PDF 时序图indices/每个指数一个 CSV命名规则是站号_指数名.cvsplots/每个指数一张 JPEG 趋势图文件名与 indices 里一一对应。注意indices目录里文件的扩展名是.cvs而不是.csv这是程序里写死的拼写不是笔误也不是你下载坏了。双击打不开是正常的用文本编辑器打开或者手动改后缀即可。工作目录里如果已经存在同名的log、indices、plots新一轮计算会直接覆盖不会提示——这也是为什么建议一个站一个目录别图省事全塞一起。3. 质量控制把哨兵值、负降水和 TminTmax 拦在计算之前RClimDex 主菜单里第二个按钮叫“Load Data and Run QC”很多人把它当成一个可跳过的前置步骤直接奔着指数计算去。这是一个典型的翻车起点。这一步不是“检查一下”而是“改数据”——它会把不合理值替换成NA并生成一份清洗后的中间文件供后续计算使用。跳过它后面的百分位指数会直接把 -99.9 当成真实低温算进第 10 百分位里出来的“冷夜日数”能翻好几倍。3.1 QC 到底改了什么四类值、两种命运质量控制做的事可以拆成两条主线。第一条是无条件替换所有-99.9转成 R 内部的NA所有负的日降水量转成NA所有“日最高温小于日最低温”的记录两个温度值一起转成NA。这三类属于物理上不可能或明确非法的值程序不问你直接动手。第二条是标记而非替换日最高温和日最低温的离群值。判定区间是[均值 - n×标准差, 均值 n×标准差]均值和标准差按“同一日历日的历史序列”算也就是 1 月 1 日只跟历年 1 月 1 日比而不是跟全年比。这个设计是对的——否则季风区一个正常的盛夏高温会被判成离群。关键差异在这里被标记的离群值不会被改成NA它只是被记录进一个清单文件。计算指数时用的仍然是原始值。换句话说QC 阶段程序只替你处理“确定错的”把“可能错的”留给你判断。这个边界划得很清楚但很多人没注意到以为标记了就等于删掉了。3.2 n 取 3 还是 4一个直接影响工作量的参数参数窗口里那个n默认是 3。按正态分布算±3σ 之外大约占 0.27%六十年逐日数据差不多六万条理论上会标出 160 条左右。但实际日温度序列不是正态的尤其日最低温在冬半年有明显偏态真正被标出来的数量经常是理论值的好几倍。手册里明确建议先用 4 起手。这是个很实在的经验值。n3 时你可能拿到几千条标记逐条核对根本不现实n4 能把量级压下来剩下的才值得人工看。n 越大漏检越多n 越小噪音越多没有普适最优我的习惯是先跑 4看标记条数如果还是超过总记录数的 1%说明这个站的温度序列本身波动极大比如高原站、内陆站可以再放到 5如果只有几十条那就逐条看。这里还有个容易忽略的点QC 窗口里有个“Station name or code”输入框很多人纠结要不要填。不用填这个参数是给后续某个环节预留的此时填了也不影响本次 QC 的文件命名——文件名始终取自你选的输入文件名。3.3 四个输出文件谁是黑匣子谁是真正被计算的QC 跑完log/目录下会多出四个 CSV 和四张 PDF。这四张 PDF 是日降水量、日最高温、日最低温、日较差的时间序列图缺测点用红点标出。它们不是装饰一眼扫过去能立刻发现某几年的数据整段缺失或者某个时段温度曲线出现台阶状跳变——后者往往是台站迁址导致的光看数字看不出来。四个 CSV 的分工要记清楚文件内容是否参与指数计算站号tempQC.csv温度类不合理值明细如 TminTmax否仅记录站号prcpQC.csv降水类不合理值明细负值否仅记录站号tepstdQC.csv所有被标记的温度离群值及日期否仅记录站号indcal.csv清洗后的完整日值序列是指数全部来源前三个是黑匣子用来追溯“程序到底动了哪些数”。第四个indcal.csv才是真正的计算底稿里面只有缺测和非法值被换成了NA被标记的离群值原封不动保留。这个设计意味着如果你认为某个离群值确实是错的必须自己去改这个文件程序不会替你决定。3.4 改完数据要不要重跑一个反复被问的问题原始输入文件在整个流程中是只读的QC 不会回写。所以流程有两种走法选错一种就白干第一种直接编辑log/站号indcal.csv把确认错误的离群值改成NA或修正值然后进入指数计算。这条路径快但缺点是可追溯性差——你改了什么、改了多少条没有留痕。第二种回到原始站号.txt修正然后重新走一遍“Load Data and Run QC”让indcal.csv重新生成。这条路径慢但每次改动都能在原始文件里查到。我一般走第二种。理由是indcal.csv是派生文件一旦你手工改过它再想复现同一批结果就得靠记忆。改原始文件虽然多花几分钟但半年后回头看源头在哪一目了然。无论走哪条有一条铁律改了数据就必须重新跑 QC否则指数计算用的还是旧的indcal.csv你会盯着一个没变化的输出怀疑程序坏了。4. 指数计算基准期、纬度符号和四类自定义阈值主菜单第三个按钮“Indices Calculation”才是真正的计算入口但它前面横着一个参数窗口。这个窗口里填错任何一项结果不会报错只会安静地错给你看。尤其是纬度那一栏填错正负号不产生任何提示但会直接改变生长期长度这类跨年指标的季节边界。4.1 参数窗口逐项拆解窗口里一共八项输入分三组。第一组是基准期起止年决定百分位指数的参考时段第二组是纬度决定南北半球的季节窗口第三组是一个降水阈值和四个温度阈值用来生成用户自定义的 Rnn、SUmm、FDmm、TRmm、IDmm 系列指数。参数项含义单位常用取值First year of base period基准期首年四位年1961Last year of base period基准期末年四位年1990Latitude台站纬度南半球填负度39.8 或 -33.9P (Rnn)日降水超过该值计一天mm10、20、25、50Upper Limit of Day HighTX 高于此值的天数SUmm℃35Lower Limit of Day HighTX 低于此值的天数FDmm℃0Upper Limit of Day LowTN 高于此值的天数TRmm℃25Lower Limit of Day LowTN 低于此值的天数IDmm℃-10基准期的选法有讲究。它必须早于或等于数据起始年之后的某个时段而且要足够长WMO 惯例是 30 年才能让第 10、第 90 百分位稳定。如果数据只有 1980 年以后硬填 1961–1990 会让基准期落在数据范围之外百分位指数算出来全是NA。纬度的作用容易被低估。生长期长度这个指数的定义里北半球按 1 月 1 日到 12 月 31 日统计南半球按 7 月 1 日到次年 6 月 30 日统计。纬度填成正数南半球的站就会被按北半球的窗口切分跨年的那一段直接被劈成两半得到一条锯齿状毫无意义的序列。4.2 27 个核心指数不是一锅端参数设好之后会弹出一个指标勾选列表默认全选。全选不是不行但代价是时间而且很多指数跟你的研究问题无关。按数据依赖关系这 27 个指数大致分三档第一档只依赖当日或当年数据比如霜冻日数 FD0TN0℃ 的年天数、夏日日数 SU25TX25℃、冰日日数 ID0TX0℃、热带夜 TR20TN20℃。这一档算得最快也不受基准期影响。第二档依赖固定阈值或年总量比如年总降水量 PRCPTOT、日降水强度 SDII、连续干旱日数 CDD、连续湿润日数 CWD、大雨日数 R10mm 与 R20mm、湿日降水百分位 R95p 与 R99p。R95p 这类指数依赖基准期内的湿日降水分布基准期一变结果全变。第三档依赖百分位和滑动窗口包括冷夜日数 TN10p、暖夜日数 TN90p、冷昼日数 TX10p、暖昼日数 TX90p以及基于这些的冷持续指数 CSDI 和暖持续指数 WSDI。这一档最慢也是 bootstrap 修正的主要受益者。如果只是想快速看一眼数据质量先只勾第一档几分钟出结果。确认流程通了再回到主菜单重新进“Indices Calculation”补算第三档——注意不用重跑 QCQC 结果还在直接算就行。4.3 百分位指数为什么非要 bootstrap这是 RClimDex 相对表格工具的核心优势值得说清楚。计算 TN10p 需要先确定“第 10 百分位”这个门槛常规做法是用整个基准期的所有日最低温算出一个固定值然后逐年统计低于它的天数比例。问题出在 inhomogeneity如果基准期是 1961–1990那么 1961 年的数据既参与了门槛的计算又被拿来跟门槛比较而 2010 年的数据只参与比较不参与门槛计算。这导致序列在基准期末尾出现一个不真实的跳变趋势被系统性高估。解决办法是“留一年法”的 bootstrap算某一年的指数时把该年从基准期中剔除用剩下 29 年重算门槛再统计这一年。这样每一年都是“外来者”序列才齐。这套重采样在表格环境里几乎无法实现但在 R 里只是几行循环。代价是计算时间——第三档指数比第一档慢一到两个数量级六十年数据可能要跑几分钟到十几分钟。慢的时候别以为卡死了程序在算等就行。手册里明确提过“较慢的机器上这一步可能需要几分钟”。4.4 结果怎么读CSV 里的数、JPEG 里的线算完之后每个被勾选的指数在indices/下生成一个站号_指数名.cvs一列年份、一列指数值逗号分隔。plots/下同名 JPEG 是该序列的年度折线图图上叠了两条拟合线实线是普通最小二乘线性趋势虚线是局部加权回归LOWESS平滑。线性拟合的统计量直接印在图上方。这里要提醒一点图上那条线性趋势的显著性是程序按最简单的最小二乘给的没有考虑序列自相关。年尺度的气候指数序列经常有一阶自相关直接读 p 值会偏乐观。要发文章的话趋势显著性得自己另做比如用去趋势预白化之后再拟合。程序图上的数字只适合快速判断方向不适合直接引用。另外indices目录里的值是“每年一个”但像 TXx月最高温的年内最大值这类指数程序实际是先算月值再取年极值。如果你的研究需要月尺度序列得自己从indcal.csv重算输出文件只给到年。5. 避坑清单五类翻车现场的现象、成因与处置以下五条都是这套流程上真实会卡住人的地方按“现象→原因→解决”整理。第一条程序在算百分位指数时突然停住R 控制台没有任何有效报错。现象是勾完指数点确定进度走了一段就没了窗口还在但不再响应。原因是手册第 4 节明确写过的已知缺陷当数据的起始年份与基准期起始年份完全相同时百分位指数的计算需要基准期边界之外的数据程序取不到就中断了。解决办法很土但有效——在输入文件第一行前面补一条记录日期取基准期起始年的前一天六个字段里后三个全填-99.9。比如基准期是 1961–1990数据也从 1961 开始就在文件头加一行1960 12 31 -99.9 -99.9 -99.9。第二条选完文件后弹窗报错提示格式不正确。现象是文件明明存在、也能用记事本打开界面就是不认。原因通常有三种文件带表头行、字段之间用了逗号、列数不是 6。RClimDex 的读取逻辑非常死板不猜分隔符也不跳表头。解决方式是用第 2 章那两行awk先扫一遍把表头删掉、把逗号换成空格、把多出来的列比如站点号、质量标志位砍掉。第三条南半球的站生长期长度序列全是 0 或者锯齿状。现象是指数能算出来但数值明显不对跨年那一段完全断裂。原因是纬度填了正数。程序按纬度正负判断半球进而决定 GSL 的统计窗口是 1 月–12 月还是 7 月–次年 6 月。解决方式是回到参数窗口把纬度改成负值重算该指数。注意这一步不需要重跑 QCQC 结果与纬度无关。第四条改了indcal.csv里的错误值重算指数结果一点没变。现象是明明把某个离谱的低温改成了NA输出序列纹丝不动。原因是指数计算读的是indcal.csv但很多人改的是tempQC.csv或者tepstdQC.csv——后两个只是记录文件不参与计算。还有一种情况是改对了文件但没有重新点“Indices Calculation”程序用的还是上一轮内存里的数据。解决方式确认改的是log/站号indcal.csv改完退出主菜单重新进计算入口或者干脆回原始文件改完重跑 QC。第五条离群值标记文件里几千条记录误以为数据烂了把整个站弃用。现象是tepstdQC.csv打开一看密密麻麻直觉判断这份数据不可用。原因是 n 取了默认值 3对偏态明显的日最低温序列来说这个门槛太紧标出的多数是正常的寒冬极值。解决方式是把 n 调到 4 甚至 5 重跑 QC观察标记数量是否落到合理区间同时抽几条标记记录回原始数据里对照同期的邻站数据确认是真错还是正常极值。只要不是成片连续标记通常不影响指数计算——毕竟这些值根本没被替换。6. 批量跑站与结果自检把 27 个指数收进总表之前先过三道校验这套工具是交互式的一站一窗口几十个站手工点下来既慢又容易漏。常见做法是分两层处理前置的数据体检和参数生成用脚本批量化后置的结果汇总与校验也用脚本批量化中间那层真正需要点鼠标的 QC 和计算保持手工但把顺序固定下来。第一道校验在算之前把每个站的六列格式、日期连续性、缺测比例一次性扫出来生成一张清单。缺测率超过 20% 的站直接标记出来这类站的百分位指数不可信。files - list.files(raw, pattern \\.txt$, full.names TRUE) chk - do.call(rbind, lapply(files, function(f) { x - read.table(f, header FALSE, na.strings -99.9, col.names c(year,month,day,prcp,tmax,tmin)) data.frame(station basename(f), n nrow(x), miss_prcp mean(is.na(x$prcp)), miss_tmax mean(is.na(x$tmax)), start x$year[1], end x$year[nrow(x)]) })) subset(chk, miss_prcp 0.2 | miss_tmax 0.2) # 需要人工确认的站第二道校验在算之后拿输出指数做内部一致性检查。R10mm的天数不可能超过当年的降水日数FD0与SU25之和不可能超过 365TXx必须大于等于当年TX的均值。这些约束任何一条被打破说明输入或参数出了问题不用往下看趋势。第三道校验是交叉复算从indcal.csv独立算一遍PRCPTOT跟indices/站号_PRCPTOT.cvs对差值应当严格为 0。qc - read.csv(log/21946indcal.csv) out - read.csv(indices/21946_PRCPTOT.cvs) my - tapply(qc$prcp, qc$year, sum, na.rm TRUE) max(abs(my - out[, 2])) # 应为 0非 0 说明中间有人动过文件这套检查用不了十分钟但能拦住九成的“结果看着怪”的问题。早年我做一批站的指数汇总图省事跳过了第三道直接把二十几个 CSV 拼在一起等到写分析时才发现其中一个站的基准期填成了数据起始年百分位指数整列缺失前面所有图白画。从那以后我每次批量跑完都强制走一遍这三道校验宁可慢半小时也不返工一周。希望帮到你。本文还有配套的精品资源点击获取