如何用R的terra::extract()从多NetCDF文件提取特定日期多年SST点数据

问题描述 #

从NASA海洋色彩项目下载了4640个L3m 4km分辨率的海表温度（SST）NetCDF文件，单个文件对应的SpatRaster信息如下：

class       : SpatRaster 
dimensions  : 766, 709, 1  (nrow, ncol, nlyr)
resolution  : 0.04165021, 0.04165796  (x, y)
extent      : 24.61, 54.14, 5.8, 37.71  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : AQUA_MODIS.20120101.L3m.DAY.SST.x_sst.nc:sst 
varname     : sst (Sea Surface Temperature) 
name        :      sst 
unit        : degree_C 

数据时间跨度为2012至2024年，覆盖红海82km区域。需求是提取特定日期的海表温度点数据，但现有代码存在两个核心问题：

• 仅提取了第一个NetCDF文件的数据，未遍历全部文件；
• 无法将数据框中的Dates列与对应日期的SST点数据关联。

原代码如下：

library(terra)
library(lubricate)

#Download Sea Surface Temperature 
folder_SST<="~/Documents/GIS_Data/SST/requested_files"
files_SST <-list.files(folder_SST, pattern='*.nc', full.names ="TRUE")

#Loop through all the file paths and the AQUA MODIS netcdf4 files to gain access to them:
for (file in files_SST) {
  nc_SST=open.nc(files_SST)
  print.nc(nc_SST)
}

#Select the sea surface temperature values from AQUA MODIS files 
SSTs <- terra::rast(files_SST[1], "sst")

#Make the dataframe a spatial object of class = "sf" with a CRS of 4326
Ds_Points <- st_as_sf(x=MyDf, 
                      coords = c("Longitude_E_DD", "Latitude_N_DD"), 
                      crs = 4326)

#Format the dates in the dataframe 
#Get dates
MyDf <- MyDf %>% mutate(Dates = as.Date(as.character(MyDf$Date), format = "%d/%m/%Y"))

#Extract the SST values from the AQUA MODIS files 
SSTs_Data <- terra::extract(SSTs, Ds_Points)

期望输出格式：

Survey_Number       Dates Longitude_E_DD Latitude_N_DD    SST
              1   2012-08-01       33.89083      27.26778 23.635
              2   2012-06-02       33.86782      27.40854 23.640
              3   2012-02-07       33.86230      27.44623 23.690
              4   2012-02-12       33.88653      27.26957 23.635
              5   2012-02-13       33.88766      27.26848 23.635
              6   2012-02-14       33.85000      27.36111 23.780
              7   2012-02-15       33.86177      27.41302 23.640

解决方案 #

原代码问题分析 #

1. 循环逻辑错误：open.nc(files_SST)应改为open.nc(file)，且该循环仅打印文件信息，未处理数据；
2. 仅加载第一个文件：files_SST[1]限制了只读取第一个NetCDF文件；
3. 操作顺序颠倒：先创建sf对象再处理日期，导致日期更新未同步到空间对象；
4. 包名拼写错误：lubricate应为lubridate。

优化后代码 #

核心思路：从文件名提取日期，匹配观测点的Dates列，按需加载对应日期的SST文件，避免遍历全部4640个文件以提升效率。

library(terra)
library(lubridate)
library(sf)
library(dplyr)
library(stringr)

# 设置SST文件文件夹路径
folder_SST <- "~/Documents/GIS_Data/SST/requested_files"
files_SST <- list.files(folder_SST, pattern = "*.nc", full.names = TRUE)

# 构建文件路径与日期的映射表
file_date_map <- tibble(
  file_path = files_SST,
  # 从文件名提取8位日期字符串（匹配AQUA_MODIS.20120101...格式）
  date_str = str_extract(basename(file_path), "\\d{8}"),
  file_date = as.Date(date_str, format = "%Y%m%d")
)

# 预处理观测点数据：先格式化日期，再转为sf空间对象
MyDf <- MyDf %>%
  mutate(Dates = as.Date(as.character(Date), format = "%d/%m/%Y")) %>%
  st_as_sf(coords = c("Longitude_E_DD", "Latitude_N_DD"), crs = 4326)

# 按日期分组匹配并提取SST值
final_result <- MyDf %>%
  group_by(Dates) %>%
  group_map(function(group_data, group_key) {
    # 获取当前日期对应的SST文件路径
    target_file <- file_date_map %>%
      filter(file_date == group_key$Dates) %>%
      pull(file_path)
    
    # 无对应文件时，SST赋值为NA
    if (length(target_file) == 0) {
      return(group_data %>% mutate(SST = NA_real_) %>% st_drop_geometry())
    }
    
    # 加载对应日期的SST栅格并提取点值
    sst_raster <- terra::rast(target_file, "sst")
    sst_extract <- terra::extract(sst_raster, group_data)
    
    # 合并结果并移除几何列
    group_data %>%
      mutate(SST = sst_extract$sst) %>%
      st_drop_geometry()
  }) %>%
  bind_rows()

# 调整列顺序以匹配期望输出
final_result <- final_result %>%
  select(Survey_Number, Dates, Longitude_E_DD, Latitude_N_DD, SST)

# 查看最终结果
print(final_result)

关键说明 #

• 日期匹配：通过str_extract从文件名提取日期，实现观测点日期与SST文件的精准匹配；
• 按需加载：仅加载观测点对应日期的SST文件，避免不必要的内存占用；
• 缺失值处理：若某日期无对应SST文件，自动将该日期的SST值设为NA；
• 格式对齐：最后调整列顺序，与期望输出完全匹配。

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