读取 GRIB2 文件时元数据丢失

Question

我想将 GRIB2 文件读入 R，但无法安装 wgrib2（经过几个小时的努力），这意味着 rNOMADS 不是一个选项。没关系，因为raster 和rgdal 包都可以读取GRIB2 文件。我遇到的问题是在读取文件时层的名称被剥离。

这是一个例子。

# Load libraries
library(raster)
library(rgdal)

# Name of file
file_name <- "https://dd.weather.gc.ca/model_gem_regional/coupled/gulf_st-lawrence/grib2/00/001/CMC_coupled-rdps-stlawrence-ocean_latlon0.02x0.03_2020120100_P001.grib2"

# Load as raster brick
b <- brick(file_name)

# Get layer names
names(b)
# [1] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.1" 
# [2] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.2" 
# [3] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.3" 
# [4] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.4" 
# [5] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.5" 
# [6] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.6" 
# [7] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.7" 
# [8] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.8" 
# [9] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.9" 
#[10] "CMC_coupled.rdps.stlawrence.ocean_latlon0.02x0.03_2020120100_P001.10"

如您所见，名称只是通用默认值。接下来，我尝试了rgdal。

# Load using rgdal
r <- readGDAL(file_name)

# Get names
names(r)

# [1] "band1"  "band2"  "band3"  "band4"  "band5"  "band6"  "band7"  "band8" 
# [9] "band9"  "band10"

再一次，默认名称。但是，如果我使用命令行实用程序 ncl_convert2nc 将 GRIB2 文件转换为 NetCDF，然后使用 ncdf4 读入 NetCDF 文件 – 如果可以的话，我不想在我的工作流程中包含一个额外的转换步骤避免 - 肯定存在变量名称。

# [1] "UOGRD_P0_L160_GLL0" "VOGRD_P0_L160_GLL0" "ICEC_P0_L1_GLL0"   
# [4] "ICETK_P0_L1_GLL0"   "UICE_P0_L1_GLL0"    "VICE_P0_L1_GLL0"   
# [7] "ICETMP_P0_L1_GLL0"  "ICEPRS_P0_L1_GLL0"  "CICES_P0_L1_GLL0"  
#[10] "WTMP_P0_L1_GLL0"   

问题：在使用rgdal 或raster 读取GRIB2 文件时，有没有办法提取或保留变量/层名称？

---

PS 我需要从文件中获取变量名称的原因是，当加载（例如）@987654334 时，层 不 与指定的层顺序匹配 on the website @。从变量值可以看出这一点。虽然我可以使用从上面显示的 NetCDF 文件中收集的变量名称，但如果层的顺序发生更改，这会破坏我的包。

Answer 1

您可以使用terra 包代替raster。

file_name <- "https://dd.weather.gc.ca/model_gem_regional/coupled/gulf_st-lawrence/grib2/00/001/CMC_coupled-rdps-stlawrence-ocean_latlon0.02x0.03_2020120100_P001.grib2"
b <- basename(file_name)
if (!file.exists(b)) download.file(file_name, b, mode="wb")
library(terra)
r <- rast(b)
r
#class       : SpatRaster 
#dimensions  : 325, 500, 10  (nrow, ncol, nlyr)
#resolution  : 0.03, 0.02  (x, y)
#extent      : -71.015, -56.015, 45.49, 51.99  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +R=6371229 +no_defs 
#source      : CMC_coupled-rdps-stlawrence-ocean_latlon0.02x0.03_2020120100_P001.grib2 
#names       : 0[-] SFC="Ground or water surface", 0[-] SFC="Ground or water surface", 0[-] SFC="Ground or water surface", 0[-] SFC="Ground or water surface", 0[m] DBSL="Depth below sea level", 0[m] DBSL="Depth below sea level", ... 

但是变量名和你的不匹配

names(r)
# [1] "0[-] SFC=\"Ground or water surface\"" "0[-] SFC=\"Ground or water surface\"" "0[-] SFC=\"Ground or water surface\""
# [4] "0[-] SFC=\"Ground or water surface\"" "0[m] DBSL=\"Depth below sea level\""  "0[m] DBSL=\"Depth below sea level\"" 
# [7] "0[-] SFC=\"Ground or water surface\"" "0[-] SFC=\"Ground or water surface\"" "0[-] SFC=\"Ground or water surface\""
#[10] "0[-] SFC=\"Ground or water surface\""

您可以将名称设置为元数据的其他部分

 nms <- trimws(grep("GRIB_ELEMENT=", desc(b), value=TRUE))
 names(r) <- gsub("GRIB_ELEMENT=", "", nms)

r
#class       : SpatRaster 
#dimensions  : 325, 500, 10  (nrow, ncol, nlyr)
#resolution  : 0.03, 0.02  (x, y)
#extent      : -71.015, -56.015, 45.49, 51.99  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +R=6371229 +no_defs 
#source      : CMC_coupled-rdps-stlawrence-ocean_latlon0.02x0.03_2020120100_P001.grib2 
#names       : ICEC, ICETK, UICE, VICE, UOGRD, VOGRD, ..

names(r)
#[1] "ICEC"   "ICETK"  "UICE"   "VICE"   "UOGRD"  "VOGRD"  "WTMP"   "ICET"   "ICEPRS" "CICES" 

我可以更改terra 的行为，使其使用“GRIB_ELEMENT”（如果有意义，请告诉我）。但我不清楚如何获得您显示的名称。例如，下面是第一层的 GDAL 元数据。你显示ICEC_P0_L1_GLL0。所有名称都有P0 和GLL0，所以至少对于这个文件来说，这些似乎是多余的。但是L1 指的是什么？

d <-desc(b) 
d[35:46]
# [1] "    GRIB_COMMENT=Ice cover [Proportion]"                                                                                                                                                
# [2] "    GRIB_DISCIPLINE=10(Oceanographic_Products)"                                                                                                                                         
# [3] "    GRIB_ELEMENT=ICEC"                                                                                                                                                                  
# [4] "    GRIB_FORECAST_SECONDS=3600 sec"                                                                                                                                                     
# [5] "    GRIB_IDS=CENTER=54(Montreal) SUBCENTER=0 MASTER_TABLE=4 LOCAL_TABLE=0 SIGNF_REF_TIME=1(Start_of_Forecast) REF_TIME=2020-12-01T00:00:00Z PROD_STATUS=0(Operational) TYPE=1(Forecast)"
# [6] "    GRIB_PDS_PDTN=0"                                                                                                                                                                    
# [7] "    GRIB_PDS_TEMPLATE_ASSEMBLED_VALUES=2 0 2 56 56 0 0 1 1 1 0 0 255 -127 -2147483647"                                                                                                  
# [8] "    GRIB_PDS_TEMPLATE_NUMBERS=2 0 2 56 56 0 0 0 1 0 0 0 1 1 0 0 0 0 0 255 255 255 255 255 255"                                                                                          
# [9] "    GRIB_REF_TIME=  1606780800 sec UTC"                                                                                                                                                 
#[10] "    GRIB_SHORT_NAME=0-SFC"                                                                                                                                                              
#[11] "    GRIB_UNIT=[Proportion]"                                                                                                                                                             
#[12] "    GRIB_VALID_TIME=  1606784400 sec UTC"    

我看到“UOGRD”和“VOGRD”有L160，并且在“GRIB_PDS_TEMPLATE_ASSEMBLED_VALUE”和“GRIB_PDS_TEMPLATE_NUMBERS”中的数字是160，而其他的数字是1。

元数据结构在here 中进行了描述，但我可以使用一些指导来了解从元数据中提取的内容。