suppressPackageStartupMessages({
  # from CRAN
  library(dplyr)
  library(sf)
  library(stars)
  library(readr)
  
  library(rerddap)
  library(ncdf4)
  library(leaflet)
  library(ggplot2)
  
  # from Github
  library(xyzt)
  library(ghrsst)
})

1) ERDDAP#

rerddap

Get a listing of servers that rerddap knows about#

s <- servers()

s

s |> 
  filter(short_name == "OOI")

What services are available in the server?#

ds <- ed_datasets(url="http://erddap.dataexplorer.oceanobservatories.org/erddap/")

paste("We have", sum(ds$tabledap!=""), "tabledap", sum(ds$griddap!=""), "griddap, and", sum(ds$wms!=""), "wms.", sep=" ")

Let’s query all the datasets that have the standard_name of sea_water_practical_salinity.#

out <- ed_search(query = "sea_water_practical_salinity", 
                 which = "tabledap",
                 url="http://erddap.dataexplorer.oceanobservatories.org/erddap/")

out

paste("Found", length(out$alldata), "datasets", sep=" ")

Let us narrow our search to deployments that are within a lon/lat/time extent.#

leaflet

lon <- c(-72, -69)
lat <- c(38, 41)

leaflet() |>
  addTiles() |>
  setView(lng=sum(lon)/2, 
          lat=sum(lat)/2, 
          zoom=6) |>
  addRectangles(lng1 = lon[1],
                lng2 = lon[2],
                lat1 = lat[1],
                lat2 = lat[2])

out_bb <- ed_search_adv(standard_name= "sea_water_practical_salinity",
                        which = "tabledap",
                        url="http://erddap.dataexplorer.oceanobservatories.org/erddap/",
                        minLat=38,
                        maxLat=41,
                        minLon=-72,
                        maxLon=-69,
                        minTime = "2016-07-10T00:00:00Z",
                        maxTime = "2017-02-10T00:00:00Z")

dataset_ids <- out_bb$info$dataset_id

out_bb

paste("Found", length(out_bb$alldata), "Datasets:", sep=" ")
dataset_ids

dataset_ids[1]

info(dataset_ids[1], 
     url="http://erddap.dataexplorer.oceanobservatories.org/erddap/")

zz <- tabledap(x = dataset_ids[1],
               fields = c("z", 
                          "latitude", 
                          "longitude", 
                          "sea_water_practical_salinity_profiler_depth_enabled", 
                          "sea_water_temperature_profiler_depth_enabled", 
                          "time"),
               url = "http://erddap.dataexplorer.oceanobservatories.org/erddap")

head(zz)

summary(zz)

zz_sub <- zz |>
  tibble() |>
  mutate(time = as.POSIXct(time),
         z = as.numeric(z),
         temperature = as.numeric(sea_water_temperature_profiler_depth_enabled),
         salinity = as.numeric(sea_water_practical_salinity_profiler_depth_enabled)) |>
  select(-sea_water_temperature_profiler_depth_enabled, -sea_water_practical_salinity_profiler_depth_enabled) |>
  filter(between(time, as.POSIXct("2016-01-01 00:00:00"), as.POSIXct("2016-12-31 00:00:00")),
         z == min(z))

summary(zz_sub)

plot(x=zz_sub$time, y=zz_sub$salinity)

plot(x=zz_sub$salinity, 
     y=zz_sub$temperature)

ggplot(data=zz_sub, aes(x=salinity, y=temperature)) +
  geom_point()

2) OPeNDAP#

NetCDF4 Cheat Sheet

Bigelow Homegrown R Packages#

Intro Video

xyzt ghrsst

url <- "https://opendap.jpl.nasa.gov/opendap/OceanTemperature/ghrsst/data/GDS2/L4/GLOB/JPL/MUR/v4.1/2022/221/20220809090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc"

X <- nc_open(url)

Retrieve the variables contained in the netcdf object#

mur_vars(X)

Retrieve the spatial resolution stated in the global metadata#

mur_res(X)

Read in CalCOFI Station Data#

x <- xyzt::read_calcofi()

x

Change x to be an sf object containing points#

points <- x |> 
  as_POINT()

covars <- ghrsst::extract(points, 
                          X, 
                          varname = "analysed_sst")

(y <- dplyr::bind_cols(x, covars))

y

summary(y)

Change x to be an sf object containing a bounding box#

bbox <- x |>
  as_BBOX()

Extract the bounding box#

layer <- ghrsst::extract(bbox, X, 
                          varname = "analysed_sst")

layer

Plot the layer we extracted and the points on top#

par(mfrow = c(1,2))
plot(layer, attr = 'analysed_sst', col = sf.colors(n=6), axes = TRUE, reset = FALSE)
plot(sf::st_geometry(points), add = TRUE, col = "black", pch = 19, cex = 1)

3) WMS#

Learning Objectives:#

  • Add a WMS layer to an interactive map. (“GEBCO - General Bathymetric Chart of the Oceans”)

leaflet() |>
  addTiles() |>
  setView(lng = -88,
          lat = 30,
          zoom = 5) |>
  addWMSTiles(
    "https://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv?",
    layers = "GEBCO_latest",
    options = WMSTileOptions(format = "image/png", transparent = TRUE),
    attribution="Weather data © 2012 IEM Nexrad")