如何用R提取NetCDF文件中各网格单元的最大底层深度数据并计算多年均值

我有一个包含lon、lat、time、depth四个维度，以及两个海洋环境变量的全球NetCDF文件。我需要提取每个非NA网格单元最大底层深度对应的数据，并计算其多年均值。目前我的代码只能处理表层（level=1，对应-1m深度）的数据，不知道如何实现最大底层深度的提取，请问如何用R（而非CDO）完成该操作？

当前代码如下：

library(raster)
library(ncdf4)

setwd("E:/Modélisation - Eric Goberville/PhD/Environmental data/Raw data/Western_MedSea/Merged")

# CMCC files
rm(list=ls())
graphics.off()

# Load .nc file and access information
Info.Data <- ncdf4::nc_open("med-cmcc-tem-rean-m_01_2000_01_2021.nc")
Data_date <- ncvar_get(Info.Data, "time") # 查看时间范围
head(Data_date)
# Data_date <- as.POSIXct("1900-01-01 00:00:00", tz = "UTC") + as.difftime(Data_date, units = "mins")# 转换为标准日期格式
print(Data_date)
Data_depth <- ncvar_get(Info.Data, "depth") # 获取深度层信息

# Extract data
names(Info.Data[["var"]])
Data_raster <- brick("med-cmcc-tem-rean-m_01_2000_01_2021.nc", varname = "bottomT", level = 1)
# Data_raster

# Calculate annual mean from monthly means
# 快速方法，已知时间范围
# 创建索引矩阵，按年分组（每月一组，每年12个月）
ID <- cbind(seq(1, dim(Data_date), by = 12), seq(12, dim(Data_date), by = 12))

Data_raster_annual <- stack() # 空栈存储年平均值
for (i in 1:nrow(ID)) {
  Data_raster_annual <- stack(Data_raster_annual,
                              calc(Data_raster[[ID[i,1]:ID[i,2]]], mean)) # 计算当前12个月的单元格均值
}
# Calculate climatology (from annual means)
Raster_climatology <- calc(Data_raster_annual, mean)
# Raster_climatology

# Mapping of your environmental parameters 
graphics.off()
my.colors = colorRampPalette(c("#5E85B8","#EDF0C0","#C13127")) 
plot(Raster_climatology, main = "Bottom temperature",col=my.colors(1000),axes=FALSE, box=FALSE) 

由于NetCDF是4维数据（lon, lat, time, depth），raster包的brick无法直接处理，因此我们需要用ncdf4读取完整数组后，逐网格识别最大有效深度层，提取对应数据再转换为栅格进行均值计算。以下是完整实现代码：

library(raster)
library(ncdf4)

setwd("E:/Modélisation - Eric Goberville/PhD/Environmental data/Raw data/Western_MedSea/Merged")

rm(list=ls())
graphics.off()

# 1. 读取NetCDF文件及维度信息
Info.Data <- ncdf4::nc_open("med-cmcc-tem-rean-m_01_2000_01_2021.nc")
lon <- ncvar_get(Info.Data, "lon")
lat <- ncvar_get(Info.Data, "lat")
Data_date <- ncvar_get(Info.Data, "time")
Data_depth <- ncvar_get(Info.Data, "depth")
# 获取各维度长度
n_lon <- length(lon)
n_lat <- length(lat)
n_time <- length(Data_date)
n_depth <- length(Data_depth)

# 2. 读取完整的4维bottomT变量（lon, lat, time, depth）
bottomT_array <- ncvar_get(Info.Data, "bottomT")
ncdf4::nc_close(Info.Data) # 关闭文件释放内存

# 3. 识别每个网格的最大有效深度层索引
# 假设Data_depth是从小到大排列（表层到深层），索引越大深度越深
max_depth_idx_matrix <- matrix(NA, nrow = n_lat, ncol = n_lon)

for (i in 1:n_lon) {
  for (j in 1:n_lat) {
    # 提取当前网格的所有时间-深度数据
    grid_data <- bottomT_array[i, j, , ]
    # 标记每个深度层是否存在有效数据（至少一个时间步非NA）
    has_valid_data <- apply(grid_data, 2, function(col) any(!is.na(col)))
    # 记录最深的有效深度层索引
    if (any(has_valid_data)) {
      max_depth_idx_matrix[j, i] <- max(which(has_valid_data))
    }
  }
}

# 4. 提取每个网格对应最大深度层的时间序列
bottomT_max_depth_matrix <- matrix(NA, nrow = n_lon * n_lat, ncol = n_time)
cell_idx <- 1

for (i in 1:n_lon) {
  for (j in 1:n_lat) {
    idx <- max_depth_idx_matrix[j, i]
    if (!is.na(idx)) {
      # 提取该深度层的所有时间步数据
      bottomT_max_depth_matrix[cell_idx, ] <- bottomT_array[i, j, , idx]
    }
    cell_idx <- cell_idx + 1
  }
}

# 5. 转换为RasterBrick，方便后续均值计算
raster_template <- raster(ncols = n_lon, nrows = n_lat, 
                          xmn = min(lon), xmx = max(lon), 
                          ymn = min(lat), ymx = max(lat), 
                          crs = "+proj=longlat +datum=WGS84")
bottomT_max_depth_brick <- brick(raster_template, nl = n_time)
values(bottomT_max_depth_brick) <- bottomT_max_depth_matrix

# 6. 计算年平均值与多年气候态均值
ID <- cbind(seq(1, n_time, by = 12), seq(12, n_time, by = 12))
Data_raster_annual <- stack()

for (i in 1:nrow(ID)) {
  annual_mean <- calc(bottomT_max_depth_brick[[ID[i,1]:ID[i,2]]], mean, na.rm = TRUE)
  Data_raster_annual <- stack(Data_raster_annual, annual_mean)
}

# 多年平均气候态
Raster_climatology <- calc(Data_raster_annual, mean, na.rm = TRUE)

# 7. 可视化
graphics.off()
my.colors <- colorRampPalette(c("#5E85B8","#EDF0C0","#C13127")) 
plot(Raster_climatology, main = "最大底层深度多年平均温度", col = my.colors(1000), axes = FALSE, box = FALSE)

关键说明 #

1. 4维数组处理：直接读取完整4维变量，避免raster包对维度的限制；
2. 最大深度层识别：逐网格检查每个深度层的有效性，取最深的非NA层索引；
3. 内存优化：若文件过大，可按经度/纬度分块循环处理，避免一次性加载全量数据；
4. 兼容性：转换为RasterBrick后，可复用你原有的均值计算与可视化逻辑。

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