Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added code to extract ensemble forecasts #25

Draft
wants to merge 5 commits into
base: main
Choose a base branch
from

Conversation

dlebauer
Copy link
Collaborator

@dlebauer dlebauer commented Dec 13, 2022

This code

  • downloads forecasts
  • downscales using interpolation
  • extracts 30 ensemble members of forecasts for each of the 28 azmet stations
  • combines these into a dataframe
  • creates long version and long version with summary stats
  • makes some exploratory plots

@dlebauer
Copy link
Collaborator Author

image

@Aariq
Copy link
Member

Aariq commented Dec 15, 2022

Thanks David! I'll review this and then either implement it as part of the targets workflow or possibly deploy it on RStudio Connect on its own and have it save the data as a pin that the targets workflow can use.

data(station_info, package = 'azmetr')

library(gdalcubes)
# remotes::install_github("neon4cast/gefs4cast")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This one doesn't seem to exist, so I assume it's the eco4cast/gefs4cast one

Suggested change
# remotes::install_github("neon4cast/gefs4cast")

```{r}
library(dplyr)
## set path for gdal binary
Sys.setenv("GDAL_BIN"="/usr/bin/")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My gdal was in usr/local/bin so I didn't have to set this. I think usr/local/bin is the default homebrew install location. Ironically, my problem was with my location of bash! This package assumes it's at /usr/bin/bash and mine is /bin/bash. I submitted an issue (eco4cast/gefs4cast#7) and I'll try again once it's addressed.

@Aariq
Copy link
Member

Aariq commented Dec 20, 2022

gefs4cast is buggy, can't get running for now. Opened issue: eco4cast/gefs4cast#7

Copy link
Member

@Aariq Aariq left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had a lot of issues with gdal and geos system libraries and eventually gave up getting it to work on my laptop and instead tried out rocker/geospatial. Everything except the animated plot works with that image, so if we do implement GEFS-based QA we can use that as a base image for deployment.


library(gdalcubes)
# remotes::install_github("neon4cast/gefs4cast")
#remotes::install_github('eco4cast/gefs4cast')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The renv way to do this:

Suggested change
#remotes::install_github('eco4cast/gefs4cast')
#renv::install('eco4cast/gefs4cast')

```{r}
library(dplyr)
## set path for gdal binary
Sys.setenv("GDAL_BIN"="/usr/bin/")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

New default in dev version of gefs4cast is to respect PATH so this probably shouldn't be set.

# may need to adjust for time zones?
# yesterday's forecast starts today?

gefsdir <- fs::dir_create(file.path("gefs", date))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gefs_cog() will create the dir.

Suggested change
gefsdir <- fs::dir_create(file.path("gefs", date))
gefsdir <- file.path("gefs", date)

Comment on lines +45 to +49
gefs_cog(gefsdir,
ens_avg = FALSE,
max_horizon = 48, # get next 48 h. might only need 24h?
date = date) |>
system.time()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

these are big files, even with compression. I think we would quickly blow through our GitHub Actions allocation with this and would need to figure out another way to deploy this.

```{r}
## create statistical summaries by station, variable, and time

ens_stats <- ens_long %>%
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is ens_long? Looks like a pivot_longer() step is missing. I think probably:

Suggested change
ens_stats <- ens_long %>%
library(tidyr)
ens_long <-
ens_forecast %>%
select(-geometry) %>%
pivot_longer(APCP:VGRD, names_to = "var")
ens_stats <- ens_long %>%


# above plot but saved as one file per station, using a for loop
plots <- list()
for(stn in unique(ens_stats$meta_station_id)){
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
for(stn in unique(ens_stats$meta_station_id)){
for(stn in unique(ens_stats$meta_station_name)){

ucl_99 = quantile(value, 0.995),
lcl_99 = quantile(value, 0.005)) %>%
ungroup() %>%
mutate(id = paste0(meta_station_name, var))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I needed to add a mutate to parse time to get the plots to work.

Suggested change
mutate(id = paste0(meta_station_name, var))
mutate(id = paste0(meta_station_name, var), time = lubridate::ymd_hms(time))

Comment on lines +178 to +181
raster_cube(gefs_cube, v) |>
gdalcubes::select_bands("TMP") |>
animate( col = viridisLite::viridis, nbreaks=100,
fps=10, save_as = "temp.gif")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Requires packge gifski which requires some system library that isn't in rocker/geospatial. Cool idea, but probably won't implement in the qa/qc report.

@Aariq Aariq closed this Jan 3, 2023
@Aariq Aariq reopened this Jan 3, 2023
@Aariq Aariq removed their assignment Jan 3, 2023
@Aariq Aariq marked this pull request as draft February 3, 2023 21:14
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants