#' Read a field from a grib file
#'
#' @param filename The grib file name.
#' @param parameter The parameter to read. Standard HARP names are used.
#' @param ... Arguments for \code{get_grib_param_info} for level and levtype.
#'
#' @return A geofield object with 2d array and projection information
#' @export
#'
#' @examples
#' file_name <- "/lustre/storeB/users/andrewts/mepsr_data/grib/fc2017052600+001grib_fp"
#' model_geofield <- read_grib(file_name, "t2m")
#' model_geofield <- read_grib(file_name, "t", level = 500)
#' model_geofield <- read_grib(file_name, "topo")
read_grib <- function(filename, parameter, ...) {
if (!requireNamespace("Rgrib2", quietly = TRUE)) {
stop("Rgrib2 must be installed to read grib files")
}
param_info <- get_grib_param_info(parameter, ...)
if (is.na(param_info$short_name)) {
stop("Unknown parameter: ", parameter)
}
grib_info <- Rgrib2::Gopen(filename)
grib_position <- grib_info %>%
dplyr::filter(
shortName == param_info$short_name,
indicatorOfParameter == param_info$param_number,
indicatorOfTypeOfLevel == param_info$level_type,
level %in% param_info$level_number
) %>%
dplyr::arrange(.data$level)
if (length(unique(grib_position$shortName)) > 1) stop("More than one shortName found")
if (length(unique(grib_position$indicatorOfParameter)) > 1) stop("More than one indicatorOfParameter found")
if (length(unique(grib_position$indicatorOfTypeOfLevel)) > 1) stop("More than one indicatorOfTypeOfLevel found")
grib_position <- dplyr::select(grib_position, level, position)
if (nrow(grib_position) == 1) {
Rgrib2::Gdec(filename, grib_position$position)
} else {
cat("Reading:", filename, "\n")
show_level <- function(.level, num_levels) {
cat("Level:", .level, "of", num_levels, "\r")
flush.console()
}
all_levels <- purrr::map2(
grib_position$level,
grib_position$position,
~{show_level(.x, nrow(grib_position)); Rgrib2::Gdec(filename, .y)}
)
cat("\n")
# Note: attributes are dropped - wait for new meteogrid for multidimension meteogrid objects.
# In the meantime return a list with domain information - read_members will only ever ask for one level.
domain_data <- meteogrid::DomainExtent(all_levels[[1]])
x <- seq(domain_data$x0, domain_data$x1, domain_data$dx)
y <- seq(domain_data$y0, domain_data$y1, domain_data$dy)
proj4_string <- paste0(
"+", paste(
meteogrid::proj4.list2str(attr(all_levels[[1]], "domain")$projection), collapse = " +"
)
)
list(
model_data = simplify2array(all_levels),
x = x,
y = y,
z = param_info$level_number,
member = NA,
proj4_string = proj4_string,
parameter = parameter,
filename = filename,
dimensions = c("x", "y", grib_level_types$description[grib_level_types$id == param_info$level_type])
)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.