ERA-5: add more metadata (#4), work on time series (#5)

.github/workflows/pkgdown.yaml

@@ -46,4 +46,4 @@ jobs:
         run: |
           git config --local user.email "[email protected]"
           git config --local user.name "GitHub Actions"
-          Rscript -e 'pkgdown::deploy_to_branch(new_process = FALSE)'
+          Rscript -e 'Rscript -e 'kwb.pkgbuild::deploy_to_branch_with_extra_files(vignettes_file_pattern_to_copy = "\\\.json$|\\\.html$|\\\\\.gif$")''

DESCRIPTION

@@ -33,9 +33,13 @@ Imports:
     tidyr
 Suggests: 
     covr,
+    cubelyr,
     knitr,
+    kwb.pkgbuild,
     raster,
-    rmarkdown
+    rmarkdown,
+    stars
 Remotes: 
-    github::kwb-r/kwb.utils
+    github::kwb-r/kwb.utils,
+    github::kwb-r/kwb.pkgbuild
 VignetteBuilder: knitr

NAMESPACE

@@ -15,4 +15,5 @@ importFrom(parallel,detectCores)
 importFrom(parallel,makeCluster)
 importFrom(parallel,parLapply)
 importFrom(parallel,stopCluster)
+importFrom(stringr,str_remove)
 importFrom(tidyr,pivot_wider)

R/.time-series.R

@@ -0,0 +1,55 @@
+library(kwb.satellite)
+
+grib_file <- "reanalysis-era5-single-levels_total_precipitation_2010-2020.grib"
+
+meta <- kwb.satellite::get_metadata_era5(grib_file)
+
+preci_stars <- stars::read_stars(grib_file)
+names(preci_stars) <- "era5_grib"
+
+
+preci_stars_cube <- cubelyr::as.tbl_cube(preci_stars)
+str(preci_stars_cube$dims)
+
+nbands <- length(preci_stars_cube$dims$band)
+has_data <- sapply(seq_len(nbands), function(band) {sum(preci_stars_cube$mets[[1]][,,band])>0})
+max_data <- sapply(seq_len(nbands), function(band) {max(preci_stars_cube$mets[[1]][,,band])})
+
+bands <- which(has_data)[1:100]
+
+val_max <- ceiling(max(max_data)*1000)
+
+max_col <- 5
+
+cols <- data.frame(breaks = c(0, seq_len(max_col)*val_max/max_col),
+                   col_name = heat.colors(max_col))
+
+animation::saveGIF(expr = for(band in bands) {
+label <- sprintf("%s [mm] (%s UTC)", meta$variable[band], meta$time_forecast[band])
+preci_stars %>%
+  dplyr::slice(index = band, along = "band") %>%
+  ### Convert units from [m] to [mm]
+  dplyr::mutate(era5_grib = era5_grib * 1000) %>%
+  raster::plot(zlim=c(0,max(max_data)*1000), main = label)
+},
+movie.name = "era5-precipitation.gif",
+interval = 0.1
+)
+
+preci_stars_points <- sf::st_as_sf(preci_stars, as_points = TRUE, merge = FALSE)
+
+
+#summary <- rgdal::GDALinfo(grib_file, silent = TRUE)
+
+# preci_sf <- sf::st_read(grib_file)
+# preci <- rgdal::readGDAL(grib_file)
+# head(preci@data)
+# class(preci)
+# raster::image(preci, zlim = c(0, 0.01))
+# ggplot2::ggplot(preci)
+# image(preci)
+# preci_df <- as(preci, "SpatialGridDataFrame")
+#
+# preci@bbox
+#
+# plot(preci$band1)

R/get_metadata_era5.R

@@ -8,6 +8,7 @@
 #' @importFrom tidyr pivot_wider
 #' @importFrom gdata startsWith
 #' @importFrom dplyr case_when mutate
+#' @importFrom stringr str_remove
 #' @seealso Code taken from https://gis.stackexchange.com/a/360652
 #'
 get_metadata_era5 <- function(grib_file) {
@@ -22,6 +23,7 @@ get_metadata_era5 <- function(grib_file) {
                              "GRIB_ELEMENT",
                              "GRIB_FORECAST_SECONDS",
                              "GRIB_REF_TIME",
+                             "GRIB_VALID_TIME",
                              "GRIB_SHORT_NAME",
                              "GRIB_UNIT"),
                            collapse = "|")
@@ -44,6 +46,7 @@ get_metadata_era5 <- function(grib_file) {
   is_grib_e <- raw_starts_with('GRIB_ELEMENT')
   is_grib_f <- raw_starts_with('GRIB_FORECAST_SECONDS')
   is_grib_r <- raw_starts_with('GRIB_REF_TIME')
+  is_grib_v <- raw_starts_with('GRIB_VALID_TIME')
   is_grib_s <- raw_starts_with('GRIB_SHORT_NAME')
   is_grib_u <- raw_starts_with('GRIB_UNIT')
 
@@ -55,17 +58,19 @@ get_metadata_era5 <- function(grib_file) {
         is_grib_e ~ sub(".*=", "", raw_),
         is_grib_f ~ gsub(".*=(.+) sec.*", "\\1", raw_),
         is_grib_r ~ gsub(".*= (.+) sec.*", "\\1", raw_),
+        is_grib_v ~ gsub(".*= (.+) sec.*", "\\1", raw_),
         is_grib_s ~ sub(".*=", "", raw_),
         is_grib_u ~ sub(".*=", "", raw_) %>%
                     stringr::str_remove("\\[") %>%
                     stringr::str_remove("\\]")
       ),
       column = dplyr::case_when(
         is_band ~ 'band',
-        is_grib_c ~ 'variable',
+        is_grib_c ~ 'variable_unit',
         is_grib_e ~ 'variable_short',
         is_grib_f ~ 'forecast_seconds',
-        is_grib_r ~ 'time',
+        is_grib_r ~ 'time_ref',
+        is_grib_v ~ 'time_valid',
         is_grib_s ~ 'short_name',
         is_grib_u ~ 'unit'
       )
@@ -86,8 +91,11 @@ get_metadata_era5 <- function(grib_file) {
     values_from = "content"
   ) %>%
     dplyr::mutate(
-      time = as.POSIXct(as.numeric(time), origin = "1970-01-01", tz = "UTC")
-    )
+      time_ref = as.POSIXct(as.numeric(time_ref), origin = "1970-01-01", tz = "UTC"),
+      time_valid = as.POSIXct(as.numeric(time_valid), origin = "1970-01-01", tz = "UTC"),
+      time_forecast = time_ref + as.numeric(forecast_seconds),
+      variable = stringr::str_remove(variable_unit,
+                                     pattern = "\\s+\\[.*\\]$"))
 
   grib1 <- grib1[, -1]