Discovering ocean data with AQUAVIEW¶

PAM–Glider Rodeo Hackweek · hands-on tutorial (R)¶

AQUAVIEW is one common catalog over 600,000+ ocean datasets from around a hundred sources — satellite, ship and float profiles, gliders, buoys, models, biodiversity, vessel traffic, and more. Instead of tracking down data across many separate servers, clouds and portals, you discover and reach it all in one place — by website, by code, or by asking an AI.

This notebook is in two parts:

  1. Meet AQUAVIEW — what it is, and the three ways to reach it.
  2. A worked example — find a real Hawaiʻi glider, pull its track, and check the satellite against it.

By the end you will have reproduced a genuine validation result: satellite SST vs glider SST along a two-week track, with bias, RMSE and correlation.

This is the exact mirror of the Python notebook — same catalog, same glider, same numbers.

What you need: rstac, readr, dplyr, ncdf4, ggplot2. On the Glider Rodeo Hub these are already installed.

In [1]:
suppressMessages({
  library(rstac)      # STAC client
  library(readr)      # fast CSV reader
  library(dplyr)      # data wrangling
  library(ncdf4)      # NetCDF / OPeNDAP
  library(ggplot2)    # plots
})
options(repr.plot.width = 10, repr.plot.height = 4)

# --- two small display helpers, used throughout -------------------------------
# Trim a long title at a word boundary rather than mid-word.
short <- function(x, n = 44) {
  x <- ifelse(is.na(x), "", as.character(x))
  vapply(x, function(s) {
    if (nchar(s) <= n) return(s)
    cut <- substr(s, 1, n)
    sp  <- regexpr(" [^ ]*$", cut)          # last space inside the cut
    if (sp > 1) cut <- substr(cut, 1, sp - 1)
    paste0(cut, intToUtf8(8230))            # ellipsis
  }, character(1), USE.NAMES = FALSE)
}

# R right-aligns every column of a data.frame and prints a <chr> type row, which
# reads badly for text. This renders text tables left-aligned, like the Python side.
show_table <- function(df) {
  esc <- function(x) { x <- as.character(x); x[is.na(x)] <- ""
                       x <- gsub("&", "&amp;", x, fixed = TRUE)
                       x <- gsub("<", "&lt;",  x, fixed = TRUE)
                       gsub(">", "&gt;", x, fixed = TRUE) }
  cell <- "text-align:left;padding:4px 18px 4px 0;white-space:nowrap"
  hd <- paste0("<th style='", cell, ";border-bottom:1px solid #999'>",
               esc(names(df)), "</th>", collapse = "")
  bd <- paste0(apply(df, 1, function(r)
         paste0("<tr>", paste0("<td style='", cell, "'>", esc(r), "</td>",
                collapse = ""), "</tr>")), collapse = "")
  IRdisplay::display_html(paste0(
    "<table style='border-collapse:collapse;font-size:13px;font-variant-numeric:tabular-nums'>",
    "<thead><tr>", hd, "</tr></thead><tbody>", bd, "</tbody></table>"))
}

cat("packages loaded ✓\n")

packages loaded ✓

Part A · Meet AQUAVIEW¶

Meet the tool¶

This is aquaview.org — the front door. One platform, one catalog, and an Explore map, an API, and AI search over it. Today we drive it from code.

aquaview.org homepage — 'One data plane from ingest to access', with the Explore map of global data density

Under the hood it's one catalog underneath, three ways to reach it, ~100 sources feeding in — today we learn the Code door, but the same catalog is behind all three:

AQUAVIEWone common catalog · 600,000+ datasets · ~100 sourcesExplorewebsite · point & clickCodePython · R (STAC API)AIMCP · plain Englishthe same catalog underneathNOAAIOOSArgoOBISmodelsAISbuoyssatellite

Five words we'll use¶

STAC (the standard AQUAVIEW speaks) nests like this — you only need these five:

AQUAVIEW
   └─ Source / Collection      e.g. "PacIOOS"          a data provider or product family
        └─ Dataset / Item      e.g. "ww3_hawaii"       one dataset record, with location + time
             └─ Assets         nc · csv · png · …       the actual downloadable files

A search returns items; each item lists its assets. That's the whole model.

A1 · Connect¶

AQUAVIEW speaks the STAC standard, so rstac just works — one line, no sign-in for discovery.

One host, defined once. Everything below uses STAC. If AQUAVIEW moves the endpoint, this is the only line you change.

In [2]:
STAC    <- "https://service.aquaview.org/stac"
catalog <- stac(STAC)

cat("Connected to AQUAVIEW ✓\n")
cat(STAC, "\n")

Connected to AQUAVIEW ✓

https://service.aquaview.org/stac

A2 · One catalog over everything¶

AQUAVIEW federates the whole ocean-data landscape into one searchable place. Here's the size, and a taste of the range.

In [3]:
# One request for the full source list, kept for the rest of the notebook.
collections <- (catalog |> collections() |> get_request())$collections
total <- (catalog |> stac_search(bbox = c(-180, -90, 180, 90), limit = 1) |>
            get_request())$numberMatched

cat(sprintf("%s datasets   ·   %d collections   ·   one catalog, one API\n\n",
            format(total, big.mark = ","), length(collections)))

kinds <- c(
  "Ship & float profiles"         = "WOD · Argo (GADR, RG_ARGO) · EN4 · CCHDO",
  "Satellite (SST, colour, SAR)"  = "CoastWatch · GOES-R · PolarWatch · Cerulean oil-slicks",
  "Gliders & autonomous"          = "IOOS Glider DAC · Spray · Voice of the Ocean",
  "Buoys & sensors"               = "NDBC · IOOS Sensors · CDIP waves",
  "Ocean & weather models"        = "HYCOM · RTOFS · GFS · Copernicus GLORYS",
  "Biodiversity & human activity" = "OBIS species · MarineCadastre AIS vessel traffic"
)
for (k in names(kinds)) cat(sprintf("  %-30s %s\n", k, kinds[[k]]))

608,355 datasets   ·   100 collections   ·   one catalog, one API


  Ship & float profiles          WOD · Argo (GADR, RG_ARGO) · EN4 · CCHDO
  Satellite (SST, colour, SAR)   CoastWatch · GOES-R · PolarWatch · Cerulean oil-slicks
  Gliders & autonomous           IOOS Glider DAC · Spray · Voice of the Ocean
  Buoys & sensors                NDBC · IOOS Sensors · CDIP waves
  Ocean & weather models         HYCOM · RTOFS · GFS · Copernicus GLORYS
  Biodiversity & human activity  OBIS species · MarineCadastre AIS vessel traffic

A3 · Find a source without memorising IDs¶

You don't have to know collection IDs in advance. Each collection carries a title, description and keywords; we already have all of them from the single request above, so we can filter locally — no extra calls, and you can see exactly what is being matched.

In [4]:
haystack <- function(c) {
  tolower(paste(c$id,
                c$title       %||% "",
                c$description %||% "",
                paste(unlist(c$keywords), collapse = " ")))
}
HAY <- vapply(collections, haystack, character(1))
IDS <- vapply(collections, function(c) c$id, character(1))

find_sources <- function(term, n = 6) {
  hit <- IDS[grepl(tolower(term), HAY, fixed = TRUE)]
  if (!length(hit)) "(none)" else head(hit, n)
}

for (term in c("glider", "sea surface temperature", "vessel", "argo"))
  cat(sprintf("  %-24s -> %s\n", term, paste(find_sources(term), collapse = ", ")))

  glider                   -> CALOOS, CCHDO, CENCOOS, IOOS, MARACOOS, NANOOS
  sea surface temperature  -> CDIP, COASTWATCH_WC, GOES_R, NEFSC, NOAA_GDP, POLARWATCH
  vessel                   -> CERULEAN, CERULEAN_ARCHIVE, CERULEAN_VESSEL_30D, MARINECADASTRE_AIS, OET_NA176
  argo                     -> CCHDO, GADR, ISAS, RG_ARGO

Read this result carefully — it is the most useful thing in Part A.

This searches how each source describes itself. That is a good first cut, but a source's blurb does not list every dataset inside it. find_sources("wave") misses PacIOOS, even though PacIOOS serves the Hawaiʻi WaveWatch III model we use later — PacIOOS simply doesn't say "wave" in its own description.

So the rule is:

You want… Use
a provider you half-remember find_sources(), above
datasets in your study area a bounding-box search — Part B, and it is the reliable one
to browse visually Explore

When the text search comes up short, search your box instead. That is what B4 does to find the glider.

A4 · Three ways in¶

The same catalog, whichever way you prefer to work — and a result is the same STAC item whichever door you use, so the data looks identical downstream:

Way in For This tutorial
Explore map at aquaview.org/explore click-and-draw discovery, no code shown in B8
STAC API (rstac / pystac-client) reproducible, scriptable the rest of this notebook
AI / MCP chat ask in plain English shown in B8

Part B · A worked example¶

Is the satellite right, along a Hawaiʻi glider track?¶

Your gliders record passive acoustics and CTD off Hawaiʻi. Before you can interpret that, you need the ocean around the track — and you need to know whether the satellite products you plan to lean on actually agree with what the glider measured.

We'll do it in six moves: search a box → narrow → count → find the glider → load its track → match the satellite to it.

B1 · Search one place, one window¶

The core move: a bounding box around the fleet and a date range. One call reaches across every source — no per-server hunting.

R note: rstac does not pad dates for you, so give the datetime as a full RFC 3339 timestamp (2025-01-01T00:00:00Z/..), not the bare 2025-01-01/... The raw API rejects the short form.

In [5]:
HAWAII <- c(-161, 18, -154, 23)          # west, south, east, north
WHEN   <- "2025-01-01T00:00:00Z/.."      # 2025 to now

search <- catalog |> stac_search(bbox = HAWAII, datetime = WHEN, limit = 100) |> get_request()
cat(sprintf("%s datasets match in the Hawaiʻi box (2025–present)\n",
            format(search$numberMatched, big.mark = ",")))

4,228 datasets match in the Hawaiʻi box (2025–present)

B2 · From a broad catch to what you actually want¶

A wide search returns everything that intersects your box and time — that's discovery working, not failing.

Two things to notice in the table below.

  1. The first page is not a ranking. Results come back in the catalog's own order, not by relevance, so the top rows tell you nothing about the mix. (Cerulean SAR oil-slick detections lead here; they are about a sixth of the box, not the bulk of it.)
  2. Many dates are a range, not a moment. A continuously-updated model or a glider deployment has start_datetime/end_datetime and a null datetime. Reading only datetime shows a blank — so we fall back.
In [6]:
when_of <- function(p) {
  # STAC items are either a moment or a range — show whichever this one has.
  if (!is.null(p$datetime)) return(substr(p$datetime, 1, 10))
  s <- substr(p$start_datetime %||% "", 1, 10)
  e <- substr(p$end_datetime   %||% "", 1, 10)
  if (nzchar(s) || nzchar(e)) paste(s, "->", e) else ""
}

to_row <- function(f) data.frame(
  source = f$collection,
  id     = short(f$id, 26),
  when   = when_of(f$properties),
  title  = short(f$properties$title %||% "", 44),
  stringsAsFactors = FALSE
)

show_table(do.call(rbind, lapply(search$features[1:6], to_row)))   # the raw, unfiltered catch
sourceidwhentitle
CERULEAN52142482026-09-12Oil slick 5214248
CERULEAN52142472026-09-12Oil slick 5214247
CERULEAN52142462026-09-12Oil slick 5214246
CERULEAN52142452026-09-12Oil slick 5214245
CERULEAN52142952026-09-12Oil slick 5214295
CERULEAN52142942026-09-12Oil slick 5214294

So you narrow — the real discovery skill. Restrict to the sources relevant to a glider mission and pull one representative dataset from each.

Note the SST source is COASTWATCH, not COASTWATCH_WC. COASTWATCH_WC is the West Coast node — over Hawaiʻi it carries HF-radar currents and ocean colour. The satellite SST products live in plain COASTWATCH, and that is the one we validate against in B6.

In [7]:
relevant <- c("COASTWATCH", "PacIOOS", "OBIS", "MARINECADASTRE_AIS", "IOOS")

focused <- catalog |> stac_search(bbox = HAWAII, collections = relevant,
                                  datetime = WHEN, limit = 1) |> get_request()
cat(sprintf("%s datasets across the sources we care about\n\n",
            format(focused$numberMatched, big.mark = ",")))

first_of <- function(cid) {
  r <- catalog |> stac_search(collections = cid, bbox = HAWAII, datetime = WHEN, limit = 1) |>
       get_request()
  if (length(r$features)) to_row(r$features[[1]]) else NULL
}
show_table(do.call(rbind, lapply(relevant, first_of)))   # identical columns, whatever the source

229 datasets across the sources we care about
sourceidwhentitle
COASTWATCHnoaacweddymunsterdaily2020-05-05 -> 2026-06-08MUltiparameter NRT System for Tracking…
PacIOOSww3_hawaii_lon1802011-06-21 -> 2026-09-19WaveWatch III (WW3) Hawaii Regional Wave…
OBIS345c1b88-2bab-4403-ad2c-17…2026-06-29DFO Maritimes Region Cetacean Sightings
MARINECADASTRE_AISmc_dataset_679192024-07-01 -> 2026-03-01Vessel Routing Measures
IOOSsg626-20250729T14522025-07-29 -> 2025-10-29sg626-20250729T1452

B3 · What's here, by source?¶

Before downloading anything, ask how much of each kind exists in your box — one cheap count per source, no data moved. numberMatched returns the total without fetching a single file.

In [8]:
of_interest <- c(
  COASTWATCH         = "Satellite SST",
  COASTWATCH_WC      = "Ocean colour / HF radar",
  PacIOOS            = "Hawaiʻi waves & models",
  OBIS               = "Species-occurrence context",
  MARINECADASTRE_AIS = "Vessel-traffic context",
  IOOS               = "Glider Data Assembly Center"
)
count_source <- function(cid)
  (catalog |> stac_search(collections = cid, bbox = HAWAII, datetime = WHEN, limit = 1) |>
     get_request())$numberMatched

counts <- vapply(names(of_interest), count_source, numeric(1))
df <- data.frame(label = unname(of_interest[names(counts)]), n = as.integer(counts))
df$label <- factor(df$label, levels = df$label[order(df$n)])

options(repr.plot.width = 8, repr.plot.height = 3.4)
ggplot(df, aes(n, label)) +
  geom_col(fill = "#0b7d84") +
  geom_text(aes(label = format(n, big.mark = ",")), hjust = -0.15, size = 3.2, colour = "#334455") +
  labs(title = "Datasets in the Hawaiʻi glider box, by source (2025–present)",
       x = "datasets", y = NULL) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.14))) +
  theme_minimal(base_size = 11) +
  theme(panel.grid.major.y = element_blank(), plot.title = element_text(size = 11))
No description has been provided for this image

B4 · Find the glider — the way you'd actually find it¶

We don't know any glider IDs. We know a box and a time. That is enough: search the IOOS Glider Data Assembly Center inside the Hawaiʻi box and read what comes back.

In [9]:
gliders <- catalog |> stac_search(bbox = HAWAII, collections = "IOOS",
                                  datetime = WHEN, limit = 50) |> get_request()

cat(sprintf("%d glider deployments in the box since 2025\n\n", length(gliders$features)))
for (f in gliders$features)
  cat(sprintf("  %-24s %s -> %s\n", f$id,
              substr(f$properties$start_datetime, 1, 10),
              substr(f$properties$end_datetime,   1, 10)))

6 glider deployments in the box since 2025


  sg626-20250729T1452      2025-07-29 -> 2025-10-29
  sg511-20260722T1831      2026-07-22 -> 2026-09-14
  sg511-20250729T1532      2025-07-29 -> 2025-09-09
  sg148-20260722T1939      2026-07-22 -> 2026-09-14
  ng457-20260713T0000      2026-07-23 -> 2026-09-14
  ng344-20260713T0000      2026-07-23 -> 2026-09-14

We'll use sg626-20250729T1452 — a completed three-month Seaglider deployment, so the numbers in this notebook stay reproducible. (The 2026 deployments are still flying; their data grows every day. Swapping one in is the exercise at the end.)

Every item carries its variable list, so you can check it has what you need before downloading.

In [10]:
GLIDER <- "sg626-20250729T1452"
g <- catalog |> collections("IOOS") |> items(GLIDER) |> get_request()

cat(g$properties$title %||% g$id, "\n")
cat("period :", substr(g$properties$start_datetime, 1, 10), "->",
                substr(g$properties$end_datetime,   1, 10), "\n")

v <- unlist(g$properties$`aquaview:variables`)
cat(sprintf("\n%d variables. The ones we want:\n", length(v)))
for (w in c("time", "latitude", "longitude", "depth", "temperature", "salinity"))
  cat(sprintf("   %-12s %s\n", w, if (w %in% v) "yes" else "NO"))

cat("\nOther instruments on board:",
    paste(grep("^(sbe43|aa5013|bb2f)", v, value = TRUE), collapse = ", "), "\n")
cat("\nTabular assets:",
    paste(intersect(c("csv", "nc", "json"), names(g$assets)), collapse = ", "), "\n")

sg626-20250729T1452 

period : 2025-07-29 -> 2025-10-29 

60 variables. The ones we want:

   time         yes
   latitude     yes
   longitude    yes
   depth        yes
   temperature  yes
   salinity     yes

Other instruments on board: aa5013_o2, bb2f_470, bb2f_695, bb2f_700, sbe43_o2 

Tabular assets: csv, nc, json

B5 · Load the track¶

The csv asset points at the IOOS Glider DAC's own ERDDAP. AQUAVIEW standardises discovery and metadata; it does not hide the source — so we append ERDDAP's native tabledap query to pick columns and a time window, and only that slice is transferred.

Always constrain the time window. Without it you would ask for the entire three-month deployment.

ERDDAP puts the units on row 2, so we skip it and name the columns ourselves.

In [11]:
T0 <- "2025-08-01T00:00:00Z"; T1 <- "2025-08-15T00:00:00Z"

url <- paste0(g$assets$csv$href,
              "?time,latitude,longitude,depth,temperature,salinity",
              "&time>=", T0, "&time<=", T1)

track <- read_csv(url, skip = 2, show_col_types = FALSE,
                  col_names = c("time", "latitude", "longitude",
                                "depth", "temperature", "salinity"),
                  col_types = cols(time = col_datetime(format = "%Y-%m-%dT%H:%M:%SZ"),
                                   .default = col_double())) |>
         filter(!is.na(latitude), !is.na(longitude), !is.na(temperature))

cat(sprintf("%s measurements over %d days\n", format(nrow(track), big.mark = ","),
            as.integer(difftime(max(track$time), min(track$time), units = "days"))))
cat(sprintf("  latitude    %.3f -> %.3f\n", min(track$latitude),  max(track$latitude)))
cat(sprintf("  longitude   %.3f -> %.3f\n", min(track$longitude), max(track$longitude)))
cat(sprintf("  depth       %.1f -> %.1f m\n", min(track$depth),   max(track$depth)))
cat(sprintf("  temperature %.2f -> %.2f C\n", min(track$temperature), max(track$temperature)))
head(as.data.frame(track))

84,176 measurements over 13 days

  latitude    18.827 -> 20.994

  longitude   -158.135 -> -156.176

  depth       -0.2 -> 907.6 m

  temperature 4.61 -> 29.38 C
A data.frame: 6 × 6
timelatitudelongitudedepthtemperaturesalinity
<dttm><dbl><dbl><dbl><dbl><dbl>
12025-08-01 04:35:2120.9936-158.13450.6027.390NaN
22025-08-01 04:35:2120.9936-158.13450.7627.373NaN
32025-08-01 04:35:2120.9936-158.13450.7227.368NaN
42025-08-01 04:35:2120.9936-158.13450.6727.358NaN
52025-08-01 04:35:2120.9936-158.13451.0127.393NaN
62025-08-01 04:35:2120.9936-158.13451.4127.386NaN
In [12]:
options(repr.plot.width = 6, repr.plot.height = 4.6)
ggplot(track, aes(longitude, latitude, colour = as.numeric(time))) +
  geom_point(size = 0.35) +
  annotate("point", x = track$longitude[1], y = track$latitude[1],
           shape = 21, size = 3.5, colour = "#b5642f", stroke = 1.1) +
  annotate("text", x = track$longitude[1], y = track$latitude[1],
           label = " start", hjust = 0, size = 3.3, colour = "#b5642f") +
  scale_colour_viridis_c(guide = "none") +
  coord_fixed() +
  labs(title = paste0(GLIDER, " — track, 1–15 Aug 2025"), x = "longitude", y = "latitude") +
  theme_minimal(base_size = 11)
No description has been provided for this image
In [13]:
options(repr.plot.width = 9, repr.plot.height = 4.2)
ggplot(track, aes(time, depth, colour = temperature)) +
  geom_point(size = 0.35) +
  scale_y_reverse() +
  scale_colour_distiller(palette = "RdYlBu", name = "°C") +
  labs(title = "Temperature section", x = NULL, y = "depth (m)") +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))
No description has been provided for this image

A two-week southbound transect from off Oʻahu down past Maui, diving to ~900 m. The section shows the warm mixed layer sitting on a sharp thermocline — exactly the structure that decides how sound travels, which is why it matters for passive acoustics.

B6 · Now check the satellite against it¶

This is the payoff, and it is the same comparison the CoastWatch glider tutorial makes — except we found both datasets through one catalog.

Step 1 — the satellite product. We know from B3 that SST lives in COASTWATCH. Find the blended analysis by searching inside our box rather than guessing an ID.

In [14]:
sst_hits <- catalog |> stac_search(bbox = HAWAII, collections = "COASTWATCH",
                                   datetime = WHEN, limit = 100) |> get_request()
for (f in sst_hits$features) {
  t <- f$properties$title %||% ""
  if (grepl("Blended", t) && grepl("Sea-Surface Temperature", t))
    cat(sprintf("%-26s %s\n", f$id, substr(t, 1, 64)))
}

noaacwBLENDEDsstDaily      Sea-Surface Temperature, NOAA Geo-polar Blended Analysis Night O
noaacwBLENDEDsstDNDaily    Sea-Surface Temperature, NOAA Geo-polar Blended Analysis Day+Nig
noaacwBLENDEDsstDLDaily    Sea-Surface Temperature, NOAA Geo-polar Blended Analysis Diurnal
noaacwBLENDEDCsstDaily     Sea-Surface Temperature, NOAA Geo-polar Blended Analysis Night O

noaacwBLENDEDsstDaily is NOAA's Geo-polar Blended daily SST — geostationary plus polar-orbiting satellites merged onto a 5 km grid, 2002 to present. It is the night-only product, which is the right choice here: daytime skin warming makes the satellite look warmer than the top ten metres a glider actually samples, and night values compare far more fairly. (...DNDaily is the day+night version if you want to see that effect for yourself.)

Step 2 — open it lazily. The nc asset is a materialising ERDDAP endpoint: requesting it with no constraints asks for the whole 24-year record and ERDDAP refuses with 413 Payload Too Large. Strip the .nc extension and you get the griddap OPeNDAP endpoint, which ncdf4 opens lazily — metadata only, nothing transferred until you ask for values.

In [15]:
item    <- catalog |> collections("COASTWATCH") |> items("noaacwBLENDEDsstDaily") |> get_request()
opendap <- sub("\\.nc$", "", item$assets$nc$href)      # <- the one change that matters

nc <- nc_open(opendap)                                 # lazy: metadata only
cat(opendap, "\n")
cat("  grid :", paste(sapply(nc$dim, function(d) paste0(d$name, "=", d$len)), collapse = "  "), "\n")
cat("  units:", ncatt_get(nc, "analysed_sst", "units")$value, "\n")

https://coastwatch.noaa.gov/erddap/griddap/noaacwBLENDEDsstDaily 

  grid : latitude=3600  longitude=7200  time=8720 

  units: kelvin

Step 3 — subset to the track, then read. ncdf4 reads by index, so we find the index ranges that cover the glider's own bounding box (plus a small pad) and its two-week window, then read only that hyperslab. That turns a 24-year global grid into a few hundred kilobytes.

In [16]:
PAD <- 0.25

lat_all <- ncvar_get(nc, "latitude")
lon_all <- ncvar_get(nc, "longitude")
tim_all <- as.Date(as.POSIXct(ncvar_get(nc, "time"), origin = "1970-01-01", tz = "UTC"))

i <- which(lat_all >= min(track$latitude)  - PAD & lat_all <= max(track$latitude)  + PAD)
j <- which(lon_all >= min(track$longitude) - PAD & lon_all <= max(track$longitude) + PAD)
k <- which(tim_all >= as.Date(substr(T0, 1, 10)) & tim_all <= as.Date(substr(T1, 1, 10)))

sst <- ncvar_get(nc, "analysed_sst",                   # <- the only network transfer
                 start = c(min(j), min(i), min(k)),
                 count = c(length(j), length(i), length(k))) - 273.15   # grid is in kelvin
nc_close(nc)

sub_lon <- lon_all[j]; sub_lat <- lat_all[i]; sub_tim <- tim_all[k]
cat(sprintf("subset: lon %d x lat %d x time %d  ->  %.0f kB\n",
            length(j), length(i), length(k), object.size(sst) / 1024))

subset: lon 49 x lat 53 x time 15  ->  305 kB

Step 4 — put them on the same footing. The glider samples continuously down to 900 m; the satellite sees one skin value per day. So we take the glider's shallowest 10 m and average it per day — the same choice the CoastWatch tutorial makes.

In [17]:
daily <- track |>
  filter(depth <= 10) |>
  mutate(day = as.Date(time)) |>
  group_by(day) |>
  summarise(lat = mean(latitude), lon = mean(longitude),
            glider_sst = mean(temperature), n = n(), .groups = "drop")

cat(sprintf("%s shallow measurements -> %d daily surface values\n",
            format(sum(track$depth <= 10), big.mark = ","), nrow(daily)))
head(as.data.frame(daily))

1,942 shallow measurements -> 14 daily surface values
A data.frame: 6 × 5
daylatlonglider_sstn
<date><dbl><dbl><dbl><int>
12025-08-0120.95329-158.091327.33631140
22025-08-0220.84741-157.983227.36510143
32025-08-0320.72226-157.874027.49277133
42025-08-0420.59244-157.744227.45597134
52025-08-0520.43775-157.611927.34807141
62025-08-0620.28226-157.495626.99283124

Step 5 — match each day to the satellite. For every glider day, take the satellite pixels within ±0.1° of where the glider actually was and average them. This is the match-up.

In [18]:
satellite_at <- function(day, lat, lon, halfwidth = 0.1) {
  ti <- which.min(abs(as.numeric(sub_tim - day)))                 # nearest day
  jj <- which(sub_lon >= lon - halfwidth & sub_lon <= lon + halfwidth)
  ii <- which(sub_lat >= lat - halfwidth & sub_lat <= lat + halfwidth)
  mean(sst[jj, ii, ti], na.rm = TRUE)
}

daily$sat_sst <- mapply(satellite_at, daily$day, daily$lat, daily$lon)
daily$diff    <- daily$glider_sst - daily$sat_sst

bias <- mean(daily$diff)
rmse <- sqrt(mean(daily$diff^2))
corr <- cor(daily$glider_sst, daily$sat_sst)

cat(sprintf("  mean bias    %+.3f C   (glider minus satellite)\n", bias))
cat(sprintf("  RMSE         %.3f C\n", rmse))
cat(sprintf("  correlation  %.3f\n", corr))
cat(sprintf("  n            %d days\n", nrow(daily)))
as.data.frame(daily |> select(day, lat, lon, glider_sst, sat_sst, diff))

  mean bias    -0.039 C   (glider minus satellite)

  RMSE         0.111 C

  correlation  0.941

  n            14 days
A data.frame: 14 × 6
daylatlonglider_sstsat_sstdiff
<date><dbl><dbl><dbl><dbl><dbl>
2025-08-0120.95329-158.091327.3363127.35813-0.021811019
2025-08-0220.84741-157.983227.3651027.30250 0.062598512
2025-08-0320.72226-157.874027.4927727.48438 0.008386577
2025-08-0420.59244-157.744227.4559727.38625 0.069722896
2025-08-0520.43775-157.611927.3480727.32500 0.023072448
2025-08-0620.28226-157.495626.9928326.97000 0.022833621
2025-08-0720.13718-157.388626.8571227.00813-0.151009857
2025-08-0819.98241-157.262726.8243727.00063-0.176252680
2025-08-0919.79673-157.149026.9381527.03250-0.094345647
2025-08-1019.60806-157.021427.0548126.96625 0.088562604
2025-08-1119.43096-156.870327.3629627.28000 0.082964565
2025-08-1219.22049-156.688427.1976927.26000-0.062315623
2025-08-1319.04210-156.462926.7265226.96000-0.233475367
2025-08-1418.88933-156.238026.8399827.00250-0.162523590
In [19]:
options(repr.plot.width = 8, repr.plot.height = 4)
long <- rbind(
  data.frame(day = daily$day, sst = daily$glider_sst,
             src = paste0("Glider ", sub("-.*", "", GLIDER), " (<=10 m)")),
  data.frame(day = daily$day, sst = daily$sat_sst,
             src = "Satellite (Geo-polar Blended)")
)
ggplot(long, aes(day, sst, colour = src, shape = src, linetype = src)) +
  geom_line(linewidth = 0.7) + geom_point(size = 2.1) +
  scale_colour_manual(values = c("#0b7d84", "#b5642f"), name = NULL) +
  scale_shape_manual(values = c(16, 15), name = NULL) +
  scale_linetype_manual(values = c("solid", "dashed"), name = NULL) +
  labs(title = "Glider vs satellite along the track",
       x = NULL, y = "sea-surface temperature (°C)") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top", axis.text.x = element_text(angle = 30, hjust = 1))
No description has been provided for this image
In [20]:
options(repr.plot.width = 5, repr.plot.height = 4.6)
lims <- range(c(daily$glider_sst, daily$sat_sst)) + c(-0.15, 0.15)
ggplot(daily, aes(sat_sst, glider_sst)) +
  geom_abline(slope = 1, intercept = 0, colour = "#c9c2b4") +
  geom_point(colour = "#0b7d84", size = 2.6) +
  coord_fixed(xlim = lims, ylim = lims) +
  labs(title = sprintf("bias %+.2f · RMSE %.2f · r %.2f", bias, rmse, corr),
       x = "satellite (°C)", y = "glider (°C)") +
  theme_minimal(base_size = 11)
No description has been provided for this image

That is a real validation result, reproduced from a standing start in one notebook: the blended satellite product tracks this glider to about a tenth of a degree, with essentially no bias, over two weeks and 250 km of ocean. Where the two part company — around 7–8 and 13–14 August — the glider is cooler than the satellite, which is what you would expect when the wind mixes the surface and the skin value the satellite sees stops representing the top ten metres.

(The Python notebook produces these same three numbers, to the third decimal.)

B7 · The layers a single-server tutorial can't add¶

Everything so far you could have done against one ERDDAP server. This is the part you cannot: the same box and window, asked of sources that have nothing to do with each other.

In [21]:
track_box <- c(min(track$longitude) - PAD, min(track$latitude) - PAD,
               max(track$longitude) + PAD, max(track$latitude) + PAD)
cat("track box:", paste(round(track_box, 2), collapse = ", "), "\n\n")

context <- c(
  COASTWATCH         = "satellite SST / ocean colour",
  PacIOOS            = "wave + circulation models",
  MARINECADASTRE_AIS = "vessel traffic (noise context for PAM)",
  OBIS               = "species occurrence datasets",
  CERULEAN           = "SAR oil-slick detections",
  IOOS               = "other gliders in the same water"
)
for (cid in names(context)) {
  n <- (catalog |> stac_search(collections = cid, bbox = track_box,
                               datetime = WHEN, limit = 1) |> get_request())$numberMatched
  cat(sprintf("  %5s  %-20s %s\n", format(n, big.mark = ","), cid, context[[cid]]))
}

track box: -158.38, 18.58, -155.93, 21.24 


     89  COASTWATCH           satellite SST / ocean colour
     51  PacIOOS              wave + circulation models
     26  MARINECADASTRE_AIS   vessel traffic (noise context for PAM)
     39  OBIS                 species occurrence datasets
    235  CERULEAN             SAR oil-slick detections
      6  IOOS                 other gliders in the same water

One box, six providers, one loop. For a PAM project the AIS row is the interesting one — vessel traffic is the dominant source of low-frequency noise, and it is discoverable here next to the oceanography rather than on a different agency's website.

Honest caveat, worth knowing on the day. The AIS records are in the catalog, but MarineCadastre publishes the underlying daily files about a year in arrears, and the 2025/2026 asset links in the catalog do not resolve yet. Treat AIS here as discovery — the catalog tells you what exists and when — and fetch 2023/2024 files directly from MarineCadastre if you need vessel data this week.

In [22]:
ais <- catalog |> stac_search(collections = "MARINECADASTRE_AIS", bbox = track_box,
                              datetime = "2025-08-01T00:00:00Z/2025-08-31T00:00:00Z",
                              limit = 20) |> get_request()

for (f in ais$features)
  cat(sprintf("%-40s %3d assets\n", f$properties$title %||% f$id, length(f$assets)))

monthly <- Filter(function(f) length(f$assets) > 0, ais$features)[[1]]
cat("\n", monthly$id, "->",
    paste(head(sort(names(monthly$assets)), 5), collapse = ", "), "...\n")
cat(" example href:", monthly$assets$day_04$href, "\n")

AIS Vessel Transit Counts 2025             0 assets
AIS Vessel Tracks 2025                     0 assets
AIS Vessel Traffic -- August 2025         32 assets

 marinecadastre_ais_2025_08 -> data_browse, day_01, day_02, day_03, day_04 ...

 example href: https://coast.noaa.gov/htdata/CMSP/AISDataHandler/2025/ais-2025-08-04.csv.zst

And the wave model, while we're here¶

The Hawaiʻi WaveWatch III run is in PacIOOS. Two things to get right, and they are the two mistakes everyone makes with ERDDAP:

  1. Pick the right longitude convention. PacIOOS publishes each model twice. ww3_hawaii runs 0–360° (199→206 here); ww3_hawaii_lon180 runs −180–180° and matches our box exactly.
  2. Give every axis a constraint. Thgt is [time][depth][latitude][longitude] — four axes. Supply three and ERDDAP silently reads your latitude range as depth and returns an error image.
In [23]:
wave <- catalog |> collections("PacIOOS") |> items("ww3_hawaii_lon180") |> get_request()
cat(wave$properties$title, "\n")
cat("variables:", paste(unlist(wave$properties$`aquaview:variables`), collapse = ", "), "\n")

graph <- paste0(wave$assets$png$href,                  # Thgt[time][depth][lat][lon]
                "?Thgt[(last)][(0.0)][(18.0):(23.0)][(-161.0):(-154.0)]",
                "&.draw=surface&.vars=longitude|latitude|Thgt&.land=under")

tmp <- file.path(tempdir(), "wave.png")
h   <- curl::curl_fetch_disk(graph, tmp)

# ERDDAP reports failure in the *body* with HTTP 200 - check the type, not the status code.
if (grepl("image", h$type)) {
  IRdisplay::display_png(file = tmp)
} else {
  cat("ERDDAP said:", substr(readLines(tmp, warn = FALSE)[1], 1, 300), "\n")
}

WaveWatch III (WW3) Hawaii Regional Wave Model, Lon+/-180 

variables: Tdir, Tper, Thgt, sdir, sper, shgt, wdir, wper, whgt
No description has been provided for this image

B8 · On the map, and by just asking¶

Prefer clicking? The same box and sources on the live Explore map. The screenshot links to the real, preconfigured view — note that Explore counts what is in the current viewport, so its number differs from numberMatched above, which counts the whole box.

AQUAVIEW Explore — Hawaiʻi box, 5 sources, hex-binned

Prefer asking? AQUAVIEW's MCP server is live at https://mcp.aquaview.org/mcp. Point any MCP client at it — in Claude Desktop, add:

{ "mcpServers": { "aquaview": { "url": "https://mcp.aquaview.org/mcp" } } }

then ask "what SST and vessel-traffic data is near 20°N 157°W in August 2025?" and you get the same STAC items you just found by hand. Same catalog, third door.

Your turn¶

You now have the whole loop: meet the catalog → search a place and time → narrow → count → find an instrument → load it → validate against it.

Change one thing and re-run:

  1. Another glider. Swap GLIDER for sg511-20260722T1831 — it is flying right now. You'll need to move T0/T1 into its window (B4 prints the dates). Does the satellite agree as well?
  2. Another variable. salinity is already in the track data frame. COASTWATCH also carries ocean colour — try matching chlorophyll instead of SST.
  3. Another place. Change HAWAII to a box over your own study area and re-run B1–B4. Which sources show up that you didn't expect?
  4. Tighten the match-up. We used a ±0.1° box and a daily mean. Does the RMSE change if you use the single nearest pixel instead?

For a hackathon project, the sources to reach for:

Question Source
Ocean state along the track COASTWATCH (SST, colour), PacIOOS (waves, currents)
What vessel traffic is nearby? MARINECADASTRE_AIS (discovery; files lag ~1 year)
Which species are recorded here? OBIS
Sound-speed / CTD context the glider's own temperature + salinity, WOD, GADR
Other gliders in the same water IOOS

Free account, STAC access included → aquaview.org