library(extractr) # devtools::install_local(force = T) # devtools::load_all()
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#>      (status 2 uses the sf package in place of rgdal)

Get regions of interests

if (!require("librarian"))
  install.packages("librarian")
#> Loading required package: librarian

librarian::shelf(
  dplyr, here, mapview, sf)

fk_roi_geo <- "https://github.com/marinebon/extract-app/raw/main/data/fk_roi.geojson"
fk_sta_geo <- "https://github.com/marinebon/extract-app/raw/main/data/fk_sta.geojson"

fk_roi_sf <- read_sf(fk_roi_geo)
fk_sta_sf <- read_sf(fk_sta_geo)
fk_sf <- bind_rows(fk_roi_sf, fk_sta_sf)

mapview(fk_sf)

Get ERDDAP dataset information

ed_info <- get_ed_info("https://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURSST41mday.html")
ed_info
#> <ERDDAP info> jplMURSST41mday 
#>  Base URL: https://coastwatch.pfeg.noaa.gov/erddap 
#>  Dataset Type: griddap 
#>  Dimensions (range):  
#>      time: (2002-06-16T00:00:00Z, 2023-08-16T00:00:00Z) 
#>      latitude: (-89.99, 89.99) 
#>      longitude: (-179.99, 180.0) 
#>  Variables:  
#>      mask: 
#>      nobs: 
#>      sst: 
#>          Units: degree_C

Extract

dir_tif <- here::here("data_tmp")
grds <- get_ed_grds(
  ed_info, ed_var="sst", ply=fk_sf, dir_tif=dir_tif, date_beg = "2021-10-01")
#> Found 22 dates between 2021-10-01 and 2023-08-16.
#> Found 22 matching tifs of 22 dates.
#> Reading existing grids ([dir_tif]/grd_sst_[date].tif) vs fetching fresh data via ERDDAP.
# Found 17 dates between 2021-10-01 and 2023-03-16

Summarize

ts_csv <- here::here("data_tmp/ts.csv")
ts <- grds_to_ts(
  grds, fxns = c("mean", "sd"), ts_csv, verbose = T)
#> Reading ts_csv
#> All grds_dates found in ts_csv$date, so returning vs recalculating.
head(ts)
#> # A tibble: 6 × 4
#>   lyr            date        mean    sd
#>   <chr>          <date>     <dbl> <dbl>
#> 1 sst_2021.10.16 2021-10-16  29.0 0.116
#> 2 sst_2021.11.16 2021-11-16  27.2 0.945
#> 3 sst_2021.12.16 2021-12-16  26.6 1.12 
#> 4 sst_2022.01.16 2022-01-16  25.7 1.25 
#> 5 sst_2022.02.16 2022-02-16  25.0 1.64 
#> 6 sst_2022.03.16 2022-03-16  26.1 0.948

Plot

plot_ts(ts_csv, "mean", "sd")