如何用R的terra::rast()将多份AQUA-MODIS分箱NetCDF转为SpatRaster

从Earthdata网站下载了2012至2024年埃及北部红海区域的3422份AQUA-MODIS叶绿素a（chlorophyll-a）NetCDF分箱文件。

使用terra::rast()（或等效函数）堆叠所有文件生成SpatRaster，计算所有文件的叶绿素a平均水平。

尝试过ncdf4、terra、raster、RNetCDF等包均未成功，具体操作及报错如下：

1. 读取文件代码：

library(terra)
library(ncdf4)
library(raster)
library(RNetCDF)

# 设置文件路径
folder<-"~/Documents/GIS_Data/CHL-A/"

# 获取文件夹内所有NC文件列表
files <-list.files(folder, pattern='*.nc', full.names ="TRUE")

2. 遍历文件代码：

for (file in files) {
     nc=open.nc(file)
     print.nc(nc)
}

3. 堆叠尝试及错误：

• 尝试1：chla1<-terra::rast(files)，报错：Error: [rast] cannot open this file as a SpatRaster: /Documents/GIS_Data/CHL-A/AQUA_MODIS.20120101.L3b.DAY.CHL.x.nc，附加警告：文件格式不被支持（GDAL error 4）。
• 尝试2：chla1<-terra::rast(nc)，报错：Error in methods::as(x, "SpatRaster") : no method or default for coercing “NetCDF” to “SpatRaster”。
• 尝试3：chl4<-raster::stack(files)，报错：Error in ncvar_type_to_string(rv$precint) : Error, unrecognized type code of variable supplied: -1。

解决步骤 #

这些AQUA-MODIS的L3b海洋NetCDF文件属于NASA Ocean Color专用格式，GDAL默认加载时无法自动识别目标变量，需明确指定叶绿素a变量名（通常为chlor_a，可通过print.nc(nc)输出确认），以下是可行方案：

方案1：指定变量名用terra批量加载堆叠 #

library(terra)

# 设置路径并获取文件列表
folder <- "~/Documents/GIS_Data/CHL-A/"
files <- list.files(folder, pattern = "\\.nc$", full.names = TRUE)

# 定义单文件加载函数
load_chla <- function(file) {
  # 明确指定读取'chlor_a'变量
  r <- rast(file, var = "chlor_a")
  # 提取文件名中的日期并设置为栅格时间属性
  date_str <- sub("AQUA_MODIS\\.(\\d{8})\\..*", "\\1", basename(file))
  time(r) <- as.Date(date_str, format = "%Y%m%d")
  return(r)
}

# 批量加载并堆叠为SpatRaster
chla_stack <- lapply(files, load_chla) |> rast()

# 计算时间序列平均值
chla_mean <- mean(chla_stack, na.rm = TRUE)

方案2：排查GDAL NetCDF驱动支持 #

若仍报错，先检查GDAL是否启用NetCDF驱动：

gdalDrivers() |> subset(name == "NetCDF")

若结果为空，需重新安装带NetCDF支持的GDAL版本：

• Linux：通过包管理器安装libgdal-dev和libnetcdf-dev后重装terra
• Windows/macOS：使用conda安装gdal和terra，或直接安装CRAN提供的预编译二进制包

方案3：用ncdf4手动读取构建SpatRaster #

如果上述方法失效，可手动读取变量和地理信息构建栅格：

library(terra)
library(ncdf4)

load_chla_manual <- function(file) {
  nc <- nc_open(file)
  # 读取叶绿素a数据、经纬度
  chla_data <- ncvar_get(nc, "chlor_a")
  lon <- ncvar_get(nc, "lon")
  lat <- ncvar_get(nc, "lat")
  nc_close(nc)
  
  # 构建栅格并调整维度匹配地理坐标系
  r <- rast(chla_data, 
            xmin = min(lon), xmax = max(lon), 
            ymin = min(lat), ymax = max(lat), 
            crs = "EPSG:4326")
  r <- t(r) |> flip("y")
  
  # 设置时间属性
  date_str <- sub("AQUA_MODIS\\.(\\d{8})\\..*", "\\1", basename(file))
  time(r) <- as.Date(date_str, format = "%Y%m%d")
  return(r)
}

chla_stack <- lapply(files, load_chla_manual) |> rast()
chla_mean <- mean(chla_stack, na.rm = TRUE)

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