如何在R中按日期从netCDF文件提取对应位置的SST值？

需求背景 #

我是R语言新手，现有一批NetCDF格式的每日海表温度（SST）文件，文件名格式为TERRA_MODIS.20191201.L3m.DAY.SST.x_sst.nc，其中20191201是yyyymmdd格式的日期。同时我有一个包含浮标编号、日期、经纬度的数据框，需要为每一行匹配对应日期NetCDF文件中对应坐标的SST值，添加为SST列。

数据框示例：

example <- read.csv('example.csv', sep=';', dec='.')
head(example)
#   Buoy       Date Longitude Latitude
# 1     1 2019-12-01  50.29614 43.92681
# 2     1 2019-12-02  50.23525 43.89244
# 3     1 2019-12-03  50.19717 43.88295

解决方案步骤 #

1. 加载依赖包并设置文件路径 #

首先加载处理NetCDF和数据操作的包，指定存放NetCDF文件的目录：

library(ncdf4)
library(dplyr)

# 设置NetCDF文件所在目录，替换为你的实际路径
nc_dir <- "./sst_netcdf_files/"

2. 标准化数据框的日期格式 #

将数据框中的Date列转换为日期格式，并生成yyyymmdd格式的字符串，用于匹配文件名：

example <- example %>%
  mutate(
    Date = as.Date(Date),
    date_str = format(Date, "%Y%m%d") # 生成与文件名匹配的日期字符串
  )

3. 编写提取SST的函数 #

编写一个函数，输入日期字符串、经度、纬度，返回对应NetCDF文件中最接近该坐标的SST值：

extract_sst <- function(date_str, lon_target, lat_target, nc_dir) {
  # 拼接对应日期的NetCDF文件名
  nc_filename <- paste0("TERRA_MODIS.", date_str, ".L3m.DAY.SST.x_sst.nc")
  nc_path <- file.path(nc_dir, nc_filename)
  
  # 检查文件是否存在
  if (!file.exists(nc_path)) {
    warning(paste("文件不存在：", nc_path))
    return(NA)
  }
  
  # 打开NetCDF文件
  nc <- nc_open(nc_path)
  
  # 获取经纬度和SST数据
  lon <- ncvar_get(nc, "lon")
  lat <- ncvar_get(nc, "lat")
  sst_data <- ncvar_get(nc, "sst")
  
  # 关闭NetCDF文件（必须执行，避免资源占用）
  nc_close(nc)
  
  # 找到最接近目标经纬度的网格索引
  lon_idx <- which.min(abs(lon - lon_target))
  lat_idx <- which.min(abs(lat - lat_target))
  
  # 返回对应位置的SST值
  return(sst_data[lon_idx, lat_idx])
}

4. 批量提取SST并添加到数据框 #

使用rowwise()逐行处理数据框，调用上述函数提取SST：

example <- example %>%
  rowwise() %>%
  mutate(
    SST = extract_sst(date_str, Longitude, Latitude, nc_dir)
  ) %>%
  ungroup() %>%
  select(-date_str) # 移除临时的date_str列

5. 查看结果 #

处理完成后查看数据框：

print(example)

输出示例：

Buoy Date       Longitude Latitude   SST
   <int> <date>        <dbl>    <dbl> <dbl>
 1     1 2019-12-01     50.3     43.9    13
 2     1 2019-12-02     50.2     43.9    16
 3     1 2019-12-03     50.2     43.9     2
 4     1 2019-12-04     50.2     43.9    10
 5     1 2019-12-05     50.3     43.9     8
 6     2 2019-12-01     51.7     46.5    18
 7     2 2019-12-02     51.8     46.6     4
 8     2 2019-12-03     51.8     46.5    17
 9     2 2019-12-04     51.8     46.6    20
10     2 2019-12-05     51.9     46.3    13
11     3 2019-12-01     51.5     41.3    14
12     3 2019-12-02     50.8     41.6    18
13     3 2019-12-03     51.4     41.8     8
14     3 2019-12-04     49.5     42.3     7
15     3 2019-12-05     49.5     42.6     2

注意事项 #

• 确保NetCDF文件目录路径正确，且所有需要的文件都存在
• 如果浮标坐标不在NetCDF网格点上，上述函数会取最接近的网格点的SST值；如果需要更精确的插值（如双线性插值），可以用akima包的interp()函数实现
• 处理大量文件时，逐行读取可能较慢，可考虑先按日期分组，批量读取文件后再匹配所有对应日期的浮标数据，提升效率

内容的提问来源于stack exchange，提问作者Bowe22