This tutorial was created to help users download, view, summarise and map coral cover predictions generated by the Virtual Reef Diver program. The software used in the tutorial is R version 4.0.3.

Download example datasets

Virtual Reef Diver model predictions:

Model predictions can be downloaded from the Virtual Reef Diver website at https://vrd-predictions.s3-ap-southeast-2.amazonaws.com/index.html. The file is named pred.mean_2016.zip.

Management Areas of the Great Barrier Reef:

Boundaries of the Great Barrier Reef management areas, in a shapefile format, can be downloaded at: https://geoportal.gbrmpa.gov.au/datasets/c5640ed04c594762a9639802dbd56a22_0/data?geometry=117.416%2C-24.882%2C179.116%2C-10.267.

Unzip the files and save them to your working directory.

Clipping the rasters with a shapefile

In this example, 2016 predictive values in coral cover are extracted for the management region of Cairns/Cooktown. # R commands

Set working directory:

# The working directory is the folder contained the saved files. 
#setwd('your/directory')

Loading R packages:

library(dplyr)
library(tidyr)
library(raster)
library(leaflet)
library(sf)
library(ggplot2)

# Suppress warnings
options(warn=-1)

Reading the raster file and shapefile:

# Read raster 
coral_pred <- raster("pred.mean_2016.tif")    

# Read shapefile and transform the coordinates using the projection from raster
mngt_spat <- st_read("Management_Areas_of_the_Great_Barrier_Reef_Marine_Park.shp",quiet = TRUE) %>% st_transform(crs(coral_pred))
unique(mngt_spat$AREA_DESCR)
## [1] Mackay/Capricorn Management Area      Townsville/Whitsunday Management Area
## [3] Cairns/Cooktown Management Area       Far Northern Management Area         
## 4 Levels: Cairns/Cooktown Management Area ... Townsville/Whitsunday Management Area
# Keep the Cairns Cooktown Management area only
mngt_sub <- subset(mngt_spat, AREA_DESCR=="Cairns/Cooktown Management Area")

Cropping the raster using the extent from the management area:

coral_pred_region <- crop(coral_pred, extent(mngt_sub))

Plotting the predictions on a map:

pal <- colorNumeric("Spectral", domain = values(coral_pred_region),
  na.color = "transparent")

p.map <- leaflet(option=leafletOptions(zoomControl = FALSE)) %>% 
  addTiles() %>%
  addRasterImage(coral_pred_region, colors = pal, opacity = 0.9) %>%
  addLegend(pal = pal, values = values(coral_pred_region),
    title = "Coral cover")

p.map

Visualising the regional distribution of coral cover and estimated mean (red line) for the year 2016:

coral_pred_df <- as.data.frame(coral_pred_region) 

ggplot(coral_pred_df, aes(x= pred.mean_2016)) + theme_bw() +
    geom_histogram(binwidth=.05, colour="black", fill="white")+
    geom_vline(aes(xintercept=mean(pred.mean_2016, na.rm=T)),   
               color="red", linetype="dashed", size=1) + xlab("Predicted regional coral cover")