extractr

Load packages and set up the environment.

library(extractr)

suppressPackageStartupMessages({
  library(dplyr)
  library(here)
  library(mapview)
  library(readr)
  library(sf)
  library(terra)
})
options(readr.show_col_types = F)

Get area of interest (AoI)

Load area of interest (AoI) from this sanctuaries.geojson for the Florida Keys National Marine Sanctuary (FKNMS).

aoi_geo <- "https://raw.githubusercontent.com/noaa-onms/onmsR/master/data-raw/sanctuaries.geojson"

aoi <- read_sf(aoi_geo) |> 
  filter(nms == "FKNMS")
(bb <- st_bbox(aoi))
#>      xmin      ymin      xmax      ymax 
#> -83.14989  24.30041 -80.06647  25.65046

mapView(aoi)

Dataset: sea surface temperature

Get info

ERDDAP dataset:

• ERDDAP - Sea Surface Temperature, NOAA Coral Reef Watch Daily Global 5km Satellite SST (CoralTemp), 1985-present, Daily - Data Access Form

ed_url  <- "https://coastwatch.noaa.gov/erddap/griddap/noaacrwsstDaily.html"
(ed <- ed_info(ed_url))
#> <ERDDAP info> noaacrwsstDaily 
#>  Base URL: https://coastwatch.noaa.gov/erddap 
#>  Dataset Type: griddap 
#>  Dimensions (range):  
#>      time: (1985-01-01T12:00:00Z, 2025-05-20T12:00:00Z) 
#>      latitude: (-89.975, 89.975) 
#>      longitude: (-179.975, 179.975) 
#>  Variables:  
#>      analysed_sst: 
#>          Units: degree_C 
#>      sea_ice_fraction: 
#>          Units: 1

Get dimensions

var <- "analysed_sst"

dims <- ed_dims(ed)
names(dims)
#> [1] "time"      "latitude"  "longitude"

# let's limit to the most recent 10 times
(times <- tail(dims[["time"]], 10))
#>  [1] "2025-05-11 12:00:00 UTC" "2025-05-12 12:00:00 UTC"
#>  [3] "2025-05-13 12:00:00 UTC" "2025-05-14 12:00:00 UTC"
#>  [5] "2025-05-15 12:00:00 UTC" "2025-05-16 12:00:00 UTC"
#>  [7] "2025-05-17 12:00:00 UTC" "2025-05-18 12:00:00 UTC"
#>  [9] "2025-05-19 12:00:00 UTC" "2025-05-20 12:00:00 UTC"

Extract dataset from AoI

d_csv <- here("data_tmp/sst_timeseries.csv")
r_tif <- here("data_tmp/sst_raster.tif")

(d <- ed_extract(
  ed        = ed,
  var       = var,
  sf_zones  = aoi,
  fld_zones = "nms",
  zonal_fun = "mean",
  zonal_csv = d_csv,
  rast_tif  = r_tif,
  time_min  = times[1],
  verbose   = TRUE))
#> Downloading 1 requests, up to 59 time slices each
#> Fetching request 1 of 1 (2025-05-11 to 2025-05-20) ~ 18:39:48 UTC
#> # A tibble: 10 × 4
#>    nms   lyr                               mean time               
#>    <chr> <chr>                            <dbl> <dttm>             
#>  1 FKNMS analysed_sst|2025-05-11 12:00:00  27.5 2025-05-11 12:00:00
#>  2 FKNMS analysed_sst|2025-05-12 12:00:00  27.4 2025-05-12 12:00:00
#>  3 FKNMS analysed_sst|2025-05-13 12:00:00  27.4 2025-05-13 12:00:00
#>  4 FKNMS analysed_sst|2025-05-14 12:00:00  27.6 2025-05-14 12:00:00
#>  5 FKNMS analysed_sst|2025-05-15 12:00:00  27.9 2025-05-15 12:00:00
#>  6 FKNMS analysed_sst|2025-05-16 12:00:00  28.2 2025-05-16 12:00:00
#>  7 FKNMS analysed_sst|2025-05-17 12:00:00  28.3 2025-05-17 12:00:00
#>  8 FKNMS analysed_sst|2025-05-18 12:00:00  28.3 2025-05-18 12:00:00
#>  9 FKNMS analysed_sst|2025-05-19 12:00:00  28.6 2025-05-19 12:00:00
#> 10 FKNMS analysed_sst|2025-05-20 12:00:00  28.7 2025-05-20 12:00:00

Plot time series

d <- read_csv(d_csv)
head(d)
#> # A tibble: 6 × 4
#>   nms   lyr                               mean time               
#>   <chr> <chr>                            <dbl> <dttm>             
#> 1 FKNMS analysed_sst|2025-05-11 12:00:00  27.5 2025-05-11 12:00:00
#> 2 FKNMS analysed_sst|2025-05-12 12:00:00  27.4 2025-05-12 12:00:00
#> 3 FKNMS analysed_sst|2025-05-13 12:00:00  27.4 2025-05-13 12:00:00
#> 4 FKNMS analysed_sst|2025-05-14 12:00:00  27.6 2025-05-14 12:00:00
#> 5 FKNMS analysed_sst|2025-05-15 12:00:00  27.9 2025-05-15 12:00:00
#> 6 FKNMS analysed_sst|2025-05-16 12:00:00  28.2 2025-05-16 12:00:00

plot_ts(d, label = "Surface Temperature (ºC)")

Map raster

r <- rast(r_tif)
names(r)
#>  [1] "analysed_sst|2025-05-11 12:00:00" "analysed_sst|2025-05-12 12:00:00"
#>  [3] "analysed_sst|2025-05-13 12:00:00" "analysed_sst|2025-05-14 12:00:00"
#>  [5] "analysed_sst|2025-05-15 12:00:00" "analysed_sst|2025-05-16 12:00:00"
#>  [7] "analysed_sst|2025-05-17 12:00:00" "analysed_sst|2025-05-18 12:00:00"
#>  [9] "analysed_sst|2025-05-19 12:00:00" "analysed_sst|2025-05-20 12:00:00"

lyr <- names(r)[1]
plet(r[lyr], tiles = "Esri.OceanBasemap")

Dataset: sea surface salinity

Get info

ERDDAP dataset:

• ERDDAP - Sea Surface Salinity, Miras SMOS, Near Real-Time, Global 0.25°, 2010-present, 3 Day Composite - Data Access Form

ed_url  <- "https://coastwatch.noaa.gov/erddap/griddap/noaacwSMOSsss3day.html"
ed <- ed_info(ed_url)
ed
#> <ERDDAP info> noaacwSMOSsss3day 
#>  Base URL: https://coastwatch.noaa.gov/erddap 
#>  Dataset Type: griddap 
#>  Dimensions (range):  
#>      time: (2010-06-03T12:00:00Z, 2025-05-13T12:00:00Z) 
#>      altitude: (0.0, 0.0) 
#>      latitude: (-89.875, 89.875) 
#>      longitude: (-179.875, 179.875) 
#>  Variables:  
#>      sss: 
#>          Units: PSU 
#>      sss_dif: 
#>          Units: PSU

Get dimensions

var <- "sss"

dims <- ed_dims(ed)
names(dims)
#> [1] "time"      "altitude"  "latitude"  "longitude"

# let's limit to the most recent 10 times
(times <- tail(dims[["time"]], 10))
#>  [1] "2025-04-16 12:00:00 UTC" "2025-04-19 12:00:00 UTC"
#>  [3] "2025-04-22 12:00:00 UTC" "2025-04-25 12:00:00 UTC"
#>  [5] "2025-04-28 12:00:00 UTC" "2025-05-01 12:00:00 UTC"
#>  [7] "2025-05-04 12:00:00 UTC" "2025-05-07 12:00:00 UTC"
#>  [9] "2025-05-10 12:00:00 UTC" "2025-05-13 12:00:00 UTC"

Extract dataset from AoI

d_csv <- here("data_tmp/sss_timeseries.csv")
r_tif <- here("data_tmp/sss_raster.tif")

(d <- ed_extract(
  ed        = ed,
  var       = var,
  sf_zones  = aoi,
  fld_zones = "nms",
  zonal_fun = "mean",
  zonal_csv = d_csv,
  rast_tif  = r_tif,
  time_min  = times[1],
  verbose   = TRUE))
#> Downloading 1 requests, up to 1282 time slices each
#> Fetching request 1 of 1 (2025-04-16 to 2025-05-13) ~ 18:39:55 UTC
#> # A tibble: 10 × 4
#>    nms   lyr                      mean time               
#>    <chr> <chr>                   <dbl> <dttm>             
#>  1 FKNMS sss|2025-04-16 12:00:00  36.7 2025-04-16 12:00:00
#>  2 FKNMS sss|2025-04-19 12:00:00  37.0 2025-04-19 12:00:00
#>  3 FKNMS sss|2025-04-22 12:00:00  36.1 2025-04-22 12:00:00
#>  4 FKNMS sss|2025-04-25 12:00:00  36.9 2025-04-25 12:00:00
#>  5 FKNMS sss|2025-04-28 12:00:00  34.8 2025-04-28 12:00:00
#>  6 FKNMS sss|2025-05-01 12:00:00  36.7 2025-05-01 12:00:00
#>  7 FKNMS sss|2025-05-04 12:00:00  38.0 2025-05-04 12:00:00
#>  8 FKNMS sss|2025-05-07 12:00:00  36.3 2025-05-07 12:00:00
#>  9 FKNMS sss|2025-05-10 12:00:00  36.5 2025-05-10 12:00:00
#> 10 FKNMS sss|2025-05-13 12:00:00  35.7 2025-05-13 12:00:00

Plot time series

d <- read_csv(d_csv)
head(d)
#> # A tibble: 6 × 4
#>   nms   lyr                      mean time               
#>   <chr> <chr>                   <dbl> <dttm>             
#> 1 FKNMS sss|2025-04-16 12:00:00  36.7 2025-04-16 12:00:00
#> 2 FKNMS sss|2025-04-19 12:00:00  37.0 2025-04-19 12:00:00
#> 3 FKNMS sss|2025-04-22 12:00:00  36.1 2025-04-22 12:00:00
#> 4 FKNMS sss|2025-04-25 12:00:00  36.9 2025-04-25 12:00:00
#> 5 FKNMS sss|2025-04-28 12:00:00  34.8 2025-04-28 12:00:00
#> 6 FKNMS sss|2025-05-01 12:00:00  36.7 2025-05-01 12:00:00

plot_ts(d, label = "Surface Salinity (PSU)")

Map raster

r <- rast(r_tif)
names(r)
#>  [1] "sss|2025-04-16 12:00:00" "sss|2025-04-19 12:00:00"
#>  [3] "sss|2025-04-22 12:00:00" "sss|2025-04-25 12:00:00"
#>  [5] "sss|2025-04-28 12:00:00" "sss|2025-05-01 12:00:00"
#>  [7] "sss|2025-05-04 12:00:00" "sss|2025-05-07 12:00:00"
#>  [9] "sss|2025-05-10 12:00:00" "sss|2025-05-13 12:00:00"

lyr <- names(r)[1]
plet(r[lyr], tiles = "Esri.OceanBasemap")