Extracting SST and spatial data • dhw

library(dhw)
library(sf)
library(tidyverse)
library(terra)

Extract GBR spatial data

Use this function to download the GBR Reefs dataset via eAtlas. Several versions of the GBRMPA shp file exist, this version from eAtlas include the Torres Straits:

This dataset consists of a shapefile of the reefs, islands, sand banks, cays and rocks of the whole Great Barrier Reef (GBR) including Torres Strait. This dataset is an extension of the mapping in the GBR Marine Park to include Torres Strait. The Torres Strait region was mapped at a scale of 1:50,000 (Lawrey, E. P., Stewart M., 2016) and these new features are referred to as the “Torres Strait Reef and Island Features” dataset. The Complete GBR Reef and Island Features dataset integrates the “Torres Strait Reef and Island Features” dataset with the existing “GBR Features” (Great Barrier Reef Marine Park Authority, 2007) to create a single composite dataset of the whole Great Barrier Reef.

for metadata see here:

Complete Great Barrier Reef (GBR) Island and Reef Feature boundaries including Torres Strait Version 1b (NESP TWQ 3.13, AIMS, TSRA, GBRMPA)

The function returns either “hull” (convex hull around reefs), “base” (basic with no changes), or “combined” (polygons merged to reef ID)


gbr_reefs_combined <- download_gbr_spatial(return="combined")
gbr_reefs <- download_gbr_spatial(return="base")

gbr_reefs

GBR_hull <- download_gbr_spatial(return="hull", crs = "EPSG:7844")

Set up a boundary to mask the raster data by buffering a 1km concave hull surrounding the GBR Reefs:


gbr_reefs_border <- gbr_reefs |> 
  st_make_valid() |> 
  dplyr::filter(FEAT_NAME %in% c("Reef", "Terrestrial Reef")) 

gbr_reefs_hull <- gbr_reefs_border |> 
  concaveman::concaveman() |> 
  st_make_valid()

gbr_reefs_hull_buffered <- gbr_reefs_hull |> 
  st_buffer(1000) 

ggplot() + theme_bw() +
  geom_sf(data=gbr_reefs_hull_buffered,  fill="turquoise3") +
  geom_sf(data=gbr_reefs_hull, fill="lightblue") +
  geom_sf(data=gbr_reefs_border, fill="black")

Downloading and extracting SST data

NOAA provide many options for downloading SST data. ERDAPP via the rerdapp library is very useful for smaller requests, but (in my experience) downloading long-time series or large spatial extents of data can result in signficant timeouts.

Below are three helper functions that directly download CRW, OISST, and ERA5 datasets at global scales:

Extract NOAA CoralTemp spatial data

download_CoralTemp() pulls global data from the ncei https as follows:




download_CoralTemp(url = "https://www.ncei.noaa.gov/thredds-ocean/fileServer/crw/5km/v3.1/nc/v1.0/daily/",
                   start_date = "2025-01-01",
                   end_date = "2025-03-31",
                   dest_dir = "/Volumes/Extreme_SSD/dhw/CRW/CRW_SST/",
                   variable = "sst",
                   mc.cores = 10)

download_CoralTemp(url = "https://www.ncei.noaa.gov/thredds-ocean/fileServer/crw/5km/v3.1/nc/v1.0/daily/",
                   start_date = "2025-01-01",
                   end_date = "2025-03-31",
                   dest_dir = "/Volumes/Extreme_SSD/dhw/CRW/CRW_SSTA/",
                   variable = "ssta",
                   quiet = FALSE,
                   mc.cores = 10)

download_CoralTemp(url = "https://www.ncei.noaa.gov/thredds-ocean/fileServer/crw/5km/v3.1/nc/v1.0/daily/",
                   start_date = "2025-01-01",
                   end_date = "2025-03-31",
                   dest_dir = "/Volumes/Extreme_SSD/dhw/CRW/CRW_HS/",
                   variable = "hs",
                   mc.cores = 10)

download_CoralTemp(url = "https://www.ncei.noaa.gov/thredds-ocean/fileServer/crw/5km/v3.1/nc/v1.0/daily/",
                   start_date = "2025-01-01",
                   end_date = "2025-03-31",
                   dest_dir = "/Volumes/Extreme_SSD/dhw/CRW/CRW_DHW/",
                   variable = "dhw",
                   quiet = FALSE,
                   mc.cores = 10)

the full global dataset for four variables (sst, ssta, hs, dhw) is 520.78GB in individual ~10-11MB .nc files.

process_CoralTemp does three functions: 1) crop-mask to a polygon shape (here the GBR hull), 2) combine rast files using c(), 3) output to an .rds file for saving. The function will be paralellised (later) as it is slow


process_CoralTemp(input = "/Volumes/Extreme_SSD/dhw/CRW/CRW_SST",
                  polygon = GBR_hull, crs = "EPSG:7844",
                  combinedfilename = "/Volumes/Extreme_SSD/dhw/summaries/GBR_CoralTemp_full.rds",
                  crop=TRUE, mask=TRUE, downsample=FALSE)

GBR_CoralTemp_full <- readRDS( "/Volumes/Extreme_SSD/dhw/summaries/GBR_CoralTemp_full.rds")


process_CoralTemp(input = "/Volumes/Extreme_SSD/dhw/CRW/CRW_DHW",
                  polygon = GBR_hull, crs = "EPSG:7844",
                  combinedfilename = "/Volumes/Extreme_SSD/dhw/summaries/GBR_CoralTemp_DHW_full.rds",
                  crop=TRUE, mask=TRUE, downsample=FALSE, mc.cores = 10)

Alternatively, cdo offers a much faster workflow using cdo mergetime *.nc outfile and cropped to the GBR extent usingremapbil,target_grid_file.nc input.nc regridded_output.nc (note: CDO does not directly support cropping to irregular polygons, so remapbil sets to the extent (bbox).

Extract OISST

download_OISST() pulls global data from the ncei https as follows:


download_OISST(url = "https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/",
               start_date = "2025-01-01",
               end_date = "2026-01-31",
               dest_dir = "/Volumes/Extreme_SSD/dhw/OISST/",
               mc.cores = 10)

the full global dataset is 26.39GB in individual ~10-11MB .nc files.

process_OISST does three functions: 1) crop-mask to a polygon shape (here the GBR hull), 2) combine rast files using c(), 3) output to an .rds file for saving. The function will be paralellised (later) as it is slow Update (2/4/25) to run with parallel::mclapply() if mc.cores >1 - processing 15918 files across 10 cores now takes 88 seconds


process_OISST(input = "/Volumes/Extreme_SSD/dhw/OISST/",
              polygon = GBR_hull, crs = "EPSG:7844",
              combinedfilename = "/Volumes/Extreme_SSD/dhw/summaries/GBR_OISST_2026.rds",
              crop=TRUE, mask=TRUE, downsample=FALSE, mc.cores=10)

terra::rast("/Volumes/Extreme_SSD/dhw/summaries/GBR_OISST_full_2025.rds")

Extract ERA5

ERA5 doesn’t have a direct https link as far as I’m aware. The download_ERA5() function uses the ecmwfr interface to ‘ECMWF’ and ‘CDS’ Data Web Services to access the daily SST data via Copernicus. Once downloaded, the process_ERA5() function combines the outputs to a single rast file.


download_ERA5(
    start_year = 1940,
    end_year = 2025,
    ecmwfr_key = "381246a1-68c9-4e99-aab2-44c36d6da73d",
    timeout=60,
    dest_dir = "/Volumes/Extreme_SSD/dhw/ERA5/"
  )
  


process_ERA5(input = "/Volumes/Extreme_SSD/dhw/ERA5/", units = "celsius",
              polygon = GBR_hull, crs = "EPSG:7844",
              combinedfilename = "/Volumes/Extreme_SSD/dhw/summaries/GBR_ERA5_full.rds",
              crop=TRUE, mask=TRUE, downsample=FALSE)
terra::rast("/Volumes/Extreme_SSD/dhw/summaries/GBR_ERA5_full.rds")