R

Every R package that reads rasters through GDAL can read these COGs from a /vsicurl/ URL, fetching only the tiles it needs. See The data for the grid, units and status fields.

  1. A URL for a date
  2. terra
  3. gdalraster
  4. vapour
  5. stars
  6. The index
  7. Many days
  8. raadtools
  9. Virtual stores

A URL for a date

mur_url <- function(date, var = "analysed_sst") {
  date <- as.Date(date)
  sprintf(
    "/vsicurl/https://data.source.coop/ausantarctic/ghrsst-mur-v2/%s/%s090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1_%s.tif",
    format(date, "%Y/%m/%d"), format(date, "%Y%m%d"), var
  )
}
dsn <- mur_url("2025-06-01")

Let GDAL skip listing the remote folder when it opens a file (set this before reading):

terra::setGDALconfig("GDAL_DISABLE_READDIR_ON_OPEN", "EMPTY_DIR")

terra

terra applies the COG’s scale and offset, so analysed_sst comes back in degrees Celsius:

library(terra)
sst <- rast(dsn)
sst

# A region: only the tiles that cover it are read
crop(sst, ext(130, 150, -55, -42))

# A coarse whole-world view: read from an overview
plot(sst, maxcell = 1e6)

Reprojecting with by_util = TRUE uses GDAL’s warper directly, which reads the overview that matches the output resolution, instead of the full grid:

target <- rast(
  ext(c(-1, 1, -1, 1) * 6e5), res = 1000,
  crs = "+proj=laea +lon_0=147 +lat_0=-45"
)
project(sst, target, by_util = TRUE)

gdalraster

gdalraster exposes GDAL’s API closely:

library(gdalraster)
ds <- new(GDALRaster, dsn, TRUE)
ds$info()
ds$getScale(1); ds$getOffset(1)   # 0.001 and 25: degrees Celsius
ds$close()

# Warp a region to a local projection, written to a file
warp(dsn, tf <- tempfile(fileext = ".tif"),
     t_srs = "+proj=laea +lon_0=147 +lat_0=-45",
     cl_arg = c("-te", -6e5, -6e5, 6e5, 6e5, "-tr", 1000, 1000, "-r", "bilinear"))

vapour

vapour reads straight to R vectors through the warper, for any extent, size and CRS in one call. Give one dimension and 0 for the other so GDAL derives it from the extent:

library(vapour)
v <- gdal_raster_data(
  dsn,
  target_ext = c(-6e5, 6e5, -6e5, 6e5),
  target_dim = c(1024, 0),
  target_crs = "+proj=laea +lon_0=147 +lat_0=-45",
  resample = "bilinear"
)
str(v)

The values come back as stored; apply the scale and offset shown by gdal raster info (see Values and units).

stars

library(stars)
sstars <- read_stars(dsn, proxy = TRUE)
sstars[st_bbox(c(xmin = 130, xmax = 150, ymin = -55, ymax = -42), crs = st_crs(4326))]

The index

The STAC GeoParquet index lists every day, its status and its asset URLs. With arrow:

library(arrow)
library(dplyr)

catalog <- read_parquet(
  "https://data.source.coop/ausantarctic/ghrsst-mur-v2/ghrsst-mur-v2.parquet",
  col_select = c("datetime", "ghrsst:status", "assets")
)
count(catalog, `ghrsst:status`)

final <- catalog |>
  filter(`ghrsst:status` == "final") |>
  mutate(href = assets$analysed_sst$href)

Or with DuckDB, from R:

library(duckdb)
con <- dbConnect(duckdb())
dbGetQuery(con, "
  SELECT datetime, assets.analysed_sst.href AS href
  FROM 'https://data.source.coop/ausantarctic/ghrsst-mur-v2/ghrsst-mur-v2.parquet'
  WHERE \"ghrsst:status\" = 'final' AND datetime >= '2024-01-01'
  ORDER BY datetime")

A single day’s STAC item is plain JSON:

item <- jsonlite::fromJSON(
  "https://data.source.coop/ausantarctic/ghrsst-mur-v2/2025/06/01/20250601090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.stac-item.json"
)
item$properties$`ghrsst:status`
item$assets$analysed_sst$href

Many days

terra opens many COGs as the layers of one SpatRaster without reading pixels, then reads only the tiles each operation touches. A year of SST at a point, from final days:

dates <- seq(as.Date("2024-01-01"), as.Date("2024-12-31"), by = "1 day")
cube <- rast(mur_url(dates))
time(cube) <- dates

series <- extract(cube, cbind(147, -45), raw = TRUE)

Opening each file reads its header, so a year is a few hundred small requests before any pixels. Filter dates to final days from the index first if you need only finals.

raadtools

At the Australian Antarctic Division, raadfiles::ghrsst_daily_files() lists these same COG URLs, and raadtools::readghrsst() reads them. See Background.

Virtual stores

A GDAL multidimensional VRT over the daily COGs, ghrsst-mur-v2.vrt, opens the archive as one time series array per variable, with no per-file step: terra::rast(x, md = TRUE), or gdalraster’s multidimensional functions. See Virtual stores.


Back to top

Code: Apache-2.0 licence (see the repository). Data: GHRSST MUR v4.1 from NASA JPL, distributed by PO.DAAC, repackaged without changing values.

This site uses Just the Docs, a documentation theme for Jekyll.