R语言处理NC栅格数据：计算指定区域逐时间步长的参数平均值

错误原因分析 #

• 第一种方案错误：你通过cbind(0:30,40:70)生成的是31个离散经纬度点对，不是0_{30经度、40}70纬度覆盖的所有网格点，提取结果仅对应这31个点的逐时间步数值，转置后自然会出现多列结果，不符合区域平均需求。
• 第二种方案错误：使用多边形提取时你没有在extract步骤指定空间平均函数，提取结果是该多边形覆盖所有网格的逐时间步数值矩阵，后续使用lapply对整行求平均相当于对所有时间步的数值再次求平均，最终只能得到一个全局均值。

最优实现方案 #

推荐先裁剪栅格到目标区域，再直接对每个时间步（每个图层）求区域均值，运算效率最高：

library(ncdf4)
library(raster)
library(rgdal)

# 读取NC数据，level=1对应1000hPa等压面，按需调整即可
r_brick <- brick('hgt.2006.nc', "hgt", level = 1)

# 构造目标区域范围：经度0~30，纬度40~70
target_extent <- extent(c(0, 30, 40, 70))
# 裁剪栅格到目标区域
r_crop <- crop(r_brick, target_extent)

# 逐时间步计算区域平均值，得到长度为365的一维向量
hgt_mean <- cellStats(r_crop, stat = mean, na.rm = TRUE)

# 整理为数据框输出
dfx <- data.frame(day = 1:365, hgt = hgt_mean)

注意：如果后续处理跨180度经线的区域，需要先将NCEP的0_{360经度转换为-180}180经度再做裁剪，避免范围匹配错误

多边形方案修正版 #

如果需要沿用多边形提取逻辑，修改如下：

library(ncdf4)
library(raster)
library(rgdal)

r_brick <- brick('hgt.2006.nc', "hgt", level = 1)
# 构造闭合矩形多边形
cds <- rbind(c(0,40), c(0,70), c(30,70), c(30, 40), c(0,40))
polys <- spPolygons(cds)

# 提取时直接指定空间平均函数，返回每个多边形每个时间步的均值
v <- extract(r_brick, polys, fun = mean, na.rm = TRUE)
# 提取结果为1行365列矩阵，转换为向量即可
hgt_mean <- as.vector(v)

dfx <- data.frame(day = 1:365, hgt = hgt_mean)

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