noaa.qmd

---
title: "Weather"
---

***DRAFT***

As a reference case, we can compare the NOAA [Global Ensemble Forecast System (GEFS)](https://www.ncei.noaa.gov/products/weather-climate-models/global-ensemble-forecast) predictions to measurements made on the ground at individual NEON sites using the EFI standard and cyberinfrastructure tools. Keep in mind these are not the usual NOAA forecasts you see in daily life, GEFS are: 31 ensemble models at 1.0 degree lat/long squares resolution that extend a full 30 days in the future. We sample these at the NEON sites and compare to NEON measurements.

```{r, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message=FALSE, warning=FALSE)
knitr::opts_chunk$set(eval=FALSE)
```

```{r setup}
library(ggiraph)
library(patchwork)
library(tidyverse)
library(lubridate)
library(neon4cast)
library(score4cast)
library(arrow)
library(glue)
library(thematic)
thematic_rmd(bg="white", fg="black", accent="blue")
source("R/plot-utils.R")
dashboard <- yaml::read_yaml("dashboard.yml") # Some configurable options

options(arrow.pull_as_vector = TRUE)

```



```{r}
neon_bucket <- function(table, endpoint = "https://sdsc.osn.xsede.org") {
  tbl =  stringr::str_replace(table, "-(DP\\d)", "/\\1") |> stringr::str_split_1("/")
  path = file.path("neon4cast-neonstore", tbl[[2]], tbl[[1]])
  bucket = paste0("bio230014-bucket01/", path)
  s3 <- arrow::S3FileSystem$create(endpoint_override = endpoint,
                                   access_key = Sys.getenv("OSN_KEY"),
                                   secret_key = Sys.getenv("OSN_SECRET"))
  s3_dir <- arrow::SubTreeFileSystem$create(bucket, s3)
  s3_dir
}

#library(neonstore)
#db <- neon_db()
#tables <- db |> DBI::dbListTables()


```


```{r}
# taat has no data after feb 28??
# FIXME replace with NEON-direct method

# temperature target
taat <- neon_bucket("TAAT_30min-basic-DP1.00003.001")

cutoff <- as.character(Sys.Date() - 70L)

sites <- arrow::open_dataset(taat) |> distinct(siteID) |> pull(siteID)
sites_subset <- sites[1:6]

temp <- arrow::open_dataset(taat) |> 
  mutate(startDateTime = as.character(startDateTime)) |>
  filter(startDateTime > cutoff, siteID %in% sites_subset) |>
  select(datetime = endDateTime, 
         observation = tempTripleMean, site_id = siteID) |>
  mutate(variable = "TMP", datetime = as_datetime(datetime)) |>
  collect()

```


```{r}

rh_bucket <- neon_bucket("RH_30min-basic-DP1.00098.001")
rh <- arrow::open_dataset(rh_bucket) |> 
  mutate(startDateTime = as.character(startDateTime)) |>
  filter(startDateTime > cutoff, siteID %in% sites_subset) |>
  select(datetime = endDateTime, 
         observation = RHMean, site_id = siteID) |>
  mutate(variable = "RH", datetime = as_datetime(datetime)) |>
  collect()
```


## Most recent forecasts

Forecasts are generated daily (actually 4 times daily in GEFS), while much NEON data is released only monthly. We obtain the date for the most recent observation (typically, the last day of the previous month.)  The most recent forecast for which we would now have observations to compare against the full 35 day forecast horizon would thus be the one made 35 days before this.  Thus the NEON latency and the forecast horizon mean that the most recent fully scored forecast will tend to be one made somewhere between 1 and 2 months ago:

```{r}
# reference_datetime for a fully-scored forecast
ref <- temp |> summarise(ref = max(datetime))  |> pull(ref) |> lubridate::as_date() -35L
ref
```

Lets access that forecast from the NOAA snapshots.


```{r}
s3 <- S3FileSystem$create(endpoint_override = "data.ecoforecast.org")
path <- SubTreeFileSystem$create(glue("neon4cast-drivers/noaa/gefs-v12/stage1/0/{ref}"), s3)
noaa_ds <- open_dataset(path) 
```



```{r}
noaa <- noaa_ds |> 
  filter(site_id %in% sites_subset) |>
  mutate(datetime = as_datetime(datetime),
         model_id = "noaa_gefs")
```







::: panel-tabset
## Temperature

```{r}
tmp_scores <- noaa |> 
  filter(variable == "TMP") |> 
  collect() |>
  score4cast::score(temp) |> 
  filter(!is.na(observation))

tmp_scores |> filter(site_id %in% sites_subset)  |> forecast_plots()
```

## RH

```{r}
rh_scores <- noaa |> 
  filter(variable == "RH") |> 
  collect() |>
  score4cast::score(rh) |> 
  filter(!is.na(observation)) 


rh_scores  |> filter(site_id %in% sites_subset)  |> forecast_plots()
```


:::

## Skill scores

::: panel-tabset

## Temperature

```{r}
tmp_scores |> ggplot(aes(datetime, crps)) + geom_line() + facet_wrap(~site_id)
```

```{r}

tmp_scores |> group_by(datetime) |> 
  summarise(crps = mean(crps, na.rm=TRUE),
            logs = mean(logs), na.rm=TRUE) |>
  pivot_longer(c("crps", "logs"), names_to="metric", values_to = "score") |>
  ggplot(aes(datetime, score)) + geom_line() + facet_wrap(~metric, scales="free")

```
```{r}
tmp_scores |> 
  group_by(site_id) |>
summarise(crps = mean(crps, na.rm=TRUE),
            logs = mean(logs), na.rm=TRUE) |>
  pivot_longer(c("crps", "logs"), 
               names_to="metric", values_to = "score") |>
  ggplot(aes(site_id, score, fill=metric)) + 
  geom_col(position="dodge") + facet_wrap(~metric, scales="free")

```


## Humidity

```{r}
rh_scores |> ggplot(aes(datetime, crps)) + geom_line() + facet_wrap(~site_id)
```

```{r}

rh_scores |> group_by(datetime) |> 
  summarise(crps = mean(crps, na.rm=TRUE),
            logs = mean(logs), na.rm=TRUE) |>
  pivot_longer(c("crps", "logs"), names_to="metric", values_to = "score") |>
  ggplot(aes(datetime, score)) + geom_line() + facet_wrap(~metric, scales="free")

```
```{r}
rh_scores |> 
  group_by(site_id) |>
summarise(crps = mean(crps, na.rm=TRUE),
            logs = mean(logs), na.rm=TRUE) |>
  pivot_longer(c("crps", "logs"), 
               names_to="metric", values_to = "score") |>
  ggplot(aes(site_id, score, fill=metric)) + 
  geom_col(position="dodge") + facet_wrap(~metric, scales="free")

```
:::