diff --git a/app/functions.R b/app/functions.R
new file mode 100644
index 0000000..a43a143
--- /dev/null
+++ b/app/functions.R
@@ -0,0 +1,80 @@
+get_prism_r <- function(dates, var){ # dates = dates_then
+ d <- d_prism_r |>
+ filter(
+ date %in% !!dates,
+ variable == !!var) |>
+ select(path_tif, lyr)
+
+ rast(unique(d$path_tif)) |>
+ subset(d$lyr) |>
+ mean() |>
+ project(leaflet:::epsg3857)
+}
+
+get_sst_r <- function(dates){ # dates = dates_then
+ lyrs <- d_sst_r |>
+ filter(
+ date %in% !!dates) |>
+ pull(lyr)
+
+ r_sst |>
+ subset(lyrs) |>
+ mean() |>
+ project(leaflet:::epsg3857)
+}
+
+map_then_now <- function(
+ r_then,
+ r_now,
+ tiles = providers$CartoDB.DarkMatter,
+ lgnd_then,
+ lgnd_now,
+ var_lbl,
+ palette = "Spectral",
+ palette_rev = TRUE){
+
+ vals <- c(values(r_now, na.rm=T), values(r_then, na.rm=T))
+ pal <- colorNumeric(
+ palette, vals, reverse = palette_rev, na.color = "transparent")
+
+ leaflet() |>
+ addMapPane("left", zIndex = 0) |>
+ addMapPane("right", zIndex = 0) |>
+ addProviderTiles(
+ tiles,
+ options = pathOptions(pane = "left"),
+ group = "base",
+ layerId = "base_l") |>
+ addProviderTiles(
+ tiles,
+ options = pathOptions(pane = "right"),
+ group = "base",
+ layerId = "base_r") |>
+ addRasterImage(
+ r_then, colors = pal, opacity = 0.8, project = F,
+ options = leafletOptions(pane = "left"),
+ group = "r_then") |>
+ addRasterImage(
+ r_now, colors = pal, opacity = 0.8, project = F,
+ options = leafletOptions(pane = "right"),
+ group = "r_now") |>
+ addLayersControl(overlayGroups = c("r_then", "r_now")) |>
+ addSidebyside(
+ layerId = "sidecontrols",
+ leftId = "base_l",
+ rightId = "base_r") |>
+ addControl(
+ HTML(lgnd_then),
+ position = "topleft") |>
+ addControl(
+ HTML(lgnd_now),
+ position = "topright") |>
+ addPolygons(
+ data = tbeptools::tbsegshed,
+ # label = ~long_name,
+ color="white", weight = 2, fillOpacity=0) |>
+ addLegend(
+ pal = pal,
+ values = vals,
+ title = var_lbl)
+}
diff --git a/app/global.R b/app/global.R
new file mode 100644
index 0000000..847ab7c
--- /dev/null
+++ b/app/global.R
@@ -0,0 +1,45 @@
+# TODO:
+# - [ ] style ggplotly https://rstudio.github.io/thematic/#shiny
+# - [ ] show sources for prism, sst, etc with About and links per card
+# - [ ] convert prism units: ºC -> ºF, mm -> in
+# - [ ] sel_t_var implement in map_prism_temp w/ avg
+# - [ ] swiper maps: precipitation, ocean temperature
+# - [ ] temp: heat index using dewpoint (`tdmean` -> humidity)
+# - [ ] prism: show `version`, `date_updated`
+# - [ ] swiper maps: leaflet-proxy-map updates to layers so can zoom / pan without refresh
+# - [ ] time series cards / panels
+
+# devtools::install_local(here("../tbeptools"), force = T)
+# devtools::load_all(here::here("../tbeptools"))
+librarian::shelf(
+ bsicons, bslib, dplyr, glue, here, leaflet, leaflet.extras2, lubridate, plotly,
+ readr, scales, shiny, slider,
+ tbep-tech/tbeptools,
+ terra, thematic, tidyr)
+source(here("app/functions.R"))
+options(readr.show_col_types = F)
+
+light <- bs_theme(preset = "flatly")
+dark <- bs_theme(preset = "darkly")
+
+dir_prism <- here("data/prism")
+prism_csv <- here("data/prism.csv")
+sst_tif <- here("data/sst/tb_sst.tif")
+sst_csv <- here("data/sst/tb_sst.csv")
+
+d_prism_r <- read_prism_rasters(dir_prism)
+yrs_prism <- range(year(d_prism_r$date))
+now_prism <- max(d_prism_r$date)
+d_prism_z <- read_csv(prism_csv)
+
+r_sst <- rast(sst_tif)
+d_sst_r <- tibble(
+ lyr = names(r_sst)) |>
+ separate(lyr, c("var", "date"), sep = "\\|", remove = F) |>
+ mutate(
+ date = as.Date(date))
+yrs_sst <- range(year(d_sst_r$date))
+now_sst <- max(d_sst_r$date)
+d_sst_z <- read_csv(sst_csv) |>
+ mutate(
+ date = as.Date(time))
diff --git a/app/server.R b/app/server.R
new file mode 100644
index 0000000..94144f6
--- /dev/null
+++ b/app/server.R
@@ -0,0 +1,230 @@
+function(input, output, session) {
+
+ # dark_mode ----
+ observe(session$setCurrentTheme(
+ if (isTRUE(input$dark_mode)) dark else light ))
+
+ # map_prism_temp ----
+ output$map_prism_temp <- renderLeaflet({
+
+ # DEBUG
+ # input <- list(
+ # sld_md = as.Date("2024-07-22"),
+ # sld_yrs_now = c(2024, 2024),
+ # sld_yrs_then = c(1981, 2001))
+
+ var = "tmax"
+ var_lbl = "Temperature (°C)"
+ md = sprintf("%02d-%02d", month(input$sld_t_md), day(input$sld_t_md))
+ yrs_now = input$sld_t_yrs_now[1]:input$sld_t_yrs_now[2]
+ yrs_then = input$sld_t_yrs_then[1]:input$sld_t_yrs_then[2]
+ dates_now = as.Date(glue("{yrs_now}-{md}"))
+ dates_then = as.Date(glue("{yrs_then}-{md}"))
+
+ if (any(dates_now > now_prism))
+ dates_now[dates_now > now_prism] <- dates_now[dates_now > now_prism] - years(1)
+
+ yrs_now_rng <- year(dates_now)
+ yrs_then_rng <- year(dates_then)
+ if (length(yrs_now_rng) > 2)
+ yrs_now_rng <- range(yrs_now_rng)
+ if (length(yrs_then_rng) > 2)
+ yrs_then_rng <- range(yrs_then_rng)
+
+ r_now <- get_prism_r(dates_now, var)
+ r_then <- get_prism_r(dates_then, var)
+
+ lgnd_now <- glue(
+ "Now
+ {format(input$sld_t_md, '%b %d')},
+ {paste(yrs_now_rng, collapse = ' to ')}")
+ lgnd_then <- glue(
+ "Then
+ {format(input$sld_t_md, '%b %d')},
+ {paste(yrs_then_rng, collapse = ' to ')}")
+
+ map_then_now(
+ r_then,
+ r_now,
+ tiles = ifelse(
+ isTRUE(input$dark_mode),
+ providers$CartoDB.DarkMatter,
+ providers$CartoDB.Positron),
+ lgnd_then,
+ lgnd_now,
+ var_lbl)
+
+ })
+
+ # map_prism_rain ----
+ output$map_prism_ppt <- renderLeaflet({
+
+ # DEBUG
+ # input <- list(
+ # sld_md = as.Date("2024-07-22"),
+ # sld_yrs_now = c(2024, 2024),
+ # sld_yrs_then = c(1981, 2001))
+
+ var = "ppt"
+ var_lbl = "Rain (mm)"
+ md = sprintf("%02d-%02d", month(input$sld_p_md), day(input$sld_p_md))
+ yrs_now = input$sld_p_yrs_now[1]:input$sld_p_yrs_now[2]
+ yrs_then = input$sld_p_yrs_then[1]:input$sld_p_yrs_then[2]
+ dates_now = as.Date(glue("{yrs_now}-{md}"))
+ dates_then = as.Date(glue("{yrs_then}-{md}"))
+
+ if (any(dates_now > now_prism))
+ dates_now[dates_now > now_prism] <- dates_now[dates_now > now_prism] - years(1)
+
+ yrs_now_rng <- year(dates_now)
+ yrs_then_rng <- year(dates_then)
+ if (length(yrs_now_rng) > 2)
+ yrs_now_rng <- range(yrs_now_rng)
+ if (length(yrs_then_rng) > 2)
+ yrs_then_rng <- range(yrs_then_rng)
+
+ r_now <- get_prism_r(dates_now, var)
+ r_then <- get_prism_r(dates_then, var)
+
+ lgnd_now <- glue(
+ "Now
+ {format(input$sld_p_md, '%b %d')},
+ {paste(yrs_now_rng, collapse = ' to ')}")
+ lgnd_then <- glue(
+ "Then
+ {format(input$sld_p_md, '%b %d')},
+ {paste(yrs_then_rng, collapse = ' to ')}")
+
+ map_then_now(
+ r_then,
+ r_now,
+ tiles = ifelse(
+ isTRUE(input$dark_mode),
+ providers$CartoDB.DarkMatter,
+ providers$CartoDB.Positron),
+ lgnd_then,
+ lgnd_now,
+ var_lbl,
+ palette = "Blues",
+ palette_rev = F)
+
+ })
+
+ # map_sst ----
+ output$map_sst <- renderLeaflet({
+
+ # DEBUG
+ # input <- list(
+ # sld_o_md = as.Date("2024-07-22"),
+ # sld_o_yrs_now = c(2024, 2024),
+ # sld_o_yrs_then = c(1981, 2001))
+
+ var_lbl = "Sea Surface Temperature (°C)"
+ md = sprintf("%02d-%02d", month(input$sld_o_md), day(input$sld_o_md))
+ yrs_now = input$sld_o_yrs_now[1]:input$sld_o_yrs_now[2]
+ yrs_then = input$sld_o_yrs_then[1]:input$sld_o_yrs_then[2]
+ dates_now = as.Date(glue("{yrs_now}-{md}"))
+ dates_then = as.Date(glue("{yrs_then}-{md}"))
+
+ if (any(dates_now > now_sst))
+ dates_now[dates_now > now_sst] <- dates_now[dates_now > now_sst] - years(1)
+
+ yrs_now_rng <- year(dates_now)
+ yrs_then_rng <- year(dates_then)
+ if (length(yrs_now_rng) > 2)
+ yrs_now_rng <- range(yrs_now_rng)
+ if (length(yrs_then_rng) > 2)
+ yrs_then_rng <- range(yrs_then_rng)
+
+ r_now <- get_sst_r(dates_now)
+ r_then <- get_sst_r(dates_then)
+
+ lgnd_now <- glue(
+ "Now
+ {format(input$sld_o_md, '%b %d')},
+ {paste(yrs_now_rng, collapse = ' to ')}")
+ lgnd_then <- glue(
+ "Then
+ {format(input$sld_o_md, '%b %d')},
+ {paste(yrs_then_rng, collapse = ' to ')}")
+
+ map_then_now(
+ r_then,
+ r_now,
+ tiles = ifelse(
+ isTRUE(input$dark_mode),
+ providers$CartoDB.DarkMatter,
+ providers$CartoDB.Positron),
+ lgnd_then,
+ lgnd_now,
+ var_lbl)
+
+ })
+
+ # plot_sst ----
+ output$plot_sst <- plotly::renderPlotly({
+
+ bay_segment <- "BCB"
+
+ d <- d_sst_z |>
+ filter(bay_segment == !!bay_segment) |>
+ mutate(
+ year = year(time),
+ date = sprintf(
+ "%d-%02d-%02d",
+ year(today()), month(time), day(time) ) |>
+ as.POSIXct(),
+ color = case_when(
+ year == year(today()) ~ "red",
+ year == year(today()) - 1 ~ "orange",
+ .default = "gray") ) |>
+ select(time, year, date, color, val) |>
+ arrange(year, date, val)
+
+ yrs <- as.character(unique(d$year))
+ colors <- setNames(rep("darkgray", length(yrs)), yrs)
+ colors[as.character(year(today()))] <- "red"
+ colors[as.character(year(today()) - 1)] <- "orange"
+
+ d <- d |>
+ group_by(year) |>
+ mutate(
+ val_sl = slider::slide_mean(
+ val, before = 3L, after = 3L, step = 1L,
+ complete = F, na_rm = T),
+ txt_date = as.Date(time),
+ txt_val = round(val_sl, 2) ) |>
+ select(-time) |>
+ ungroup()
+
+ # TODO: darkly theme w/ bslib
+ g <- ggplot(
+ d,
+ aes(
+ x = date,
+ y = val_sl,
+ group = year,
+ color = factor(year),
+ date = txt_date,
+ value = txt_val)) + # frame = yday
+ geom_line(
+ # aes(text = text),
+ alpha = 0.6) +
+ scale_colour_manual(
+ values = colors) +
+ theme(legend.position = "none") +
+ scale_x_datetime(
+ labels = date_format("%b %d")) +
+ labs(
+ x = "Day of year",
+ y = "Temperature ºC")
+ g
+ # x, y, alpha, color, group, linetype, size
+
+ # add color theming
+ # https://rstudio.github.io/thematic/articles/auto.html
+
+ ggplotly(g, tooltip=c("date","value"))
+ })
+
+}
diff --git a/app/ui.R b/app/ui.R
new file mode 100644
index 0000000..172fa84
--- /dev/null
+++ b/app/ui.R
@@ -0,0 +1,197 @@
+page_navbar(
+ title = "Tampa Bay Climate Change Indicators",
+ theme = dark,
+
+ nav_panel(
+ title = tagList(
+ "Air Temperature", bs_icon("thermometer-half")),
+
+ card(
+ card_header(
+ "Air Temperature - Map of Then vs Now",
+ popover(
+ title = "Settings",
+ bs_icon("gear", class = "ms-auto"),
+ selectInput(
+ "sel_t_var",
+ "Variable",
+ c("Max" = "tmax",
+ "Average" = "tmean",
+ "Min" = "tmin"),
+ selected = "tmax"),
+ sliderInput(
+ "sld_t_md",
+ "Month and day of year",
+ min = as.Date(glue("{year(now_prism)}-01-01")),
+ value = now_prism,
+ max = as.Date(glue("{year(now_prism)}-12-31")),
+ timeFormat = "%b %d",
+ animate = T),
+ sliderInput(
+ "sld_t_yrs_then",
+ "Then year(s)",
+ min = yrs_prism[1],
+ value = c(yrs_prism[1], yrs_prism[1]+20),
+ max = yrs_prism[2] - 1,
+ step = 1,
+ animate = T,
+ sep = ""),
+ sliderInput(
+ "sld_t_yrs_now",
+ "Now year(s)",
+ min = yrs_prism[1] + 1,
+ value = c(yrs_prism[2], yrs_prism[2]),
+ max = yrs_prism[2],
+ step = 1,
+ animate = T,
+ sep = "")) ),
+
+ leafletOutput("map_prism_temp") ) ),
+
+ nav_panel(
+ title = tagList(
+ "Rain", bs_icon("cloud-rain-fill")),
+
+ card(
+ card_header(
+ "Precipitation - Map of Then vs Now",
+ popover(
+ title = "Settings",
+ bs_icon("gear", class = "ms-auto"),
+ sliderInput(
+ "sld_p_md",
+ "Month and day of year",
+ min = as.Date(glue("{year(now_prism)}-01-01")),
+ value = now_prism,
+ max = as.Date(glue("{year(now_prism)}-12-31")),
+ timeFormat = "%b %d",
+ animate = T),
+ sliderInput(
+ "sld_p_yrs_then",
+ "Then year(s)",
+ min = yrs_prism[1],
+ value = c(yrs_prism[1], yrs_prism[1]+20),
+ max = yrs_prism[2] - 1,
+ step = 1,
+ animate = T,
+ sep = ""),
+ sliderInput(
+ "sld_p_yrs_now",
+ "Now year(s)",
+ min = yrs_prism[1] + 1,
+ value = c(yrs_prism[2], yrs_prism[2]),
+ max = yrs_prism[2],
+ step = 1,
+ animate = T,
+ sep = "")) ),
+
+ leafletOutput("map_prism_ppt") ) ),
+
+ nav_panel(
+ title = tagList(
+ "Sea Level", bs_icon("water")) ),
+
+ nav_panel(
+ title = tagList(
+ "Ocean Temperature", bs_icon("thermometer-low")),
+
+ navset_card_underline(
+ nav_panel(
+ "Map",
+
+ card(
+ card_header(
+ "Sea Surface Temperate - Map of Then vs Now",
+ popover(
+ title = "Settings",
+ bs_icon("gear", class = "ms-auto"),
+ sliderInput(
+ "sld_o_md",
+ "Month and day of year",
+ min = as.Date(glue("{year(now_prism)}-01-01")),
+ value = now_prism,
+ max = as.Date(glue("{year(now_prism)}-12-31")),
+ timeFormat = "%b %d",
+ animate = T),
+ sliderInput(
+ "sld_o_yrs_then",
+ "Then year(s)",
+ min = yrs_prism[1],
+ value = c(yrs_prism[1], yrs_prism[1]+20),
+ max = yrs_prism[2] - 1,
+ step = 1,
+ animate = T,
+ sep = ""),
+ sliderInput(
+ "sld_o_yrs_now",
+ "Now year(s)",
+ min = yrs_prism[1] + 1,
+ value = c(yrs_prism[2], yrs_prism[2]),
+ max = yrs_prism[2],
+ step = 1,
+ animate = T,
+ sep = "")) ),
+
+ leafletOutput("map_sst") ) ),
+
+ nav_panel(
+ "Plot",
+
+ card(
+ card_header(
+ "Sea Surface Temperate - Time Series",
+ # popover(
+ # title = "Settings",
+ # bs_icon("gear", class = "ms-auto"),
+ # sliderInput(
+ # "sld_o_md",
+ # "Month and day of year",
+ # min = as.Date(glue("{year(now_prism)}-01-01")),
+ # value = now_prism,
+ # max = as.Date(glue("{year(now_prism)}-12-31")),
+ # timeFormat = "%b %d",
+ # animate = T),
+ # sliderInput(
+ # "sld_o_yrs_then",
+ # "Then year(s)",
+ # min = yrs_prism[1],
+ # value = c(yrs_prism[1], yrs_prism[1]+20),
+ # max = yrs_prism[2] - 1,
+ # step = 1,
+ # animate = T,
+ # sep = ""),
+ # sliderInput(
+ # "sld_o_yrs_now",
+ # "Now year(s)",
+ # min = yrs_prism[1] + 1,
+ # value = c(yrs_prism[2], yrs_prism[2]),
+ # max = yrs_prism[2],
+ # step = 1,
+ # animate = T,
+ # sep = ""))
+ ),
+
+ plotlyOutput("plot_sst") ) )
+ )
+ ),
+
+ nav_spacer(),
+ nav_menu(
+ title = "More",
+ align = "right",
+ nav_item(
+ tags$a(
+ shiny::icon("info-circle"),
+ "About",
+ href = "https://tbep-tech.github.io/climate-change-indicators/",
+ target = "_blank") ),
+ nav_item(
+ tags$a(
+ shiny::icon("envelope", class = "fa-solid"),
+ "Contact",
+ href = "https://tbep-tech.github.io/climate-change-indicators/",
+ target = "_blank") ) ),
+
+ nav_item(
+ input_switch("dark_mode", bs_icon("moon-stars-fill"), TRUE) )
+)
diff --git a/data.qmd b/data.qmd
index 88674d6..ef7c43a 100644
--- a/data.qmd
+++ b/data.qmd
@@ -7,7 +7,7 @@ editor_options:
# Data
Themes for indicators seeking relevance to:
-- Land & Sea
+- Air & Sea
- People & Ecosystems
Icons
@@ -68,7 +68,9 @@ library(webshot2)
options(readr.show_col_types = F)
```
-## Temperature
+## Air
+
+noaacrwsstDaily
### Observed
@@ -387,7 +389,7 @@ p
```
-### prism
+## PRISM
- [HMS: Hydrologic Micro Services \| United States Environmental Protection Agency \| US EPA](https://qed.epa.gov/hms/meteorology/humidity/algorithms/)
@@ -434,37 +436,94 @@ Temporal data availability:
> Daily data are considered "early" for the current month. The previous six months are provisional data. After six months data are considered stable. Thus early data only exist for daily data, while there can be monthly (and presumably yearly) provisional data. -- `?prism::prism_archive_clean`
+### Example `tbeptools::read_importprism()`
+
+- `@example` from [`tbeptools::read_importprism()`](https://github.com/tbep-tech/tbeptools/blob/climate-change-indicators/R/read_importprism.R)
+
```{r}
-#| label: get_prism
+#| label: example_read_importprism
#| messages: false
#| warning: false
-# librarian::shelf(glue, here, leaflet, lubridate,
-# # mapview, stars,
-# prism, sf, stringr, tbep-tech/tbeptools, terra)
+# devtools::install_local("../tbeptools", force = T)
+# devtools::load_all("../tbeptools")
+
+# setup output directory and table
+dir_tif <- here::here("../tbeptools/inst/prism")
+zonal_csv <- file.path(dir_tif, "_zones.csv")
+
+# run function for Tampa Bay watersheds for first 3 days and 4 variables
+d <- read_importprism(
+ vars = c("tmin", "tmax", "tdmean", "ppt"),
+ date_beg = as.Date("1981-01-01"),
+ date_end = as.Date("1981-01-03"),
+ dir_tif = dir_tif,
+ sf_zones = tbsegshed,
+ fld_zones = "bay_segment",
+ zonal_csv = zonal_csv)
+
+# plot first of output rasters
+tifs <- list.files(dir_tif, pattern = ".tif$", full.names = T)
+basename(tifs)
+r <- terra::rast(tifs[1])
+names(r) # {data observed}_{variable}_v{version}_{date updated}
+terra::plet(
+ r[[3]],
+ main = names(r)[3],
+ col = "Spectral",
+ tiles = "CartoDB.DarkMatter")
+
+# show summary by zone
+d
+```
-# renv register libraries
-library(glue)
-library(here)
-library(leaflet)
-library(lubridate)
-# library(mapview)
-library(prism)
-# library(stars)
-library(sf)
-library(stringr)
-library(tbeptools)
-library(terra)
+### Operationalize
-dir_prism <- here("tmp/prism")
-dir.create(dir_prism, showWarnings = F)
-prism_set_dl_dir(dir_prism)
+```{r}
+#| label: operationalize_read_importprism
+#| eval: false
+
+librarian::shelf(
+ dplyr, here, sf, tbep-tech/tbeptools)
+
+# outputs
+dir_prism <- here("data/prism")
+prism_csv <- here("data/prism.csv")
+
+# Tampa Bay watershed zones, including whole bay ("TB")
+tb_zones <- tbsegshed |>
+ bind_rows(
+ tbshed |>
+ mutate(bay_segment = "TB") |>
+ select(bay_segment, geometry))
+
+# bounding box with 0.2º margin around watershed, rounded to .1
+bb <- tbshed |>
+ st_buffer(0.2) |>
+ st_bbox() |>
+ round(1)
+# c(xmin = -82.9, ymin = 27.4, xmax = -81.9, ymax = 28.4)
+# after running terra::trim() on raster:
+bb <- c(xmin = -82.9, ymin = 27.2, xmax = -81.7, ymax = 28.6)
+
+d <- read_importprism(
+ vars = c("tmin", "tmax", "tdmean", "ppt"),
+ date_beg = as.Date("1981-01-01"),
+ date_end = Sys.Date(),
+ bbox = bb,
+ dir_tif = dir_prism,
+ sf_zones = tb_zones,
+ fld_zones = "bay_segment",
+ zonal_csv = prism_csv,
+ verbose = T)
+
+# show summary by zone
+d
```
-- https://prism.nacse.org/downloads/
- - https://prism.nacse.org/documents/PRISM_downloads_web_service.pdf
- - http://services.nacse.org/prism/data/public/4km/tmin/20090405
+TODO: Github Action on daily cron to run `read_importprism()` for all days and variables
+```{r}
# Get monthly (every month) and annual 30-year normals for precipitation
# get_prism_normals(
@@ -479,7 +538,6 @@ prism_set_dl_dir(dir_prism)
# keepZip = FALSE)
-```{r}
var <- "tmax"
period <- "daily"
date <- today() - days(1) # yesterday
@@ -1296,7 +1354,9 @@ Questions:
- news summary: [A new way to quantify climate change impacts: “Outdoor days” \| MIT News \| Massachusetts Institute of Technology](https://news.mit.edu/2024/new-way-quantify-climate-change-impacts-outdoor-days-0322)
- Shiny app: [California Outdoor Days \| Eltahir Research Group](https://eltahir.mit.edu/california-outdoor-days/)
-## Sea Level Rise
+## Sea
+
+### Sea Level
Sea level rise occurs from principally two sources: 1) thermal expansion; and 2) freshwater inputs from glacial melting. Data for these trends can be obtained from NOAA's [Sea Level Trends](https://tidesandcurrents.noaa.gov/sltrends/sltrends.html) (@fig-slr-noaa)
@@ -1306,7 +1366,7 @@ Types of data:
- Level-3 products distributed through NOAA CoastWatch (Sea Level Anomaly and along-track altimetry)
2. Projected (future). modeled
-### Gauges
+#### Gauges
 zoomed into the Tampa Bay.](figures/slr_tidesandcurrents.noaa.gov-sltrends.png){#fig-slr-noaa}
@@ -1329,128 +1389,12 @@ https://tidesandcurrents.noaa.gov/sltrends/sltrends_station.shtml?id=8726724 htt
- [FAQ for Sea Level Trends | NOAA Tides & Currents](https://tidesandcurrents.noaa.gov/sltrends/faq.html#:~:text=Relative%20Sea%20Level%20Trends%20reflect%20changes,stations%2C%20called%20monthly%20mean%20sea%20level.)
> **Relative Sea Level Trends** reflect changes in local sea level over time and are typically the most critical sea level trend for many coastal applications, including coastal mapping, marine boundary delineation, coastal zone management, coastal engineering, sustainable habitat restoration design, and the general public enjoying their favorite beach. This website focuses on relative sea level trends, computed from monthly averages of hourly water levels observed at specific tide stations, called **monthly mean sea level**.
-- [Tidal Datums | NOAA Tides & Currents](https://tidesandcurrents.noaa.gov/datum_options.html)
-- Water Level Monthly Mean Data
-- Date Time Date and time of the data
-- MHHW Mean Higher-High Water
-- MHW Mean High Water
-- DTL Mean Diurnal Tide Level
-- MTL Mean Tide Level
-- MSL Mean Sea Level
-- MLW Mean Low Water
-- MLLW Mean Lower-Low Water
-- GT Great Diurnal Range
-- MN Mean Range of Tide
-- DHQ Mean Diurnal High Water Inequality
-- DLQ Mean Diurnal Low Water Inequality
-- HWI Greenwich High Water Interval (in Hours)
-- LWI Greenwich Low Water Interval (in Hours)
-- highest Highest Tide
-- lowest Lowest Tide
-- inferred A flag that when set to 1 indicates that the water level value has been inferred
-
+See `tbeptools` additions:
-```{r}
-
-
-library(glue)
-
-station_id
-slr_csv <- "https://tidesandcurrents.noaa.gov/sltrends/data/{station_id}_meantrend.csv"
-slr_csv <- "https://tidesandcurrents.noaa.gov/sltrends/data/8726724_meantrend.csv"
-
-rm(d)
-d <- suppressWarnings(
- request("https://tidesandcurrents.noaa.gov") |>
- req_url_path_append(
- "sltrends", "data",
- glue("{station_id}_meantrend.csv")) |>
- req_perform() |>
- resp_body_string() |>
- read_csv(
- skip = 5,
- name_repair = "unique_quiet",
- show_col_types = F) )
-# names(d) |> paste(collapse = ", ")
-# Year, Month, Monthly_MSL, Linear_Trend, High_Conf., Low_Conf.
-
-problems(d)
-
-```
-
-https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/stations/id.extension
-https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/stations/9414290.json https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/stations/id/resource.extension https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/stations/8726724/datums.json "units": "feet", "datums": "name": "MSL", "description": "Mean Sea Level", "value": 3.2 "min": -0.84, "mindate": "19770119", "mintime": "08:06", "max": 8.491, "maxdate": "19930313", "maxtime": "04:48",
-
-https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/ports.extension https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/ports.json
-
-```{r}
-librarian::shelf(
- dplyr, httr2, purrr, readr, tibble)
-# # register with renv
-# library(dplyr)
-# library(httr2)
-# library(readr)
-# # library(rnoaa)
-# library(tibble)
-
-d_csv <- here::here("data/slr/monthly_sea_levels.csv")
-
-d_stations <- tribble(
- ~station_id, ~station_name,
- 8726724 , "Clearwater Beach",
- 8726674 , "East Bay",
- 8726384 , "Port Manatee",
- 8726520 , "St. Petersburg")
-
-get_station_data <- function(station_id){
- httr2::request("https://api.tidesandcurrents.noaa.gov/api/prod/datagetter") |>
- httr2::req_url_query(
- station = station_id,
- begin_date = 19010101,
- end_date = lubridate::today() |> format("%Y%m%d") |> as.integer(),
- product = "monthly_mean", # maxdur = 3653 days (10 yrs)
- datum = "stnd",
- time_zone = "lst", # local time
- units = "metric", # meters
- format = "csv") |>
- httr2::req_perform() |>
- httr2::resp_body_string() |>
- readr::read_csv(
- name_repair = "unique_quiet",
- show_col_types = F)
-}
-
-d <- d_stations |>
- dplyr::mutate(
- data = purrr::map(station_id, get_station_data)) |>
- tidyr::unnest(data)
-
-write_csv(d, d_csv)
-
-# d
-# # A tibble: 1,868 × 20
-# station_id station_name Year Month Highest MHHW MHW MSL MTL MLW
-#
-# 1 8726724 Clearwater Beach 1973 5 1.55 1.38 1.26 0.96 0.96 0.658
-# 2 8726724 Clearwater Beach 1973 6 1.61 1.38 1.25 0.951 0.954 0.655
-# 3 8726724 Clearwater Beach 1973 7 1.58 1.38 1.25 0.951 0.942 0.631
-# 4 8726724 Clearwater Beach 1973 8 1.57 1.40 1.28 0.981 0.972 0.668
-# 5 8726724 Clearwater Beach 1973 9 1.54 1.43 1.33 1.04 1.03 0.732
-# 6 8726724 Clearwater Beach 1973 10 1.66 1.38 1.29 1.02 1.02 0.747
-# 7 8726724 Clearwater Beach 1973 11 1.64 1.40 1.29 1.01 1.02 0.753
-# 8 8726724 Clearwater Beach 1973 12 1.53 1.28 1.16 0.896 0.893 0.622
-# 9 8726724 Clearwater Beach 1974 1 1.53 1.30 1.16 0.872 0.866 0.573
-# 10 8726724 Clearwater Beach 1974 2 1.61 1.27 1.15 0.869 0.863 0.579
-# # ℹ 1,858 more rows
-# # ℹ 10 more variables: MLLW , DTL , GT , MN , DHQ ,
-# # DLQ , HWI , LWI , Lowest , Inferred
-# # ℹ Use `print(n = ...)` to see more rows
-# names(d) |> paste(collapse=", ")
-# station_id, station_name, Year, Month, Highest, MHHW, MHW, MSL, MTL, MLW, MLLW, DTL, GT, MN, DHQ, DLQ, HWI, LWI, Lowest, Inferred
-#
-# https://tidesandcurrents.noaa.gov/axis/webservices/waterlevelverifiedmonthly/response.jsp?stationId=8726724&beginDate=20220201&endDate=20240514&datum=MSL&unit=0&timeZone=0&format=html&Submit=Submit
+- [sealevelstations.R](https://github.com/tbep-tech/tbeptools/blob/climate-change-indicators/R/sealevelstations.R)
+- [read_importsealevels.R](https://github.com/tbep-tech/tbeptools/blob/climate-change-indicators/R/read_importsealevels.R)
+- [sealevels.Rmd](https://github.com/tbep-tech/tbeptools/blob/climate-change-indicators/vignettes/sealevels.Rmd)
-```
- [Sea Level Trends - NOAA Tides & Currents](https://tidesandcurrents.noaa.gov/sltrends/sltrends.html)
@@ -1476,18 +1420,7 @@ write_csv(d, d_csv)
- [Inundation Dashboard - NOAA Tides & Currents](https://tidesandcurrents.noaa.gov/inundationdb/)
-### Hurricanes
-
-- property damage
-
-```{r}
-
-rnoaa::coops_search()
-# swdi - Severe Weather Data Inventory (SWDI) vignette
-
-```
-
-### Sea Surface Temperature
+### Surface Temperature
- coral reefs: ecosystem health, ecotourism
@@ -1498,14 +1431,14 @@ rnoaa::coops_search()
- background: [NOAA Coral Reef Watch Tutorial](https://coralreefwatch.noaa.gov/product/5km/tutorial/crw05a_sst_product.php)
- [ERDDAP - Sea Surface Temperature, NOAA Coral Reef Watch Daily Global 5km Satellite SST (CoralTemp), 1985-present, Daily - Data Access Form](https://coastwatch.noaa.gov/erddap/griddap/noaacrwsstDaily.html)
+Output is stored in a single raster file containing about a little over 13K layers, which seems to read and render reasonably fast (whereas for the more numerous PRISM layers, it was way too slow).
```{r}
-
librarian::shelf(
sf, tbeptools, mapview)
tbshed_buf <- tbshed |>
- st_buffer(10000) |>
+ st_buffer(10000) |> # buffer by 10 km
st_bbox() |>
st_as_sfc() |>
st_as_sf()
@@ -1513,12 +1446,283 @@ tbshed_buf <- tbshed |>
# st_bbox(tbshed_buf) |> round(1)
# xmin ymin xmax ymax
# -83.0 27.3 -81.8 28.5
+# after clipping to non-NA raster values
+bbox = c(xmin = -83.0, ymin = 27.2, xmax = -82.3, ymax= 28.5)
+sst_bbox <- bbox |>
+ sf::st_bbox() |>
+ sf::st_as_sfc() |>
+ sf::st_as_sf(crs = "+proj=longlat +datum=WGS84 +no_defs")
+
+mapView(tbsegshed, col.regions="gray") +
+ mapView(tbshed_buf, col.regions = "blue") +
+ mapView(sst_bbox, col.regions = "red")
+```
-mapView(tbshed, col.regions="gray") +
- mapView(tbshed_buf)
+```{r}
+#| label: get-sst
+#| eval: false
+librarian::shelf(
+ marinebon/extractr,
+ here)
+
+ed_extract(
+ ed = ed_info("https://coastwatch.noaa.gov/erddap/griddap/noaacrwsstDaily.html"),
+ var = "analysed_sst",
+ bbox = c(xmin = -83.0, ymin = 27.2, xmax = -82.3, ymax= 28.5),
+ aoi = tbeptools::tbsegshed,
+ zonal_csv = here("data/sst/tbep_sst.csv"),
+ rast_tif = here("data/sst/tbep_sst.tif"),
+ mask_tif = F)
```
+
+```{r}
+librarian::shelf(
+ ggplot2, here, highcharter, lubridate, plotly, readr)
+
+# data("mpg", "diamonds", "economics_long", package = "ggplot2")
+# economics_long2 <- dplyr::filter(economics_long, variable %in% c("pop", "uempmed", "unemploy"))
+# hchart(economics_long2, "line", hcaes(x = date, y = value01, group = variable))
+
+sst_csv <- here("data/sst/tb_sst.csv")
+d_sst <- read_csv(sst_csv, show_col_types = F)
+table(d_sst$bay_segment)
+# BCB HB LTB MR MTB OTB TCB
+# 13169 13169 13169 13169 13169 13169 13169
+
+bay_segment <- "BCB"
+
+d <- d_sst |>
+ filter(bay_segment == !!bay_segment) |>
+ mutate(
+ year = year(time),
+ # yday = yday(time),
+ date = sprintf(
+ "%d-%02d-%02d",
+ year(today()), month(time), day(time) ) |>
+ as.POSIXct(),
+ color = case_when(
+ year == year(today()) ~ "red",
+ year == year(today()) - 1 ~ "orange",
+ .default = "gray") ) |>
+ # select(year, yday, date, color, val) |>
+ # arrange(year, yday, date, val)
+ select(time, year, date, color, val) |>
+ arrange(year, date, val)
+# table(d$yday) |> range()
+
+yrs <- as.character(unique(d$year))
+colors <- setNames(rep("darkgray", length(yrs)), yrs)
+colors[as.character(year(today()))] <- "red"
+colors[as.character(year(today()) - 1)] <- "orange"
+
+librarian::shelf(
+ slider, scales)
+
+d <- d |>
+ group_by(year) |>
+ mutate(
+ val_sl = slider::slide_mean(
+ val, before = 3L, after = 3L, step = 1L,
+ complete = F, na_rm = T),
+ txt_date = as.Date(time),
+ txt_val = round(val_sl, 2) ) |>
+ select(-time) |>
+ ungroup()
+
+# TODO: darkly theme w/ bslib
+g <- ggplot(
+ d,
+ aes(
+ x = date,
+ y = val_sl,
+ group = year,
+ color = factor(year),
+ date = txt_date,
+ value = txt_val)) + # frame = yday
+ geom_line(
+ # aes(text = text),
+ alpha = 0.6) +
+ scale_colour_manual(
+ values = colors) +
+ theme(legend.position = "none") +
+ scale_x_datetime(
+ labels = date_format("%b %d")) +
+ labs(
+ x = "Day of year",
+ y = "Temperature ºC")
+g
+# x, y, alpha, color, group, linetype, size
+
+# add color theming
+
+# https://rstudio.github.io/thematic/articles/auto.html
+
+ggplotly(g, tooltip=c("date","value"))
+
+# https://plotly-r.com/scatter-traces 3
+# https://plotly-r.com/client-side-linking 16
+
+
+# https://stackoverflow.com/questions/76435688/how-to-autoplay-a-plotly-chart-in-shiny
+d |>
+ plot_ly(
+ x = ~date,
+ y = ~val,
+ frame = ~yday,
+ # group = ~year,
+ # color = ~color,
+ type = 'scatter',
+ mode = 'lines',
+ # marker = list(size = 20),
+ showlegend = F,
+ transforms = list(
+ list(
+ type = 'groupby',
+ groups = d$year,
+ styles = list(
+ list(target = 2024, value = list(line =list(color = 'red'))),
+ list(target = 2023, value = list(line =list(color = 'orange'))),
+ list(target = 2022, value = list(line =list(color = 'darkgray'))) ) ) ) )
+
+
+ animation_button(visible = T) |>
+ onRender("
+ function(el,x) {
+ Plotly.animate(el);
+ }")
+
+library(plotly)
+library(htmlwidgets)
+
+df <- data.frame(
+ x = c(1,2,1),
+ y = c(1,2,1),
+ f = c(1,2,3)
+)
+df %>%
+ plot_ly(
+ x = ~x,
+ y = ~y,
+ frame = ~f,
+ type = 'scatter',
+ mode = 'markers',
+ marker = list(size = 20),
+ showlegend = FALSE
+ ) %>%
+ animation_button(visible = TRUE) %>%
+ onRender("
+ function(el,x) {
+ Plotly.animate(el);
+ }")
+
+
+# https://stackoverflow.com/questions/61152879/change-the-frame-label-in-plotly-animation
+DF <- data.frame(
+ year = rep(seq(1980L, 2020L), each = 12),
+ month = rep(1:12, 41),
+ month_char = rep(factor(month.abb), 41),
+ avg_depth = runif(492) )
+# with(DF, paste0(sprintf("%02d", month), " - ", month_char) )
+fig <- DF |>
+ plot_ly(
+ x = ~year,
+ y = ~avg_depth,
+ frame = ~paste0(sprintf("%02d", month), " - ", month_char),
+ type = 'bar') |>
+ animation_slider(
+ currentvalue = list(prefix = "Month: ") )
+
+fig
+
+# https://stackoverflow.com/questions/50843134/r-plotly-animated-chart-only-showing-groups-with-data-in-initial-frame
+
+dates <- 2000:2010
+countries <- c("US", "GB", "JP")
+df <- merge(dates, countries, all=TRUE)
+names(df) <- c("Date", "Country")
+
+
+df$x <- rnorm(nrow(df))
+df$y <- rnorm(nrow(df))
+
+df[1:3, c("x", "y")] <- NA
+
+
+p <- ggplot(df, aes(x, y, color = Country)) +
+ geom_point(aes(frame = Date)) + theme_bw()
+
+ggplotly(p)
+
+# other
+
+g <- ggplot(
+ d,
+ aes(
+ x = date,
+ y = val,
+ color = factor(year))) +
+ geom_line(alpha = 0.6) +
+ scale_colour_manual(
+ values = colors) +
+ theme(legend.position = "none")
+g
+
+
+p <- ggplotly(g) |>
+ onRender("
+ function(el,x) {
+ Plotly.animate(el);
+ }")
+p
+
+d |>
+ hchart(
+ "line", # "line"
+ hcaes(
+ x = date,
+ y = val,
+ group = year,
+ segmentColor = clr),
+ showInLegend = F)
+
+mpgman2 <