Data processing in R

1. Initialise data

load libraries, import raster datasets, download spatial data for GBR reefs (see dhw::download_spatial_gbr for details)

library(tidyverse)
library(sf)
library(RColorBrewer)
library(dhw)
library(foreach)

# --- Load data ---
GBR_OISST <- rast("/Users/rof011/GBR-dhw/datasets/GBR_2025_OISST_DHW.rds")
GBR_CRW_DHW <- rast("/Users/rof011/GBR-dhw/datasets/GBR_CoralTemp_DHW_full.rds")
gbr_reefs_base<- download_gbr_spatial(return = "base") |> st_transform(4326)
gbr_reefs <- download_gbr_spatial(return = "combined") |> st_transform(4326)
GBRMPpolygon <- vect(gbr_reefs |> st_transform(4326))

2. Calculate DHW for OISST

Compute the climatology baseline and extract DHW (Degree Heating Weeks) from the OISST dataset.

OISST_climatology <- create_climatology(GBR_OISST, quiet=TRUE)
GBR_OISST_DHW <- OISST_climatology$dhw

3. Calculate annual DHWmax

Compute yearly maximum DHW for both OISST and CoralTemp CRW datasets cropped to the GBRMP boundary

OISST_DHW <- summarise_raster(GBR_OISST_DHW, index = "years", fun = max, na.rm = TRUE, overwrite = TRUE) |>
  crop(GBRMPpolygon)

CRW_DHW <- summarise_raster(GBR_CRW_DHW, index = "years", fun = max, na.rm = TRUE, overwrite = TRUE) |>
  crop(GBRMPpolygon)

4. Extract DHW metrics per reef

Extract per-reef DHW values by calculating area-weighted mean DHW for each reef polygon and each year.

OISST_DHW_reefs <- foreach(i = seq_len(nlyr(OISST_DHW)), .combine = rbind) %do% {
  gbr_reefs |>
    mutate(
      date = rep(time(OISST_heat[[i]]), nrow(gbr_reefs)),
      sst = exactextractr::exact_extract(OISST_DHW[[i]], gbr_reefs, progress = FALSE, fun = "weighted_mean", weights = "area")
    )
}


CRW_DHW_reefs <- foreach(i = seq_len(nlyr(CRW_DHW)), .combine = rbind) %do% {
  gbr_reefs |>
    mutate(
      date = rep(time(OISST_heat[[i]]), nrow(gbr_reefs)),
      sst = exactextractr::exact_extract(CRW_DHW[[i]], gbr_reefs, progress = FALSE, fun = "weighted_mean", weights = "area")
    )
}

5. Format data for d3 heatmap

Extract centroids and latitude of each reef for d3 heatmap.

OISST_DHW_reefs_centroids <- OISST_DHW_reefs %>%
  st_centroid() %>%
  mutate(
    lon = st_coordinates(.)[, 1],
    lat = st_coordinates(.)[, 2]
  ) %>%
  st_drop_geometry() %>%
  as.data.frame() |>
  rename(dhw = sst) |>
  rename(year = date) |>
  mutate(dhw = round(dhw,1)) |>
  mutate(lat = round(lat,1)) |>
  select(-area, -lon) |>
  pivot_wider(id_cols=c("LABEL_ID", "GBR_NAME", "lat"), values_from="dhw", names_from="year") |>
  rename(ID = LABEL_ID, Reef=GBR_NAME)

CRW_DHW_reefs_centroids <- CRW_DHW_reefs %>%
  st_centroid() %>%
  mutate(
    lon = st_coordinates(.)[, 1],
    lat = st_coordinates(.)[, 2]
  ) %>%
  st_drop_geometry() %>%
  as.data.frame() |>
  rename(dhw = sst) |>
  rename(year = date) |>
  mutate(dhw = round(dhw,1)) |>
  mutate(lat = round(lat,1)) |>
  select(-area, -lon) |>
  pivot_wider(id_cols=c("LABEL_ID", "GBR_NAME", "lat"), values_from="dhw", names_from="year") |>
  rename(ID = LABEL_ID, Reef=GBR_NAME)


#write.csv(OISST_DHW_reefs_centroids, "/dhw3/data/OISST_DHWmax.csv")
#write.csv(CRW_DHW_reefs_centroids, "/dhw3/data/CRW_DHWmax.csv")