RClimDex实战：极端气候指数计算与数据质控全流程

简介：RClimDex 是加拿大气象部门研发的极端气候指数计算软件，已被 CCl/CLIVAR ETCCDMI 国际工作坊广泛使用，面向气候科学家、研究人员及气候变化专业学生。软件包内为 1 份 doc 格式的完整用户手册，约 600KB，系统讲解 R 环境的安装运行、RClimDex 数据质量控制、缺失值与异常值处理，以及高温日数、低温日数、极端降雨量等 27 个核心指数的计算方法。手册还包含输入数据格式、指数定义、阈值与基准期温度指数计算方式、R 在 Windows 环境下的常见故障解答，附录列有完整的气候指数清单，便于读者按需查阅和快速排错。目前已有 3836 人学习下载。对需要借助 R 平台开展极端气候事件监测、理解气候变率并基于客观数据作出气候变化分析的读者，这是一份实用且精炼的入门参考工具。

1. 极端气候指数计算：为什么我建议你直接选 RClimDex

做气候变化检测的人，十有八九绕不开那 27 个 ETCCDMI 核心指数。早年大家用 Excel 版的 ClimDex 算，后来发现基于百分位的温度指数序列存在不均匀性，需要引入 bootstrap 过程才能修正，这在 Excel 环境里根本跑不动。RClimDex 就是为了解决这个问题而生的——它把数据质量控制、27 个核心指数计算和结果可视化全部封装在一个 R 脚本里，装上 R 环境就能用。适合谁？做气候监测的业务人员、研究极端降水温度的在校学生、以及要给气象数据做批量质控和指数提取的工程师。它对输入数据格式有严格要求，但只要按规范准备好站点日值文件，从质控到出图基本一条龙。这篇笔记我会把安装、数据格式、质控参数、指数计算配置和各类踩坑记录完整过一遍。

2. 安装 R 与加载 RClimDex：环境配置和首个报错排查

2.1 先装 R：版本选择和 TclTk 依赖

RClimDex 1.0 是在 R 1.8.4 下开发的，官方手册说 1.8.4 或更高版本都能跑。我自己的经验是，别用太老的版本，也别盲目追最新——R 4.x 下跑 RClimDex 1.0 的旧脚本偶尔会遇到一些函数兼容性警告，但基本不影响核心计算。如果你用的是 Windows，直接从 CRAN 镜像下载安装程序，默认安装就会带上 TclTk 支持。TclTk 是 RClimDex 图形界面的底层依赖，R 1.9.0 之后默认包含，但如果你因为某些原因装了精简版 R，记得确认 TclTk 是否可用。

安装完 R 之后，在 Windows 桌面双击 R 图标进入 RGui。Linux 或 Unix 环境下直接敲 R 命令进入控制台。退出统一用 q()。这里有一个容易忽略的点：某些 Windows 机器上 R 会报错提示找不到 HOME 环境变量，这属于 R for Windows FAQ 里的老问题，在系统环境变量里新增一个 HOME 指向任意工作目录就能解决。

提示：安装路径不要带中文和空格。R 本身对路径兼容性尚可，但 RClimDex 读取数据文件和输出结果时会对路径字符串做拼接处理，纯英文路径能帮你省掉一堆莫名其妙的问题。

2.2 加载 RClimDex 的两种方式

RClimDex 不需要 install.packages，因为早期版本就是一个 rclimdex.r 源文件。加载方式取决于你手头文件在哪。

如果你已经把 rclimdex.r 下载到本地 D:\climdex\ 目录，在 R 控制台运行：

source("D:/climdex/rclimdex.r")

如果你联网且想拉取当时的最新版，可以直接从远程地址加载：

source("http://cccma.seos.uvic.ca/ETCCDMI/RClimDex/rclimdex.r")

加载成功的标志是 R 控制台弹出 RClimDex 主菜单，菜单项包括 Load Data and Run QC、Indices Calculation 和 Exit。如果 source 之后控制台报错，最常见的原因是文件路径写错或网络源不可达。我一般建议把 rclimdex.r 下载到本地再用，远程 source 在后面的调试过程中会因为你反复修改脚本而变得非常别扭。

提示：这里的 http 地址是 2004 年手册时代的产物。时至今日该地址未必仍然有效。如果你加载失败，先检查网络连通性，再去找随资源包附带的本地 rclimdex.r 文件。

2.3 环境变量和运行时异常排查

RClimDex 的图形界面依赖 tcltk 包。如果你加载时看到类似 "there is no package called 'tcltk'" 的报错，说明你的 R 安装不完整，需要重新安装完整版 R 或单独安装 tcltk 包。另一个高频问题是 Windows 下 RGui 菜单里找不到 "Source R code" 选项——这是因为你用的是 R Console 而不是 RGui，右键菜单或 File 菜单里才会出现这个入口。

我遇到过一次比较特殊的场景：在某台服务器上通过 SSH 跑 R，没有图形环境，RClimDex 主菜单弹不出来。这种情况下就别硬试了，RClimDex 1.0 的设计强依赖交互式图形界面，无头服务器环境建议改用其他支持命令行批处理的工具。

3. 数据格式与质量控制：这一步决定你后面能不能算出靠谱的指数

3.1 输入文件格式：严格按照五列结构准备

RClimDex 对输入数据格式要求非常严格。官方手册附录 B 给出了明确规范，我这里拆开讲。文件必须是纯文本 ASCII 格式，文件名形如 stationname.txt。每一行代表一天，共五列，列之间用空格或制表符分隔，按顺序是：年、月、日、日最高气温 TX、日最低气温 TN、日降水量 PRCP。

一个标准的示例数据长这样：

1961 01 01 12.5 3.2 0.0 1961 01 02 14.1 5.0 2.5 1961 01 03 -99.9 -99.9 -99.9

第三行是缺失日期的记录方式——用 -99.9 表示该日数据缺失。这里有一个容易踩坑的点：RClimDex 判断缺失值只看 -99.9，其他任何缺失编码（比如 -999、NA、空值）都会被视为真实观测值，从而污染后续质控和指数计算。所以准备数据时，务必把所有缺失值统一替换为 -99.9，并且确保文件中没有空行、没有逗号分隔、没有表头。

提示：如果你从数据库导出数据，最常见的格式问题就是时间字段带斜杠或横杠（比如 1961-01-01），或者降水量字段带小数点以外的特殊符号。建议在导入 RClimDex 之前先用脚本清洗一遍。

3.2 加载数据进入质控流程

在 RClimDex 主菜单选择 "Load Data and Run QC"，会弹窗让你选择数据文件。选好站点文件后，RClimDex 会自动完成两件事：把所有 -99.9 转换成 R 能识别的 NA；把所有不合理值替换为 NA。不合理值包括两种情况——日降水量小于零、日最高气温低于日最低气温。

# 假设你有一个原始数据文件 raw_data.txt，这是用 awk 做格式预处理的示例 # 把 -999 和空值统一替换为 -99.9 awk '{if ($4=="") $4="-99.9"; if ($5=="") $5="-99.9"; if ($6=="") $6="-99.9"; print}' raw_data.txt > clean_data.txt

这个命令的逻辑很简单：逐行读取原始文件，检查第 4、5、6 列是否为空，为空则填 -99.9，最后输出到新文件。你也可以顺手把非数值字符过滤掉，但我更推荐在导出环节就处理干净，awk 只做缺失值统一。这里要注意 awk 不会处理文件中可能存在的表头，如果原始数据有表头行，需要额外加一行判断，或者用 sed 直接删掉第一行：

sed '1d' raw_data.txt | awk '{if ($4=="") $4="-99.9"; print}' > clean_data.txt

3.3 质控参数设置：阈值 n 到底该填 3 还是 4

数据加载成功后，RClimDex 会弹出 "Set Parameters for Data QC" 窗口，需要你输入一个关键参数 n。这个 n 用于判断温度数据中的离群值——以某个日期多年平均值为中心，均值加减 n 倍标准差构成一个区间，日最高温或日最低温落在区间之外就被标记为潜在离群值。默认值是 3，但手册明确建议用户可以从 4 开始试。

为什么建议从 4 开始？因为 n=3 会把大量正常气候波动标记为离群值。我自己处理北方站点数据时，用 n=3 跑一遍经常标记出上百个离群点，人工核查根本看不过来。n=4 相对宽松，标记数量会降到可操作范围。如果你处理的是高寒地区或热带站点，气温日较差大的话，n=4 仍然可能标记过多，这时候可以再加大到 5。

质控参数窗口里还有 "Station name or code" 输入框，这里可以留空，因为当前版本这个字段不会用到。填了也不影响计算。点击 OK 之后程序开始跑质控，数据量大的时候在老旧电脑上可能耗时几分钟，这是正常现象，别误以为死机。

3.4 质控输出文件：四个 CSV 一个都不能少

质控完成后，RClimDex 会在当前工作目录下创建一个 log 子目录，并生成四个文件：

其中 indcal.csv 是后续计算的基础。注意一个细节：不合理值和缺失值在这里已经被替换为 NA，但被标记为离群值的温度数据并不会被自动修改或替换，只是单独列在 tepstdQC.csv 里供你人工判断。这样做是合理的，因为离群值不一定是错误值——某个站点某天出现极端高温，可能是真实的气候事件，不该被质控流程粗暴抹掉。

log 目录下还会生成四个 PDF 文件，分别是日降水量、日最高温、日最低温和日较差的时间序列图，缺失值用红点标注。你可以在这些图上快速判断数据整体质量。如果发现某个时间段红点连成一片，那大概率是观测仪器故障或数据录入问题，应该回头修正原始数据文件，而不是在 indcal.csv 里手动改。

3.5 质控后的数据修正逻辑

这里有一个很多人搞错的逻辑：RClimDex 的指数计算是基于质控后的 indcal.csv，不会改动你的原始输入文件。如果你发现原始数据里有错误的非极端值（比如某天温度录错了），你当然可以手动编辑 indcal.csv 来修正，但更规范的做法是修正原始 txt 文件后，重新执行 Load Data and Run QC。因为 indcal.csv 里的离群值标记和 PDF 图都是基于原始数据自动生成的，只改 indcal.csv 不重新质控，会导致最终指数结果和可视化输出不对应。

我自己的习惯是：第一步跑 n=4 质控，导出 tepstdQC.csv，用 Excel 或文本编辑器逐个核查标记日期；第二步把确认无误的离群值保留，确认有误的回到原始文件修正；第三步重新加载修正后的文件再跑一次质控。这个过程看似麻烦，但在数据量大的时候能帮你省掉后面指数异常排查的大量时间。

4. 27 个核心指数计算：基准期、纬度、阈值参数配置

4.1 计算前必须理解的几个参数

进入 Indices Calculation 后，首先要配置一组参数。这组参数直接决定哪些指数能算、算出来是什么含义。核心参数有四个：

基准期的选择非常重要。RClimDex 用基准期内每一天的气温分布来计算第 10 和第 90 百分位，然后统计整个时间序列中低于第 10 百分位或高于第 90 百分位的天数。基准期内数据缺失太多的话，百分位计算会不稳定。我一般要求基准期内每个日历年至少要有 80% 的有效天数，否则会考虑换基准期。

站点纬度的作用体现在 Growing Season Length（GSL）指数上。北半球 GSL 计算从 1 月 1 日到 12 月 31 日，南半球则是从 7 月 1 日到次年 6 月 30 日。这个参数填错不会报错，但算出来的 GSL 序列会完全不符合气候意义。我在处理南美站点数据时踩过这个坑，当时忘了把纬度填成负值，结果 GSL 在 7 月出现大段缺失，后来检查参数才发现问题。

4.2 用户自定义阈值：扩展 27 个核心指数之外

除了 ETCCDMI 核心指数，RClimDex 还支持自定义温度的 SUmm、FDmm、TRmm、IDmm 四类指数。这里的 mm 就是你在参数窗口里填的四个温度阈值：

• User defined Upper Limit of Day High：日最高温度大于该阈值的天数
• User defined Lower Limit of Day High：日最高温度小于该阈值的天数
• User defined Upper Limit of Day Low：日最低温度大于该阈值的天数
• User defined Lower Limit of Day Low：日最低温度小于该阈值的天数

我一般会按业务需求填四组常见阈值。比如对农业区，我会把 Day High 上限填 35 摄氏度，统计高温热害天数；把 Day Low 下限填 0 摄氏度，统计霜冻天数。对城市热岛研究，Day High 上限填 30，Day Low 上限填 25（对应热带夜）。注意这些阈值和核心指数里固定的 SU25、FD0、TR20、ID0 是独立的，互不影响，你可以同时算默认的和自定义的。

# 如果你不想走图形界面，可以尝试直接调用内部函数（需要自行修改 rclimdex.r） # 伪代码示意，实际函数名以脚本内定义为准 # indcal <- read.table("log/21946indcal.csv", header=TRUE, sep=",") # compute.indices(indcal, base.start=1961, base.end=1990, lat=31.5)

这段代码体现的是核心逻辑：读取质控后数据，传入基准期起止年和纬度，然后计算全部指数。实际场景中你大概率还是用图形界面操作，但理解了底层函数调用关系，后续调试会容易得多。

4.3 指数选择：全选还是按需勾选

参数设置好之后点击 OK，程序会进入指数选择窗口，默认全部 27 个指数都是勾选状态。如果你是第一次跑某个站点的数据，我强烈建议先全选跑一遍，让程序把整套指数序列都生成出来，你再看哪些指标在目标区域的信号确实有意义。全选计算耗时更长，但对理解站点气候特征帮助很大。

如果你已经确定研究目标，比如只关注极端降水相关的 R95p、R99p、Rx1day、Rx5day，那可以只勾选这几个，计算会快很多。这里要注意一个取舍：RClimDex 的指数文件是按指数分别保存的，勾选得越少生成的 CSV 文件越少，后续整理汇总的工作量也越小。

4.4 计算结果：CSV 序列与 JPEG 趋势图

指数计算完成后，RClimDex 会在当前目录下创建 indices 子目录，每个指数生成一个独立的 CSV 文件，文件名规则是"站点名_指数名.csv"，列之间用逗号分隔。同时创建一个 plots 子目录，存放每个指数的年度序列图，图上有两条拟合线：实线是线性最小二乘拟合，虚线是局部加权回归拟合，图上还标注了线性趋势的统计量。

# 以 TXx 为例的指数输出文件内容示意 YEAR,TXx 1961,33.5 1962,34.1 1963,32.8

CSV 文件没有表头说明，第一列是年份，第二列是指数值。你自己做后续分析时，可以直接用 read.csv 读入 R 或 pandas 读取，列名自己指定即可。不过注意有些早期版本输出的文件可能在末尾多一个逗号，用 pandas 读取时稍不注意就会多一列全 NaN，我在处理多个站点批量读取时被骗过一回。

5. 避坑与常见问题：RClimDex 典型报错和处理记录

5.1 程序启动即停止：基准期起点与数据起点重合

这是一个官方手册明确记载的已知 bug。现象是：当输入数据的起始年份和基准期起始年份相同，程序会在计算百分位温度指数时突然停止运行，控制台没有任何明确的错误提示，只是卡住不动或直接退出。原因是计算百分位指数时需要用到基准期前一年的数据作为边界条件，而数据起点恰好是基准期起点时，边界外的数据不存在，计算过程就会崩溃。

解决方法是往输入数据文件第一行前面加一条虚拟记录。假设基准期是 1961-1990，数据也从 1961 年开始，那就补一行 1960 年 12 月 31 日的记录，气温和降水全部填 -99.9：

1960 12 31 -99.9 -99.9 -99.9 1961 01 01 12.5 3.2 0.0

这条记录不会影响实际指数计算，因为它是缺失值，但提供了边界计算所需的占位。我处理过的站点里大概有三分之一会撞上这个坑，所以现在准备数据时会主动检查数据起始年和基准期起始年是否相同，相同则提前补行。

5.2 质控阶段大量标记离群值：n 值过小

现象是：跑完质控后 tepstdQC.csv 里标记了上千个离群点，人工核查根本看不过来，而且大部分离群点在气候意义上其实是正常值。原因就是你设置的 n 值太小。n=3 意味着均值加减 3 倍标准差，对于日气温这种季节性很强的变量，标准差在过渡季节会偏大，3 倍标准差仍然会把很多正常波动挡在外面。

解决方法是把 n 从 3 改到 4 或 5。但这里的 n 值调大不一定会减少标记数量——因为质控用的是逐日气候平均值和标准差，如果站点的高原特征明显、气温日较差本来就大，n=5 也可能标记出不少点。更合理的做法是先跑 n=5 看一次标记数量，再逐步调小到 n=4 对比，找到一个标记数量适中且不遗漏真正异常的平衡点。

注意：n 只影响离群值的标记，不影响不合理值的替换。降水量小于零、最高温低于最低温这种情况无论 n 是多少都会被自动替换为 NA。

5.3 输入文件报错：格式不符合规范

现象：点击 Load Data and Run QC 后，R 控制台报错，弹窗提示数据加载失败。最常见的原因是输入文件格式不正确——具体来说可能是列数不对、分隔符不统一、包含表头或空行。RClimDex 的读取逻辑是按固定列数解析的，多一列少一列都会直接报错。

解决方法是严格对照示例数据检查格式。我建议准备数据后用 head 命令先看一眼前几行：

head -5 21946.txt

确保每一行都是五列，并且列顺序是"年 月 日 TX TN PRCP"。如果你不确定原始数据里有没有隐藏字符，可以用 cat -A 查看行尾和特殊符号：

cat -A 21946.txt | head -5

如果这一行输出里有 ^M 之类的符号，说明文件是 Windows 换行符，部分 Unix 环境下 R 读取可能会出问题，用 dos2unix 转换一下即可。

5.4 指数计算结果全为 NA：基准期内数据缺失过多

现象是：计算完成后，某些百分位指数（TN10p、TX10p、TN90p、TX90p）的序列几乎全是 NA 或大片缺失。原因是基准期内每天计算百分位时需要足够的历史样本，如果基准期内的数据大量缺失，某个日期的百分位无法被估计，整个计算就返回 NA。

解决方法是检查基准期内的数据完整性。如果基准期内某个月整月缺失，那这个月的每日百分位都会是 NA。你可以考虑缩短基准期，比如从 1961-1990 改为 1971-2000，前提是后半段数据质量更好。另一个思路是接受这些 NA，在后续趋势分析时用线性插值或直接剔除缺失年份，但这样做会降低结果的可靠性。

5.5 输出文件打不开：目录或文件名中的特殊字符

现象是：RClimDex 报告计算完成，但生成的 CSV 文件在 Excel 里打开乱码，或者文件路径里找不到 indices 目录。原因有两种可能：一是站点文件名包含特殊字符，导致生成的文件名在 Windows 下非法；二是输出目录中包含中文或空格，R 的路径处理在某些版本下会出问题。

解决方法是规范站点文件名，只用字母、数字和下划线命名。目录路径避免中文，避免空格，统一用英文下划线。我一般把所有站点数据放在 D:\rclimdex_data\raw 下，输出自然也在 D:\rclimdex_data 下生成 log、indices、plots 三个子目录，全程英文路径，从未出过问题。

6. 进阶验证：用 R 脚本批量复算一个指数并做趋势显著性检验

算完指数拿到 CSV 序列只是第一步，真正要回答"这个站点的极端气候到底有没有变"这个问题，还得对指数序列做趋势估计和显著性检验。RClimDex 生成的 JPEG 图里已经画了线性拟合线并显示统计量，但那个统计量信息有限，通常只有斜率估计，没有给出置信区间和 p 值。我习惯把某个重点指数的 CSV 读进 R，用 lm 函数重新做一次线性回归，并提取完整的统计信息。

# 读取某个指数的年度序列 txx <- read.csv("indices/21946_TXX.csv", header = TRUE) # 查看数据结构，确认列名 str(txx) # 线性趋势拟合 fit <- lm(TXX ~ YEAR, data = txx) # 输出回归系数、标准误、t 值和 p 值 summary(fit) # 提取斜率并计算每十年的变化量 slope <- coef(fit)[2] decadal_change <- slope * 10 cat("每十年变化量：", decadal_change, "\n")

这段代码的逻辑分四步：先用 read.csv 读入指数序列；然后用 lm 函数以年份为自变量、指数值为因变量做最小二乘回归；再通过 summary 查看完整回归结果，包括斜率的标准误和显著性；最后计算每十年的变化量，方便在报告里统一表述。注意这里的 read.csv 默认第一行是列名——如果你的输出文件第一列叫 YEAR 第二列叫 TXx，那就直接能用；如果文件没有表头，需要加 header = FALSE 并手动指定列名。

参数说明：decadal_change 的单位取决于原指数的单位。TXx 单位是摄氏度，所以每十年变化量的单位是摄氏度/十年。如果是 Rx1day（日最大降水量），单位就是毫米/十年。对大多数 ETCCDMI 指数，用每十年变化量来表述趋势比用年斜率更直观，也更容易和其他文献对比。

我自己处理某个站点数据时发现 TXx 的线性趋势 p 值是 0.08，RClimDex 的图上线性拟合线看起来有明显上升，但严格来说没通过 0.05 显著性检验。这种情况在业务报告里就需要特别谨慎表述，不能说"显著上升"，只能说"呈上升趋势但未达显著水平"。从那以后我每次拿到指数序列，都会先跑一遍 lm 提取 p 值再下结论，而不是看着图上的斜率就写报告。希望这篇笔记能帮你少走几段弯路。

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