Метаданные теряются при чтении файла GRIB2

Я хотел бы прочитать файл 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, а затем считываю файл NetCDF с помощью ncdf4 — дополнительный шаг преобразования, который я не хочу включать в свой рабочий процесс, если его можно избежать — есть определенно присутствуют имена переменных.

# [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, показанного выше, если бы порядок слоев изменился, это сломало бы мой пакет.


r raster rgdal grib
person Lyngbakr    schedule 01.12.2020    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Вы можете использовать пакет 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"

Я вижу, что УОГРД и ВОГРД имеют L160 и имеют число 160 в GRIB_PDS_TEMPLATE_ASSEMBLED_VALUE и GRIB_PDS_TEMPLATE_NUMBERS, где у остальных 1.

Структура метаданных описана здесь, но я мог бы использовать некоторые рекомендации о том, где искать, чтобы понять что извлечь из метаданных.

person Robert Hijmans    schedule 01.12.2020
comment
Спасибо, Роберт, это чрезвычайно полезно. Я не обращал внимания на пакет terra, но, похоже, это именно то, что мне нужно. Я не совсем уверен, как появились имена, которые я показываю, поскольку они были созданы ncl_convert2nc, но имена переменных, которые вы извлекаете из GRIB_ELEMENT, на самом деле то, что мне нужно. - person Lyngbakr; 02.12.2020
comment
убедитесь, что у вас есть самая последняя версия terra для этого --- она ​​только что попала в CRAN, но все еще нуждается в компиляции (или подождите пару дней, пока CRAN сделает это) - person Robert Hijmans; 02.12.2020