基于R语言计算四维NetCDF海流变量的时间平均值并生成栅格

所需R包 #

先安装并加载必要的工具包：

install.packages(c("ncdf4", "raster", "dplyr"))
library(ncdf4)
library(raster)
library(dplyr)

步骤1：读取NetCDF数据 #

打开目标文件并提取核心数据与维度信息：

# 替换为你的NetCDF文件路径
currents_nc <- nc_open("your_currents_file.nc")

# 提取经纬度、深度、时间维度数据
lon <- ncvar_get(currents_nc, "longitude")
lat <- ncvar_get(currents_nc, "latitude")
depth <- ncvar_get(currents_nc, "depth")

# 提取海流变量（维度顺序：经度、纬度、深度、时间）
uo_data <- ncvar_get(currents_nc, "uo")
vo_data <- ncvar_get(currents_nc, "vo")

# 关闭文件连接，释放资源
nc_close(currents_nc)

步骤2：计算时间平均值 #

针对每个深度层，计算所有时间步的海流均值：

# 对uo沿时间维度（第4维）求平均，保留经纬度和深度维度
uo_time_mean <- apply(uo_data, MARGIN = c(1,2,3), FUN = mean, na.rm = TRUE)
# 对vo执行同样操作
vo_time_mean <- apply(vo_data, MARGIN = c(1,2,3), FUN = mean, na.rm = TRUE)

步骤3：转换为栅格对象 #

将各深度层的平均海流数据转换为栅格，最终生成包含所有深度层的栅格栈：

# 创建栅格模板，匹配数据的经纬度范围与坐标系
raster_template <- raster(nrows = length(lat), ncols = length(lon),
                          xmn = min(lon), xmx = max(lon),
                          ymn = min(lat), ymx = max(lat),
                          crs = "+proj=longlat +datum=WGS84")

# 构建uo的深度层栅格栈
uo_raster_stack <- stack()
for (i in 1:length(depth)) {
  # 转置数据适配栅格的行列逻辑（栅格行对应纬度，列对应经度）
  layer_data <- t(uo_time_mean[,,i])
  layer_raster <- setValues(raster_template, layer_data)
  names(layer_raster) <- paste0("uo_depth_", round(depth[i], 1))
  uo_raster_stack <- addLayer(uo_raster_stack, layer_raster)
}

# 构建vo的深度层栅格栈
vo_raster_stack <- stack()
for (i in 1:length(depth)) {
  layer_data <- t(vo_time_mean[,,i])
  layer_raster <- setValues(raster_template, layer_data)
  names(layer_raster) <- paste0("vo_depth_", round(depth[i], 1))
  vo_raster_stack <- addLayer(vo_raster_stack, layer_raster)
}

步骤4：保存栅格文件 #

将栅格栈导出为GeoTIFF格式，方便后续与水深数据叠加分析：

# 保存uo时间平均栅格
writeRaster(uo_raster_stack, "uo_time_mean_depth_layers.tif", format = "GTiff", overwrite = TRUE)
# 保存vo时间平均栅格
writeRaster(vo_raster_stack, "vo_time_mean_depth_layers.tif", format = "GTiff", overwrite = TRUE)

注意事项 #

• 若数据存在缺失值，mean函数中的na.rm = TRUE会自动忽略，可根据实际需求调整
• 若后续叠加的水深数据坐标系不同，需用projectRaster统一坐标系
• 处理超大文件时，可使用ncvar_get的start和count参数分块读取，避免内存溢出

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