我想将GRIB2文件读入R,但无法安装wgrib2
(经过几个小时的努力(,这意味着rNOMADS
不是一个选项。这没关系,因为GRIB2文件可以由raster
和rgdal
包读取。我遇到的问题是,在读取文件时,层的名称会被剥离。
下面是一个例子。
# 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我需要从文件中获取变量名的原因是,当加载(例如(raster
时,层与网站上指定的层顺序不匹配。从变量值中可以明显看出这一点。虽然我可以使用从上面显示的NetCDF文件中收集的变量名,但如果层的顺序发生变化,这将破坏我的包。
您可以使用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"
我可以改变CCD_ 12的行为;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";具有CCD_ 17并且具有"1"中的数字160;GRIB_PDS_TEMPLATE_ASSEMBLED_VALUE";以及";GRIB_ PDS_TEMPLATE_;其中其他具有1。
这里描述了元数据结构,但我可以使用一些关于在哪里查找的指导来理解从元数据中提取什么。