---
title: "Example Usage"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Example Usage}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
  
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```

```{r message = FALSE}
library(ggplot2)
library(dplyr)
library(purrr)
library(tibble)
library(ropenmeteo)
library(lubridate)
library(readr)
```

## Weather forecasts

The open-meteo project combines the the best models for each location across the globe to provide the best possible forecast.  open-meteo defines this as `model = "generic"`.

[https://open-meteo.com/en/docs]

```{r}
df <- get_forecast(latitude = 37.30,
                   longitude = -79.83,
                   forecast_days = 7, 
                   past_days = 2, 
                   model = "generic",
                   variables = c("temperature_2m"))
head(df)
```

```{r}
df |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction)) + 
  geom_line(color = "#F8766D") + 
  geom_vline(aes(xintercept = reference_datetime)) + 
  facet_wrap(~variable, scale = "free")
```

## Ensemble Weather Forecasts

Ensemble forecasts from individual models are available. 

[https://open-meteo.com/en/docs/ensemble-api]

```{r eval = FALSE}
df <- get_ensemble_forecast(
  latitude = 37.30,
  longitude = -79.83,
  forecast_days = 7,
  past_days = 2,
  model = "gfs_seamless",
  variables = c("temperature_2m"))
head(df)
```

```{r forecast-plot, eval=FALSE}
df |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction, color = ensemble)) + 
  geom_line() + 
  geom_vline(aes(xintercept = reference_datetime)) + 
  facet_wrap(~variable, scale = "free", ncol = 2)
```

Options for models and variables are at https://open-meteo.com/en/docs/ensemble-api

Note that `ecmwf_ifs04` does not include solar radiation.  

List of global model ids: 

```
icon_seamless, icon_global, gfs_seamless, gfs025, gfs05, ecmwf_ifs04, gem_global
```

### Use with the General Lake Model

We have included functions that allow the output to be used with the General Lake Model ([https://doi.org/10.5194/gmd-12-473-2019]).
Since the open-meteo models do not include longwave radiation, the package provides a function to calculate it from the cloud cover and air temperature.

GLM requires a set of variables that are provided

```{r}
df <- get_ensemble_forecast(
  latitude = 37.30,
  longitude = -79.83,
  forecast_days = 7,
  past_days = 2,
  model = "gfs_seamless",
  variables = glm_variables(product = "ensemble_forecast", 
                                        time_step = "hourly"))
head(df)
```

```{r emsemble-plot}
df |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction, color = ensemble)) + 
  geom_line() + 
  geom_vline(aes(xintercept = reference_datetime)) + 
  facet_wrap(~variable, scale = "free", ncol = 2)
```

The following converts to GLM format

```{r}
path <- tempdir()
df |> 
  add_longwave() |>
  write_glm_format(path = path)
head(read_csv(list.files(path = path, full.names = TRUE, pattern = ".csv")[1]))
```

### Converting to Ecological Forecasting Initative convention

The standard used in the NEON Ecological Forecasting Challenge is slightly different from the standard in this package.  It uses the column `parameter` for ensemble because the Challenge standard allows the flexibility to use parametric distributions (i.e., normal distribution `mean` and `sd`) in the same standard as a ensemble (or sample) forecast.  The `family` column defines the distribution (here `family` = `ensemble`).

The EFI standard also follows CF-conventions so the variable names are converted to be CF compliant. 

The output from `convert_to_efi_standard()` is the same as the output from `neon4cast::stage2()`

Learn more about `neon4cast::stage2()` here: [https://projects.ecoforecast.org/neon4cast-docs/Shared-Forecast-Drivers.html]

```{r}
df |>
  add_longwave() |>
  convert_to_efi_standard()
```

Note that `neon4cast::stage3()` is similar to

```{r}
df |>
  add_longwave() |>
  convert_to_efi_standard() |> 
  filter(datetime < reference_datetime)
```

With the only difference that the number of days is equal to the `past_days` in the call to `get_ensemble_forecast()`.  The max `past_days` from open-meteo is ~60 days.

## Historical Weather

If you need more historical days for model calibration and testing, historical data are available through open-meteo's historical weather API.

[https://open-meteo.com/en/docs/historical-weather-api] 

```{r}
df <- get_historical_weather(
  latitude = 37.30,
  longitude = -79.83,
  start_date = "2023-01-01",
  end_date = Sys.Date() - lubridate::days(1),
  variables = c("temperature_2m")) 
tail(df |> na.omit())
```

Notice the delay of ~7 days. 

```{r hist-plot}
df |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction)) + 
  geom_line(color = "#F8766D") + 
  geom_vline(aes(xintercept = lubridate::with_tz(Sys.time(), tzone = "UTC"))) + 
  facet_wrap(~variable, scale = "free")
```

## Seasonal Forecasts

Weather forecasts for up to 9 months in the future are available from the NOAA Climate Forecasting System

[https://open-meteo.com/en/docs/seasonal-forecast-api]

```{r}
df <- get_seasonal_forecast(
  latitude = 37.30,
  longitude = -79.83,
  forecast_days = 274,
  past_days = 5,
  variables = c("temperature_2m"))
head(df)
```

```{r seasonal-plot}
df |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction, color = ensemble)) + 
  geom_line() + 
  geom_vline(aes(xintercept = reference_datetime)) +
  facet_wrap(~variable, scale = "free")
```

### Downscaling from 6 hour to 1 hour time-step

The downscaling uses the GLM variables

```{r}
df <- get_seasonal_forecast(
  latitude = 37.30,
  longitude = -79.83,
  forecast_days = 30,
  past_days = 5,
  variables = glm_variables(product = "seasonal_forecast", 
                            time_step = "6hourly"))
```


```{r downscale-plot}
df |> 
  six_hourly_to_hourly(latitude = 37.30, longitude = -79.83, use_solar_geom = TRUE) |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction, color = ensemble)) + 
  geom_line() + 
  geom_vline(aes(xintercept = reference_datetime)) + 
  facet_wrap(~variable, scale = "free", ncol = 2)
```


## Climate Projections

Climate projections from different models are available through 2050.  The output is a daily time-step.

Note the units for shortwave radiation are different for the climate projection.

[https://open-meteo.com/en/docs/climate-api]

```{r}
df <- get_climate_projections(
  latitude = 37.30,
  longitude = -79.83,
  start_date = Sys.Date(),
  end_date = Sys.Date() + lubridate::years(1),
  model = "EC_Earth3P_HR",
  variables = c("temperature_2m_mean"))
head(df)
```

```{r climate-plot}
df |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction)) + 
  geom_line(color = "#F8766D") + 
  facet_wrap(~variable, scale = "free")
```

## Downloading multiple sites or models

### Multiple models

```{r}
models <- c("CMCC_CM2_VHR4","FGOALS_f3_H","HiRAM_SIT_HR","MRI_AGCM3_2_S","EC_Earth3P_HR","MPI_ESM1_2_XR","NICAM16_8S")

df <- map_df(models, function(model){
  get_climate_projections(
    latitude = 37.30,
    longitude = -79.83,
    start_date = Sys.Date(),
    end_date = Sys.Date() + lubridate::years(1),
    model = model,
    variables = c("temperature_2m_mean"))
})

```

```{r multi-climate-plot}
df |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction, color = model_id)) + 
  geom_line() +
  facet_wrap(~variable, scale = "free")
```

### Multiple sites

The download of multiple sites uses the optional `site_id` to add column that denotes the different sites.

```{r}
sites <- tibble(site_id = c("fcre", "sunp"),
                latitude = c(37.30, 43.39),
                longitude = c(-79.83, -72.05))

df <- map_df(1:nrow(sites), function(i, sites){
  get_climate_projections(
    latitude = sites$latitude[i],
    longitude = sites$longitude[i],
    site_id = sites$site_id[i],
    start_date = Sys.Date(),
    end_date = Sys.Date() + lubridate::years(1),
    model = "MPI_ESM1_2_XR",
    variables = c("temperature_2m_mean"))
},
sites)
head(df)
```

```{r multi-site-plot}
df |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction, color = site_id)) + 
  geom_line() +
  facet_wrap(~variable, scale = "free")
```

### Converting from daily to hourly time-step

Photosynthesis is non-linearly sensitive to shortwave radiation.  Therefore, the photosynthesis response to hourly radiation is different than the response to the aggregated daily mean radiation.  To address this issue, we provide a function to convert the daily sum of shortwave radiation to hourly values that uses solar geometry to impute.  Additionally, the sum of precipitation is divided by 24 hours to convert to an hourly time-step.  All other variables have their daily mean applied to each hour.

```{r}
df <- get_climate_projections(
  latitude = 37.30,
  longitude = -79.83,
  start_date = Sys.Date(),
  end_date = Sys.Date() + lubridate::years(1),
  model = "EC_Earth3P_HR",
  variables = glm_variables(product = "climate_projection", time_step = "daily"))
```

```{r daily-to-hourly-plot}
df |> 
  daily_to_hourly(latitude = 37.30, longitude = -79.83) |> 
  mutate(variable = paste(variable, unit)) |> 
  ggplot(aes(x = datetime, y = prediction)) + 
  geom_line(color = "#F8766D") + 
  facet_wrap(~variable, scale = "free", ncol = 2)
```
