如何用R读取含多维度的概率气象预报GRIB2数据集

问题背景 #

尝试用R读取一份包含以下维度的概率气象预报GRIB2文件：

• 2个变量
• 114个预报时效
• 18个经度、24个纬度
• 50个集合成员

最初通过Python脚本将GRIB文件转为NetCDF格式后再用R读取，但该方案依赖Python及特定包，无法在便携环境运行。转而尝试用R的terra包直接读取，但无法正确提取所有维度数据，尝试代码如下：

require(terra)

## Isn't it possible to get them automatically ?
lat_prev <- (rev(seq(42.875,48.625,by=0.25)))
lon_prev <- (seq(3.375,7.625,by=0.25))

latlon <- expand.grid(lon=lon_prev, lat=lat_prev)

latlons <- terra::vect(latlon, geom=c('lon','lat'), crs="+proj=longlat")
## Not sure about that...
pts <- project(latlons, "+proj=lcc +lat_0=38.5 +lon_0=262.5 +lat_1=38.5 +lat_2=38.5 +x_0=0 +y_0=0 +R=6371229 +units=m")

grib_data <- terra::rast(destfile_CF)
## gives a data frame of NaN
e1 <- extract(grib_data, pts)

核心疑问：需要提取所有维度数据，但并非所有维度都是空间点（SpatialVector类），应创建何种对象？

解决方案 #

已找到可行方案，代码如下：

require(terra)
require(dplyr)
require(data.table)
require(stats)

destfile <- "C:/Users/XXX/Documents/Example_file.grb"

## 读取GRIB文件
grib_data <- terra::rast(destfile)
print(grib_data)

## 转换为数据框
df <- as.data.frame(grib_data, xy=TRUE)

## 重命名列：成员(50) × 时间(57) × 变量(2)组合
colNames <- paste(names(grib_data), as.character(time(grib_data)), sep="_")
colnames(df) <- c("lon", "lat", colNames)

## 宽表转长表
df2 <- data.table::melt(as.data.table(df), c("lon", "lat"))

## 拆分变量与时间字段
df2$time_UTC <- sub(".*_", "", df2$variable)    
df2$variable <- sub("_.*", "", df2$variable)   

## 添加集合成员编号
df2 <- df2 %>% group_by(lon, lat, variable, time_UTC) %>% mutate(member=(1:length(value)))

## 转换为多维数组
df_array <- stats::xtabs(value~lon+lat+variable+member+time_UTC, df2, drop=F)

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