如何使用 R 语言计算各省市区县相对湿度面板数据

今天给大家分享如何使用 R 语言从 NASA MERRA-2 卫星再分析数据中计算各省市区县的相对湿度日度面板数据的方法，本课程包含如下内容：

1. 从 NASA GES DISC 下载 MERRA-2 NetCDF（.nc）格式数据；
2. 使用 terra 包读取并处理栅格数据；
3. 根据气象公式计算相对湿度（RH）；
4. 使用 terra::extract() 将栅格汇总为省市区县面板数据；
5. 合并多天结果并导出为 Excel 格式。

数据来源

本课程使用的原始数据来自 NASA MERRA-2（Modern-Era Retrospective Analysis for Research and Applications, Version 2）逐小时近地表气象数据集：

数据集：M2I1NXLFO 5.12.4
数据类型：NetCDF（.nc）格式，逐小时，覆盖全球

该数据集包含以下变量：

变量名	含义	单位
HLML	地表层高度（Surface Layer Height）	m
PS	地表气压（Surface Pressure）	Pa
QLML	近地表比湿（2m Specific Humidity）	kg/kg
SPEEDLML	近地表风速（Surface Wind Speed）	m/s
TLML	近地表气温（2m Temperature）	K

其中**比湿（specific humidity）**是指在一团湿空气中，水汽的质量与该团空气总质量的比值。如果湿空气与外界无质量交换，并且没有相变发生，那么比湿会保持不变。比湿通常以 g/g 或 g/kg 为单位，通常大气中的比湿都小于 40 g/kg，是记录大气水汽状况的重要指标。

相对湿度的计算原理

第一步：下载 MERRA-2 数据

这部分内容可以参考：

名师讲堂｜逆温数据是如何处理的？使用 R 语言完成整个过程：https://rstata.duanshu.com/#/brief/course/55b47a7a83124138b918e7475f4141d9

注册 NASA 账号并获取下载链接

首先前往 NASA GES DISC 注册账号：https://disc.gsfc.nasa.gov

然后打开数据集页面，点击右下角 Subset/Get Data，按以下设置筛选：

1. Download Method：选择 Get File Subsets using OPeNDAP；
2. Refine Date Range：选择所需时间段，例如 2025-01-01 到 2025-12-31；
3. Refine Region：输入中国经纬度范围。

中国的经纬度范围可以通过 R 代码获取：

library(sf)
sf::read_sf("2021行政区划/省.shp") -> prov
st_bbox(prov) %>% round(3) %>% paste0(collapse = ",")

1. 变量全部勾选（HLML、PS、QLML、SPEEDLML、TLML）；
2. 点击 Get Data → Download links list，下载一个 .txt 格式的下载链接文件。

注意：下载链接的有效期仅为 两天，过期后需要重新生成。

注意要选择 Mean of complete time range，否则下载得到的是小时数据。

配置 NASA 登录凭证（Mac/Linux）

在终端依次运行以下命令（将 <uid> 和 <password> 替换为你的 NASA 账号）：

cd ~
touch .netrc
echo "machine urs.earthdata.nasa.gov login <uid> password <password>" >> .netrc
chmod 0600 .netrc
touch .urs_cookies

批量下载 nc 文件

将下载链接文件放到工作目录的 data 文件夹下，然后在终端运行：

cd "你的工作目录/data/"
cat subset_M2I1NXLFO_5.12.4_xxxxx.txt | tr -d '\r' | xargs -n 1 curl --globoff -LJO -n -c ~/.urs_cookies -b ~/.urs_cookies

由于数据量较大，可以使用 R 将链接拆分成多段，同时开多个终端窗口并行下载：

library(tidyverse)

# 读取下载链接文件
read_lines("subset_M2I1NXLFO_5.12.4_20260228_092259_.txt") %>%
  as_tibble() %>%
  rename(X1 = value) -> df

# 去除已经下载好的文件，避免重复下载
df %>%
  mutate(filename = str_match(X1, "\\d{8}")[,1]) %>%
  anti_join(
    fs::dir_ls("data") %>%
      as.character() %>%
      as_tibble() %>%
      mutate(filename = str_match(value, "\\d{8}")[,1])
  ) %>%
  select(X1) -> df

# 例如按 74 条链接一组拆分
lapply(seq(1, nrow(df), by = 74), function(x){
  df %>%
    slice(x:(x+73)) %>%
    write_csv(paste0("data/downlist", x, ".txt"), col_names = F, quote = "none")
}) -> res

# 输出每组对应的 curl 命令
for (x in seq(1, nrow(df), by = 74)) {
  message("cd 你的工作目录/data/")
  paste0("cat downlist", x, ".txt | tr -d '\\r' | xargs -n 1 curl --globoff -LJO -n -c ~/.urs_cookies -b ~/.urs_cookies") %>%
    message()
  message(" ")
}

下载完成后，data 文件夹中会包含多个类似 MERRA2_400.inst1_2d_lfo_Nx.20250101.SUB.nc 的文件：

fs::dir_ls("data")

第二步：读取 nc 文件并转换为日均 tif

nc 数据可以直接使用 terra 包读取：

library(tidyverse)  # 加载数据科学包集合
library(sf)         # 处理矢量空间数据
library(terra)      # 处理栅格和 NetCDF 数据

# 获取所有 nc 文件路径
fs::dir_ls("data") -> fls

# 读取单个文件查看结构
rast(fls[1]) -> rst
rst

可以看到，每个文件有 5 层（5 个变量）。

第三步：计算相对湿度

本课程提供了一个封装好的函数 calculate_merra_rh()，位于 calculate_merra_rh.R 文件中：

source("calculate_merra_rh.R")

函数内部实现如下：

#' 从MERRA-2栅格数据计算相对湿度
#'
#' @param merra_stack 包含MERRA-2变量的SpatRaster堆栈，必须包含以下层：
#'                    "PS"(地表气压，Pa)、"QLML"(2m比湿，kg/kg)、"TLML"(2m温度，K)
#' @return 包含相对湿度(%)的SpatRaster对象
#'
#' @examples
#' # 假设你已经有一个包含MERRA-2数据的 SpatRaster 堆栈
#' # merra_data <- c(ps_rast, qlml_rast, tlml_rast) # 组合为堆栈
#' # names(merra_data) <- c("PS", "QLML", "TLML")
#' # rh_rast <- calculate_merra_rh(merra_data)
calculate_merra_rh <- function(merra_stack) {
  # 检查必要的包
  if (!requireNamespace("terra", quietly = TRUE)) {
    stop("请先安装'terra'包: install.packages('terra')")
  }

  # 验证输入数据
  required_layers <- c("PS", "QLML", "TLML")
  if (!all(required_layers %in% names(merra_stack))) {
    missing_layers <- setdiff(required_layers, names(merra_stack))
    stop(paste("缺少必要的栅格层:", paste(missing_layers, collapse = ", ")))
  }

  # 从堆栈中提取各变量
  ps <- merra_stack[["PS"]]    # 地表气压(Pa)
  q <- merra_stack[["QLML"]]   # 2m比湿(kg/kg)
  tk <- merra_stack[["TLML"]]  # 2m温度(K)

  # 计算实际水汽压(e) - 单位hPa
  # e = (q * P) / (0.622 + 0.378 * q)
  # 先将PS从Pa转换为hPa
  e <- (q * (ps / 100)) / (0.622 + 0.378 * q)

  # 计算饱和水汽压(e_s) - 单位hPa
  # 使用改进的Magnus公式(适用于水面)
  # e_s = 6.112 * exp((17.67 * T_C) / (T_C + 243.5))
  tc <- tk - 273.15  # 转换为摄氏度
  e_s <- 6.112 * exp((17.67 * tc) / (tc + 243.5))
  # 计算相对湿度(%)
  rh <- (e / e_s) * 100

  # 设置输出栅格的名称和单位
  names(rh) <- "RH"
  terra::units(rh) <- "%"

  # 限制RH在0-100%范围内(处理可能的数值误差)
  rh <- terra::clamp(rh, 0, 100)

  return(rh)
}

对单个文件计算相对湿度并可视化：

# 读取某天的
rast(fls[1]) -> rst

# 计算相对湿度
rhrst <- calculate_merra_rh(rst)
rhrst

# 绘制栅格图
plot(rhrst)

循环处理所有天的数据，生成相对湿度栅格并保存：

# 循环所有的
dir.create("rhtif")

lapply(fls, function(x){
  # 日期
  str_extract(x, "\\d{8}") -> y
  rast(x) -> rst
  rhrst <- calculate_merra_rh(rst)
  rhrst %>%
    writeRaster(paste0("rhtif/", y, ".tif"), overwrite = T)
}) -> res

第四步：读取行政区划矢量数据

将栅格数据汇总到省市区县，需要先读取行政区划的 shp 文件（本课程附件提供 2021 年行政区划数据）：

library(sf)

# 分别读取省、市、县三级行政区划
read_sf("2021行政区划/省.shp") %>%
  select(-contains("类型")) %>%
  st_transform(4326) -> prov

read_sf("2021行政区划/市.shp") %>%
  select(-contains("类型")) %>%
  st_transform(4326) -> city

read_sf("2021行政区划/县.shp") %>%
  select(-contains("类型")) %>%
  st_transform(4326) -> county

注意：如果需要处理乡镇级数据，可以额外读取乡镇 shp 或 rds 文件，代码结构相同。

第五步：提取各省市区县均值

使用 terra::extract() 对每天的原始 MERRA-2 栅格（含 HLML、PS、QLML、SPEEDLML、TLML）和相对湿度栅格分别提取各行政区划的区域均值，结果保存在 res1（原始变量）和 res2（相对湿度）两个文件夹中：

此处代码需下载讲义材料查看～

提取结果示例（省级）：

read_rds("res1/20250110_prov.rds")

第六步：合并所有天的结果

将 res1 和 res2 中每天的结果分别合并，再 left_join 到一起，得到包含所有变量的完整面板数据：

此处代码需下载讲义材料查看～

查看省级面板结果：

dfprovall

相对湿度的频率分布：

hist(dfcityall$RH,
     main = "城市相对湿度分布",
     xlab = "相对湿度（%）",
     col = "steelblue", border = "white")

第七步：保存结果

# 保存为 rds 文件
dfcityall   %>% write_rds("dfcityall.rds")
dfcountyall %>% write_rds("dfcountyall.rds")
dfprovall   %>% write_rds("dfprovall.rds")

dfprovall %>%
  arrange(省, 省代码, date) %>%
  writexl::write_xlsx("各省份相对湿度日度数据.xlsx")

# 导出城市数据（按年份拆分，避免单个文件过大）
dfcityall %>%
  mutate(year = year(date)) -> dfcityall

dir.create("各城市相对湿度日度数据")
lapply(unique(dfcityall$year), function(x){
  dfcityall %>%
    filter(year == x) %>%
    writexl::write_xlsx(paste0("各城市相对湿度日度数据/", x, ".xlsx"))
})

# 导出区县数据（按年份拆分，每年再拆成两份，否则可能装不上）
dfcountyall %>%
  mutate(year = year(date)) -> dfcountyall

dir.create("各区县相对湿度日度数据")
lapply(unique(dfcountyall$year), function(x){
  dfcountyall %>% filter(year == x) %>%
    slice(1:500000) %>%
    writexl::write_xlsx(paste0("各区县相对湿度日度数据/", x, "_part1.xlsx"))
  dfcountyall %>% filter(year == x) %>%
    slice(500001:nrow(.)) %>%
    writexl::write_xlsx(paste0("各区县相对湿度日度数据/", x, "_part2.xlsx"))
})

数据结构说明

最终面板数据包含如下变量：

变量名	含义	单位
省/市/县	行政区划名称	—
省代码/市代码/县代码	行政区划代码	—
date	日期	—
HLML	地表层高度	m
PS	地表气压	Pa
QLML	近地表比湿	kg/kg
SPEEDLML	近地表风速	m/s
TLML	近地表气温	K
RH	相对湿度	%

年度和月度数据可以使用日度数据分组汇总得到，大家可以独立完成 group_by + summarise 的汇总操作。

注意事项

1. 数据量较大：全国区县级的日度数据每天约 2800 条记录，一年约 100 万条，读取和处理时内存消耗较大，建议分年处理或使用 parallel 包多线程加速；
2. 依赖包版本：本课程使用的是 terra 包 1.7+ 版本，旧版本可能在读取某些 nc 文件时出现坐标系识别问题，如遇此类情况可手动设置：

ext(rst) <- c(73.4375, 135.3125, 3.75, 53.75)
crs(rst) <- "+proj=longlat +datum=WGS84 +no_defs"

1. 比湿单位：MERRA-2 数据中 QLML 的单位为 kg/kg（与 g/g 等价），在相对湿度公式中直接使用，无需换算；
2. 相对湿度合理性校验：计算结果已通过 terra::clamp(rh, 0, 100) 限定在 0～100% 范围内，极端值通常由数值误差引起，属于正常现象。

相关课程推荐

• 使用 R 语言处理 Merra2 数据获取各省市区县的比湿、降水量、风速和气压数据
• 使用 R 语言处理 netCDF 格式的数据
• 1980年1月1日～2025年12月31日各省市区县相对湿度、比湿、地表高度、气压、风速和气温日度数据

点击这里跳转到 RStata 短书平台获取附件：如何使用 R 语言计算各省市区县相对湿度面板数据