如何使用孟加拉国shapefile从ERA5 NetCDF数据中提取臭氧均值？

核心思路是将读取到的四维臭氧数组（经度、纬度、气压层、时间）转换为空间栅格对象，和孟加拉国空间多边形对齐投影后，通过掩膜提取、区域统计得到对应范围的均值，替换原有的经纬度硬索引逻辑。

前置注意事项 #

• 需提前确认ERA5的经度范围：ERA5原始数据常使用0_{360°的经度编码，若你的shp用的是-180}180°的WGS84坐标系，需先将ERA5的经度转换到一致范围
• 优先使用terra包替代原有的raster包，处理多维栅格效率更高、语法更简洁

修改后完整代码 #

# 加载所需包
library(terra)
library(ncdf4)

# 1. 读取孟加拉国空间多边形，转换为terra的矢量对象，统一投影为WGS84
south_asia <- vect('E:/Bangladesh.shp')
south_asia <- project(south_asia, "EPSG:4326")

# 2. 读取ERA5 nc文件
nc <- nc_open("E:/ERA5_2000_03.nc")
lat <- ncvar_get(nc, "lat")
lon <- ncvar_get(nc, "lon")
t <- ncvar_get(nc, "time")
pres <- ncvar_get(nc, "plev")
o3 <- ncvar_get(nc,'o3')
nc_close(nc) # 读取完成后关闭连接，避免资源占用

# 3. 校准经纬度顺序和范围，保证栅格方向正确
# 若ERA5经度为0~360范围，取消下一行注释转换为-180~180
# lon <- ifelse(lon > 180, lon - 360, lon)
# 按升序排列经度
lon_sort_idx <- order(lon)
lon <- lon[lon_sort_idx]
o3 <- o3[lon_sort_idx,,,]
# 按升序排列纬度（ERA5默认从北到南存储，反转后适配栅格规则）
lat_sort_idx <- order(lat)
lat <- lat[lat_sort_idx]
o3 <- o3[,lat_sort_idx,,]

# 4. 构建多维栅格对象
# 调整数组维度顺序：将[lon,lat,pres,time]转为栅格要求的[lat,lon,层]
o3_perm <- aperm(o3, c(2,1,3,4))
dim(o3_perm) <- c(length(lat), length(lon), length(pres)*length(t))
# 生成栅格并配置空间属性
o3_rast <- rast(o3_perm,
                extent = ext(min(lon), max(lon), min(lat), max(lat)),
                crs = "EPSG:4326")
# 给栅格层命名，方便后续匹配气压层和时间
layer_names <- expand.grid(pres = pres, time = as.character(as.Date(t)))
names(o3_rast) <- paste0(layer_names$pres, "_", layer_names$time)

# 5. 按多边形提取区域均值
# 先裁剪到多边形范围提升运算效率
o3_crop <- crop(o3_rast, south_asia)
# 掩膜处理，仅保留多边形范围内的像素
o3_mask <- mask(o3_crop, south_asia)
# 统计每个层的区域均值，自动忽略空值
o3_mean <- global(o3_mask, fun = "mean", na.rm = TRUE)

# 6. 整理为和原代码输出完全一致的结构：行对应气压层，列对应时间
p <- matrix(o3_mean$mean, nrow = length(pres), ncol = length(t),
            dimnames = list(pres = pres, time = as.character(as.Date(t))))

代码说明 #

• 完全替换了原有硬编码经纬度索引的逻辑，所有范围筛选严格遵循空间多边形的实际边界
• 输出的p结构和你原有代码的输出完全一致，可直接对接后续分析流程
• 若臭氧数据存在缺失值，可调整global函数的na.rm参数控制是否忽略空值

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