convert Rmd to R script

Submit_forecast/forecast_code_template.Rmd → Submit_forecast/forecast_code_template.R

@@ -1,89 +1,70 @@
----
-title: "NEON forecast challenge submission"
-output: html_document
-date: "`r Sys.Date()`"
----
-
-```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE)
-```
-
-
-```{r load-packages, echo = F, warning=F, message=F}
 ## install.packages('remotes')
-## install.packages('fpp3') # package for applying simple forecasting methods
-## install.packages('tsibble') # package for dealing with time series data sets and tsibble objects
 ## install.packages('tidyverse') # collection of R packages for data manipulation, analysis, and visualisation
 ## install.packages('lubridate') # working with dates and times
 ## remotes::install_github('eco4cast/neon4cast') # package from NEON4cast challenge organisers to assist with forecast building and submission
 
-# Load packages
+# ------ Load packages -----
 library(tidyverse)
 library(lubridate)
+#--------------------------#
 
-```
 
-```{r get-targets, message=F}
-#read in the targets data
+
+#------- Read data --------
+# read in the targets data
 targets <- read_csv('https://data.ecoforecast.org/neon4cast-targets/aquatics/aquatics-targets.csv.gz')
 
 # read in the sites data
 aquatic_sites <- read_csv("https://raw.githubusercontent.com/eco4cast/neon4cast-targets/main/NEON_Field_Site_Metadata_20220412.csv") |>
   dplyr::filter(aquatics == 1)
 
-lake_sites <- aquatic_sites %>%
-  filter(field_site_subtype == 'Lake')
+focal_sites <- aquatic_sites |> 
+  filter(field_site_subtype == 'Lake') |> 
+  pull(field_site_id)
 
 # Filter the targets
 targets <- targets %>%
-  filter(site_id %in% lake_sites$field_site_id,
+  filter(site_id %in% focal_sites,
          variable == 'temperature')
-```
+#--------------------------#
 
-This is where you could change or add meteorological variables that are used to predict the target
 
-Other variable names can be found at <https://projects.ecoforecast.org/neon4cast-docs/Shared-Forecast-Drivers.html#stage-3>
 
-```{r}
+# ------ Weather data ------
 met_variables <- c("air_temperature")
-```
-
-```{r get-NOAA-past, message = F}
 
-# Past stacked weather
-noaa_past_s3 <- neon4cast::noaa_stage3()
+# Past stacked weather -----
+weather_past_s3 <- neon4cast::noaa_stage3()
 
-noaa_past <- noaa_past_s3  |> 
-  dplyr::filter(site_id %in% lake_sites$field_site_id,
+weather_past <- weather_past_s3  |> 
+  dplyr::filter(site_id %in% focal_sites,
                 datetime >= ymd('2017-01-01'),
                 variable %in% met_variables) |> 
   dplyr::collect()
 
 # aggregate the past to mean values
-noaa_past_mean <- noaa_past |> 
+weather_past_daily <- weather_past |> 
   mutate(datetime = as_date(datetime)) |> 
   group_by(datetime, site_id, variable) |> 
   summarize(prediction = mean(prediction, na.rm = TRUE), .groups = "drop") |> 
   # convert air temperature to Celsius if it is included in the weather data
   mutate(prediction = ifelse(variable == "air_temperature", prediction - 273.15, prediction)) |> 
-  pivot_wider(names_from = variable, values_from = prediction) 
-```
+  pivot_wider(names_from = variable, values_from = prediction)
 
-```{r get-NOAA-future, message = F}
-# Future weather
+# Future weather forecast --------
 # New forecast only available at 5am UTC the next day
 forecast_date <- Sys.Date() 
 noaa_date <- forecast_date - days(1)
 
-noaa_future_s3 <- neon4cast::noaa_stage2(start_date = as.character(noaa_date))
+weather_future_s3 <- neon4cast::noaa_stage2(start_date = as.character(noaa_date))
 
-noaa_future <- noaa_future_s3 |> 
+weather_future <- weather_future_s3 |> 
   dplyr::filter(datetime >= forecast_date,
-                site_id %in% lake_sites$field_site_id,
+                site_id %in% focal_sites,
                 variable %in% met_variables) |> 
   collect()
 
-noaa_future_daily <- noaa_future |> 
+weather_future_daily <- weather_future |> 
   mutate(datetime = as_date(datetime)) |> 
   # mean daily forecasts at each site per ensemble
   group_by(datetime, site_id, parameter, variable) |> 
@@ -92,9 +73,13 @@ noaa_future_daily <- noaa_future |>
   mutate(prediction = ifelse(variable == "air_temperature", prediction - 273.15, prediction)) |> 
   pivot_wider(names_from = variable, values_from = prediction) |> 
   select(any_of(c('datetime', 'site_id', met_variables, 'parameter')))
-```
 
-```{r model-setup}
+#--------------------------#
+
+
+
+# ----- Fit model & generate forecast----
+
 # Generate a dataframe to fit the model to 
 targets_lm <- targets |> 
   pivot_wider(names_from = 'variable', values_from = 'observation') |> 
@@ -104,94 +89,75 @@ targets_lm <- targets |>
 # Loop through each site to fit the model
 forecast_df <- NULL
 
-```
-
-```{r forecast-loop}
 for(i in 1:length(lake_sites$field_site_id)) {  
 
   example_site <- lake_sites$field_site_id[i]
 
   site_target <- targets_lm |>
     filter(site_id == example_site)
-
+  
   noaa_future_site <- noaa_future_daily |> 
     filter(site_id == example_site)
 
   #Fit linear model based on past data: water temperature = m * air temperature + b
   #you will need to change the variable on the left side of the ~ if you are forecasting oxygen or chla
   fit <- lm(site_target$temperature ~ site_target$air_temperature)
   # fit <- lm(site_target$temperature ~ ....)
-    
+
   # use linear regression to forecast water temperature for each ensemble member
   # You will need to modify this line of code if you add additional weather variables or change the form of the model
   # The model here needs to match the model used in the lm function above (or what model you used in the fit)
   forecasted_temperature <- fit$coefficients[1] + fit$coefficients[2] * noaa_future_site$air_temperature
-    
+
   # put all the relevant information into a tibble that we can bind together
   curr_site_df <- tibble(datetime = noaa_future_site$datetime,
-                        site_id = example_site,
-                        parameter = noaa_future_site$parameter,
-                        prediction = forecasted_temperature,
-                        variable = "temperature") #Change this if you are forecasting a different variable
+                         site_id = example_site,
+                         parameter = noaa_future_site$parameter,
+                         prediction = forecasted_temperature,
+                         variable = "temperature") #Change this if you are forecasting a different variable
 
   forecast_df <- dplyr::bind_rows(forecast_df, curr_site_df)
   message(example_site, 'forecast run')
 
 }
-```
+my_model_id <- 'example_ID'
 
-Remember to change the model_id when you make changes to the model structure!
+#--------------------------#
 
-```{r}
-my_model_id <- 'example_ID'
-```
 
-```{r make-standard}
+#---- Covert to EFI standard ----
+
 # Make forecast fit the EFI standards
 forecast_df_EFI <- forecast_df %>%
   filter(datetime > forecast_date) %>%
   mutate(model_id = my_model_id,
          reference_datetime = forecast_date,
          family = 'ensemble',
-         parameter = as.character(parameter)) %>%
-  select(datetime, reference_datetime, site_id, family, parameter, variable, prediction, model_id)
+         duration = 'P1D',
+         parameter = as.character(parameter),
+         project_id = 'neon4cast') %>%
+  select(datetime, reference_datetime, duration, site_id, family, parameter, variable, prediction, model_id, project_id)
+#---------------------------#
+
 
-```
 
-```{r write-forecast}
+# ----- Submit forecast -----
 # Write the forecast to file
 theme <- 'aquatics'
 date <- forecast_df_EFI$reference_datetime[1]
-forecast_name_1 <- paste0(forecast_df_EFI$model_id[1], ".csv")
-forecast_file_1 <- paste(theme, date, forecast_name_1, sep = '-')
-forecast_file_1
-
-
-if (!dir.exists('Forecasts')) {
-  dir.create('Forecasts')
-}
-
-write_csv(forecast_df_EFI, file.path('Forecasts',forecast_file_1))
-```
-
-Check that forecast format is valid
-
-```{r}
-neon4cast::forecast_output_validator(file.path('Forecasts',forecast_file_1))
-```
+forecast_name <- paste0(forecast_df_EFI$model_id[1], ".csv")
+forecast_file <- paste(theme, date, forecast_name, sep = '-')
 
-Change eval = TRUE if you want to submit
+write_csv(forecast_df_EFI, forecast_file)
 
-```{r submit-forecast, eval= FALSE}
+neon4cast::forecast_output_validator(forecast_file)
 
-neon4cast::submit(forecast_file = file.path('Forecasts', forecast_file_1), ask = FALSE) # if ask = T (default), it will produce a pop-up box asking if you want to submit
 
-```
+neon4cast::submit(forecast_file =  forecast_file, ask = FALSE) # if ask = T (default), it will produce a pop-up box asking if you want to submit
+#--------------------------#
 
-```{r plot-forecast}
 forecast_df_EFI |> 
   ggplot(aes(x=datetime, y=prediction, group = parameter)) +
   geom_line() +
   facet_wrap(~site_id) +
   labs(title = paste0('Forecast generated for ', forecast_df_EFI$variable[1], ' on ', forecast_df_EFI$reference_datetime[1]))
-```