From f0bd55cbf5cc4812db85ed052de8891e01260602 Mon Sep 17 00:00:00 2001 From: brry Date: Mon, 30 May 2022 10:55:02 +0200 Subject: [PATCH] readDWD.grib2: reading with terra and stars is now possible. Closes #28 --- DESCRIPTION | 4 +-- R/readDWD.R | 69 ++++++++++++++++++++++++++++++++++---------- man/readDWD.grib2.Rd | 44 ++++++++++++++++++++-------- 3 files changed, 88 insertions(+), 29 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0a3d1c6..92989f6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,12 @@ Package: rdwd Title: Select and Download Climate Data from 'DWD' (German Weather Service) -Version: 1.5.32 +Version: 1.5.33 Date: 2022-05-30 Depends: R(>= 2.10) Imports: berryFunctions (>= 1.21.11), pbapply Suggests: RCurl, leaflet, knitr, rmarkdown, testthat, roxygen2, devtools, remotes, bit64, data.table, OSMscale, raster, R.utils, ncdf4, readr, dwdradar, - XML, sp, rgdal + XML, sp, rgdal, terra, stars Author: Berry Boessenkool Maintainer: Berry Boessenkool Description: Handle climate data from the 'DWD' ('Deutscher Wetterdienst', see diff --git a/R/readDWD.R b/R/readDWD.R index e376385..de60481 100644 --- a/R/readDWD.R +++ b/R/readDWD.R @@ -141,6 +141,8 @@ msg <- paste0("Reading ",length(file)," file", if(length(file)>1)"s", " with ",n if(any("data" %in% type)) msg <- paste0(msg, " and fread=",nt(fread,"","")) if(any(c("radar", "raster", "asc") %in% type)) msg <- paste0(msg, " and dividebyten=",nt(dividebyten,"","")) +if(any("grib2" %in% type)) + message(msg, appendLF=FALSE) else message(msg, " ...") } @@ -1195,49 +1197,86 @@ return(invisible(list(dat=rbmat, meta=rbmeta))) #' \cr #' @examples #' \dontrun{ # Excluded from CRAN checks, but run in localtests -#' # Deactivated 2021-04-08 since readDWD.grib2 -> rgdal::readGDAL -> Error: -#' # **.grib2 is a grib file, but no raster dataset was successfully identified. -#' warning("readDWD.grib2 is not tested due to unresolved problems.") -#' if(FALSE){ -#' nwp_t2m_base <- "ftp://opendata.dwd.de/weather/nwp/icon-d2/grib/03/p" +#' nwp_t2m_base <- "ftp://opendata.dwd.de/weather/nwp/icon-d2/grib/15/soiltyp" #' nwp_urls <- indexFTP("", base=nwp_t2m_base, dir=tempdir()) -#' nwp_file <- dataDWD(nwp_urls[6], base=nwp_t2m_base, dir=tempdir(), +#' # for p instead of soiltyp, icosahedral_model-level files fail with GDAL errors, +#' # see https://github.com/brry/rdwd/issues/28 +#' # regular-lat-lon_pressure-level files work with pack="terra" or "stars" +#' +#' nwp_file <- dataDWD(tail(nwp_urls,1), base=nwp_t2m_base, dir=tempdir(), #' joinbf=TRUE, dbin=TRUE, read=FALSE) -#' nwp_data <- readDWD(nwp_file, quiet=TRUE) -#' plotRadar(nwp_data, project=FALSE) +#' nwp_data <- readDWD(nwp_file) +#' terra::plot(nwp_data) # same map with sp::plot +#' addBorders() # the projection seems to be perfectly good :) #' -#' nwp_data_rgdal <- readDWD(nwp_file, toraster=FALSE) -#' sp::plot(nwp_data_rgdal) +#' # index of GRIB files +#' if(FALSE){ # indexing takes about 6 minutes! +#' grib_base <- "ftp://opendata.dwd.de/weather/nwp/icon-d2/grib" +#' grib_files <- indexFTP("", base=grib_base, dir=tempdir()) +#' for(f in unique(substr(grib_files, 1,3))) print(grib_files[which(substr(grib_files, 1,3)==f)[1]]) +#' View(data.frame(grep("regular",grib_files, value=TRUE))) #' } #' } #' @param file Name of file on harddrive, like e.g. #' cosmo-d2_germany_regular-lat-lon_single-level_2021010100_005_T_2M.grib2.bz2 +#' @param pack Char: package used for reading. +#' One of "terra" (the default), "stars" +#' or "rgdal" (for the deprecated cosmo-d2 data). +#' See [issue](https://github.com/brry/rdwd/issues/28). +#' DEFAULT: "terra" #' @param bargs Named list of arguments passed to #' [R.utils::bunzip2()], see `gargs` in [readDWD.raster()]. DEFAULT: NULL #' @param toraster Logical: convert [rgdal::readGDAL] output with [raster::raster()]? +#' Only used if pack="rgdal". #' DEFAULT: TRUE #' @param quiet Silence readGDAL completely, including warnings on #' discarded ellps / datum. #' DEFAULT: FALSE through [rdwdquiet()] -#' @param \dots Further arguments passed to [rgdal::readGDAL()], -readDWD.grib2 <- function(file, bargs=NULL, toraster=TRUE, quiet=rdwdquiet(), ...) +#' @param \dots Further arguments passed to [stars::read_stars()], +#' [rgdal::readGDAL()] or [rgdal::readGDAL()]. +readDWD.grib2 <- function( +file, +pack="terra", +bargs=NULL, +toraster=TRUE, +quiet=rdwdquiet(), +...) { checkSuggestedPackage("R.utils", "rdwd:::readDWD.grib2") -checkSuggestedPackage("rgdal" , "rdwd:::readDWD.grib2") # bunzip arguments: bdef <- list(filename=file, remove=FALSE, skip=TRUE) bfinal <- berryFunctions::owa(bdef, bargs, "filename") # unzip: bdata <- do.call(R.utils::bunzip2, bfinal) -# rgdal reading: + +# actual reading: +# TERRA --- +if(pack=="terra"){ +if(!quiet) message(" with pack='terra' ...") +checkSuggestedPackage("terra" , "rdwd:::readDWD.grib2") +out <- terra::rast(bdata, ...) +} else +# STARS --- +if(pack=="stars"){ +if(!quiet) message(" with pack='stars' ...") +checkSuggestedPackage("stars" , "rdwd:::readDWD.grib2") +out <- stars::read_stars(bdata, quiet=quiet, ...) +} else +# RGDAL --- +if(pack=="rgdal"){ +if(!quiet) message(" with pack='rgdal' ...") +checkSuggestedPackage("rgdal" , "rdwd:::readDWD.grib2") out <- if(!quiet) rgdal::readGDAL(bdata, ...) else suppressWarnings(rgdal::readGDAL(bdata, silent=TRUE, ...)) -# conversion to raster: if(toraster) { checkSuggestedPackage("raster", "rdwd:::readDWD.grib2 with toraster=TRUE") out <- raster::raster(out) } +# WRONG --- +} else +tstop("pack='",pack,"' is not a valid option.") + # Output: return(invisible(out)) } diff --git a/man/readDWD.grib2.Rd b/man/readDWD.grib2.Rd index e065ee9..cb1b9ec 100644 --- a/man/readDWD.grib2.Rd +++ b/man/readDWD.grib2.Rd @@ -4,23 +4,38 @@ \alias{readDWD.grib2} \title{read nwp forecast data} \usage{ -readDWD.grib2(file, bargs = NULL, toraster = TRUE, quiet = rdwdquiet(), ...) +readDWD.grib2( + file, + pack = "terra", + bargs = NULL, + toraster = TRUE, + quiet = rdwdquiet(), + ... +) } \arguments{ \item{file}{Name of file on harddrive, like e.g. cosmo-d2_germany_regular-lat-lon_single-level_2021010100_005_T_2M.grib2.bz2} +\item{pack}{Char: package used for reading. +One of "terra" (the default), "stars" +or "rgdal" (for the deprecated cosmo-d2 data). +See \href{https://github.com/brry/rdwd/issues/28}{issue}. +DEFAULT: "terra"} + \item{bargs}{Named list of arguments passed to \code{\link[R.utils:compressFile]{R.utils::bunzip2()}}, see \code{gargs} in \code{\link[=readDWD.raster]{readDWD.raster()}}. DEFAULT: NULL} \item{toraster}{Logical: convert \link[rgdal:readGDAL]{rgdal::readGDAL} output with \code{\link[raster:raster]{raster::raster()}}? +Only used if pack="rgdal". DEFAULT: TRUE} \item{quiet}{Silence readGDAL completely, including warnings on discarded ellps / datum. DEFAULT: FALSE through \code{\link[=rdwdquiet]{rdwdquiet()}}} -\item{\dots}{Further arguments passed to \code{\link[rgdal:readGDAL]{rgdal::readGDAL()}},} +\item{\dots}{Further arguments passed to \code{\link[stars:read_stars]{stars::read_stars()}}, +\code{\link[rgdal:readGDAL]{rgdal::readGDAL()}} or \code{\link[rgdal:readGDAL]{rgdal::readGDAL()}}.} } \value{ rgdal or raster object, depending on \code{toraster} @@ -31,19 +46,24 @@ Intended to be called via \code{\link[=readDWD]{readDWD()}}.\cr } \examples{ \dontrun{ # Excluded from CRAN checks, but run in localtests -# Deactivated 2021-04-08 since readDWD.grib2 -> rgdal::readGDAL -> Error: -# **.grib2 is a grib file, but no raster dataset was successfully identified. -warning("readDWD.grib2 is not tested due to unresolved problems.") -if(FALSE){ -nwp_t2m_base <- "ftp://opendata.dwd.de/weather/nwp/icon-d2/grib/03/p" +nwp_t2m_base <- "ftp://opendata.dwd.de/weather/nwp/icon-d2/grib/15/soiltyp" nwp_urls <- indexFTP("", base=nwp_t2m_base, dir=tempdir()) -nwp_file <- dataDWD(nwp_urls[6], base=nwp_t2m_base, dir=tempdir(), +# for p instead of soiltyp, icosahedral_model-level files fail with GDAL errors, +# see https://github.com/brry/rdwd/issues/28 +# regular-lat-lon_pressure-level files work with pack="terra" or "stars" + +nwp_file <- dataDWD(tail(nwp_urls,1), base=nwp_t2m_base, dir=tempdir(), joinbf=TRUE, dbin=TRUE, read=FALSE) -nwp_data <- readDWD(nwp_file, quiet=TRUE) -plotRadar(nwp_data, project=FALSE) +nwp_data <- readDWD(nwp_file) +terra::plot(nwp_data) # same map with sp::plot +addBorders() # the projection seems to be perfectly good :) -nwp_data_rgdal <- readDWD(nwp_file, toraster=FALSE) -sp::plot(nwp_data_rgdal) +# index of GRIB files +if(FALSE){ # indexing takes about 6 minutes! +grib_base <- "ftp://opendata.dwd.de/weather/nwp/icon-d2/grib" +grib_files <- indexFTP("", base=grib_base, dir=tempdir()) +for(f in unique(substr(grib_files, 1,3))) print(grib_files[which(substr(grib_files, 1,3)==f)[1]]) +View(data.frame(grep("regular",grib_files, value=TRUE))) } } }