added code to extract ensemble forecasts

for each station

notes/gefs_forecasts.qmd

@@ -0,0 +1,170 @@
+---
+title: "Ensemble forecasts"
+author: "David LeBauer"
+format: html
+     df-print: paged
+editor: source
+---
+
+
+## gefs4cast
+
+### gefs4cast::noaa_gefs
+
+
+```{r}
+library(dplyr)
+## set path for gdal binary 
+Sys.setenv("GDAL_BIN"="/usr/bin/")
+#remotes::install_github("cct-datascience/azmetr")
+#library(azmetr)
+data(station_info, package = 'azmetr')
+readr::write_csv(station_info %>% 
+  dplyr::mutate(site_name = meta_station_name, 
+         site_id = meta_station_id),
+  "new_station_info.csv")
+
+
+gefs4cast::noaa_gefs(date = Sys.Date(),
+  cycle = "00",
+  threads = 70,
+  gdal_ops = "",
+  s3 = arrow::s3_bucket("drivers", endpoint_override = "data.ecoforecast.org"),
+  max_horizon = 48,
+  purge = TRUE,
+  quiet = FALSE,
+  dest = tempdir(),#".",
+  locations = "new_station_info.csv",#paste0("https://github.com/eco4cast/neon4cast-noaa-download/",
+   # "raw/master/noaa_download_site_list.csv"),
+  name_pattern = "noaa/gefs-v12/stage1/{cycle_int}/{nice_date}/part-0.parquet")
+```
+
+
+### Vignette 
+
+From https://projects.ecoforecast.org/gefs4cast/articles/gefs-cog.html
+```{r}
+library(gdalcubes) 
+# remotes::install_github("neon4cast/gefs4cast")
+#remotes::install_github('eco4cast/gefs4cast')
+library(gefs4cast) 
+library(stringr)
+library(lubridate)
+library(stars)
+library(viridisLite)
+library(tmap)
+
+gdalcubes_options(parallel = TRUE)#24)
+```
+
+```{r}
+date <- Sys.Date() - lubridate::days(1)#as.Date("2022-12-05")
+
+gefsdir <- fs::dir_create("gefs")
+
+#gefs_cog(tmpdir, ens_avg=TRUE, max_horizon=240, date = date)
+
+gefs_cog(gefsdir, 
+         ens_avg = FALSE,
+         max_horizon = 48, 
+         date = date) |>
+  system.time()
+```
+
+```{r}
+# 000 file has different bands and so cannot be stacked into the collection
+
+
+gefs_cubes <- list()
+for (i in 1:30) {
+  ens <- stringr::str_pad(i, 2, pad = '0')
+  ens_id <- paste0('gep', ens)
+  .f <- fs::dir_ls(gefsdir, 
+                      type="file", 
+                      recurse = TRUE,
+                      regex = ens_id)[-1]
+  step <- .f |> str_extract("f\\d{3}\\.tif") |> str_extract("\\d{3}") 
+  datetime <- Sys.Date() + hours(step) 
+
+  fc <- stars::read_stars(.f[1]) 
+  bandnames <- 
+    stars::st_get_dimension_values(fc, 3) |>
+    stringr::str_extract("([A-Z]+):") |> str_remove(":")
+
+  gefs_cubes[[ens]] <- create_image_collection(.f,
+                                     date_time = datetime, 
+                                     band_names = bandnames)
+}
+```
+
+
+#### Cube View
+
+```{r}
+
+azmet_pts <- station_info %>% 
+  sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
+  
+azmet_box <- azmet_pts %>% 
+  sf::st_bbox() 
+
+iso <-"%Y-%m-%dT%H:%M:%S"
+ext <- list(t0 = as.character(min(datetime),iso), 
+              t1 = as.character(max(datetime),iso),
+              left = azmet_box[1], right = azmet_box[3],
+              top = azmet_box[4], bottom = azmet_box[2])
+v <- cube_view(srs = "EPSG:4326", 
+               extent = ext,
+              #dx = 0.5, dy = 0.5, # original resolution -- half-degree
+              dx = 0.1, 
+              dy = 0.1,
+              dt = "PT1H",
+              #dt = "PT3H",   # original resolution
+               aggregation = "mean", 
+               resampling = "cubicspline"
+)
+```
+
+```{r}
+
+preds <- list()
+for(ens in names(gefs_cubes)){
+  gefs_cube <- gefs_cubes[[ens]]
+  x <- raster_cube(gefs_cube, v)  |> 
+         extract_geom(sf = azmet_pts)
+  p <- azmet_pts
+  p$FID <- rownames(azmet_pts)
+  y <- merge(p, x, by = "FID")
+  preds[[ens]] <- y
+}
+ens_forecast <- dplyr::bind_rows(preds, .id = 'ens')
+```
+
+```{r}
+library(ggplot2)
+
+ggplot(data = ens_forecast) +
+  geom_line(aes(x = time, y = TMP, group = ens)) +
+  facet_wrap(~meta_station_name, ncol = 6) +
+  theme_bw()
+
+ggplot(data = ens_forecast) +
+  geom_line(aes(x = time, y = TMP, color = meta_station_name)) +
+  facet_wrap(~ens, ncol = 6) 
+
+  #geom_sf(aes(fill = TMP)) +
+  facet_wrap(~meta_station_name, ncol = 6) +
+  scale_fill_viridis_c(option = "plasma", na.value = "white") +
+  theme_void() +
+  theme(legend.position = "none")
+```
+
+```{r}
+x <- raster_cube(gefs_cube, v)  |> extract_geom(azmet_pts)
+
+raster_cube(gefs_cube, v) |>
+  gdalcubes::select_bands("TMP") |>
+  animate( col = viridisLite::viridis, nbreaks=100,
+           fps=10, save_as = "temp.gif")
+```
+