如何在R中从小时级NetCDF数据生成月度序列及栅格地图？

核心思路 #

避免一次性加载全年8760小时的全量数据，通过分时间块（或分月）读取数据，逐块计算月度总和与计数，最终得到每个像素的月度平均值。这种方法每次仅加载小批量数据（如单日24小时），内存压力与你已成功读取的单日数据一致。

具体步骤（使用ncdf4和raster包） #

1. 安装并加载依赖包 #

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

2. 打开NetCDF文件，获取维度与时间信息 #

# 打开文件（不读取全量数据）
nc <- nc_open("NEE_2022.nc")

# 获取空间维度
lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")

# 获取时间维度并解析为日期时间
time_units <- ncatt_get(nc, "time", "units")$value
time_datetime <- as.POSIXct(ncvar_get(nc, "time"), units = time_units, tz = "UTC")
# 提取每个时间点对应的月份（格式：YYYY-MM）
months <- format(time_datetime, "%Y-%m")
unique_months <- unique(months)

3. 分天读取数据，累加月度总和与计数 #

这种方式内存占用最小，仅加载单日24小时数据：

# 初始化存储月度总和、有效数据计数的数组
monthly_sum <- array(0, dim = c(length(lon), length(lat), length(unique_months)))
monthly_count <- array(0, dim = c(length(lon), length(lat), length(unique_months)))

# 遍历2022年每一天（平年365天）
for (day_idx in 0:364) {
  # 计算当天对应的时间索引范围
  start_time <- day_idx * 24 + 1
  end_time <- (day_idx + 1) * 24
  
  # 读取当天24小时的NEE数据
  nee_day <- ncvar_get(nc, "NEE", start = c(1, 1, start_time), count = c(-1, -1, 24))
  
  # 获取当天所属的月份（处理跨月边界情况）
  day_month <- unique(months[start_time:end_time])
  m_idx <- which(unique_months == day_month)
  
  # 累加当天数据的总和与有效计数
  monthly_sum[,,m_idx] <- monthly_sum[,,m_idx] + apply(nee_day, c(1,2), sum, na.rm = TRUE)
  monthly_count[,,m_idx] <- monthly_count[,,m_idx] + apply(nee_day, c(1,2), function(x) sum(!is.na(x)))
}

# 计算月度平均值（避免除以0）
monthly_avg <- monthly_sum / pmax(monthly_count, 1)

# 关闭NetCDF文件
nc_close(nc)

4. 生成月度栅格地图并保存 #

# 遍历每个月份，生成栅格并保存/绘图
for (i in seq_along(unique_months)) {
  # 创建栅格对象（注意：若原始纬度是北->南顺序，需翻转）
  r <- raster(t(monthly_avg[,,i]), 
              xmn = min(lon), xmx = max(lon), 
              ymn = min(lat), ymx = max(lat),
              crs = "+proj=longlat +datum=WGS84")
  # 翻转栅格（适配大多数GIS工具的纬度顺序）
  r <- flip(r, direction = "y")
  
  # 保存为GeoTIFF格式
  writeRaster(r, filename = paste0("NEE_monthly_avg_", unique_months[i], ".tif"), 
              format = "GTiff", overwrite = TRUE)
  
  # 绘制月度地图（可选）
  plot(r, main = paste("NEE Monthly Average -", unique_months[i]))
}

优化建议 #

• 如果分月读取数据（而非分天）内存足够（如单月700+小时数据），可直接读取整月数据后计算平均，代码更简洁：将步骤3替换为分月循环读取对应时间切片，直接用apply(nee_month, c(1,2), mean, na.rm = TRUE)得到月度平均。
• 若想加速计算，可替换apply为matrixStats包的rowMeans2/colSums2函数，效率更高。
• 确认时间变量的单位：若time_units不是默认的小时/天，需根据NetCDF文件的属性调整解析逻辑。

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