ohi logo
OHI Science | Citation policy

Summary

This script creates the Sea Surface Temperature (SST) layer for the 2016 global Ocean Health Index assessment.

Updates from previous assessment

The only difference in this layer is the reference point. For OHI 2015, each SST pressure layer was rescaled using the 99.99th quantile of the data from within each year. This has been changed to the 99.99th quantile across all years. It is a small difference, from 130 to 128 (number of anomalous weeks within a 5 year period).

Data Source

Data comes from CoRTAD version 5

See prs_sst/v2015/dataprep.R for preparation of the “annual_pos_anomalies” data.

Native Data Resolution: 4km
Description: Cortadv5_SSTA.nc = SST anomalies (weekly SST minus weekly climatological SST), weekly data for all years, degrees Kelvin Cortadv5_weeklySST.nc = SST, weekly data for all years, degrees Kelvin
Time Range: 1982 - 2012 (weekly averages across all years)
Format: NetCDF

Methods

  1. Extreme events per year based calculated as number of times SST anomoly exceeds SST Standard Deviation based on weekly values (annual_pos_anomalies data, see v2015/dataprep.R for analysis).
  2. Sum extreme events for five year periods to control for yearly variation.
  3. Change in extreme events: Subtract number of extreme events for each five year period from control period (1985-1989).
  4. Rescale “Change in extreme events” data to values between 0 and 1 by dividing by the 99.99th quantile among all years of data.

Setup

knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',message = FALSE, warning = FALSE)

#source('~/github/ohiprep/src/R/common.R')

source('../../../src/R/common.R')

library(raster)
library(RColorBrewer)
library(tidyverse)
library(rgdal)
library(doParallel)
library(foreach)

cols = rev(colorRampPalette(brewer.pal(11, 'Spectral'))(255)) # rainbow color scheme

#ocean raster at 1km
ocean = raster(file.path(dir_M, 'model/GL-NCEAS-Halpern2008/tmp/ocean.tif'))
ocean_shp <- readOGR(file.path(dir_M,'git-annex/globalprep/spatial/d2014/data'),layer = 'regions_gcs')
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/shares/ohi/git-annex/globalprep/spatial/d2014/data", layer: "regions_gcs"
## with 551 features
## It has 7 fields
## Integer64 fields read as strings:  rgn_id ant_id
land      <- ocean_shp%>%
              subset(rgn_typ %in% c('land','land-disputed','land-noeez'))

#set mollweide projection CRS
  mollCRS=crs('+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs')

# setwd if not running entire Rmd
  #setwd("globalprep/prs_sst/v2016")

Create 5 year cumulative sum of extreme events and calculate difference from historical

  l   <- list.files(file.path(dir_M,'git-annex/globalprep/prs_sst/v2015/tmp'),pattern='annual_pos_anomalies',full.names=TRUE)
  

# Get 5 year aggregates
ref_years <- c(grep(c('1985'), l), grep(c('1986'), l), grep(c('1987'), l), grep(c('1988'), l), grep(c('1989'), l))

  ref <- stack(l[ref_years])%>%sum(.) # This is the time period we are using for historical comparison (1985 - 1989)

  
  registerDoParallel(10)  
  
foreach(i = 1986:2008)%dopar%{ #i=2005
  
  print(i)
  
  yrs <- c(i:(i+4))
  
  s   <- stack(l[substr(l,81,84)%in%yrs])%>%sum(.)
  
  diff = overlay(s,ref,fun=function(x,y){x-y})%>% #calculate difference between recent 5 year cumulative sum and historical (1985-1989)
          mask(land,inverse=T)
  
  writeRaster(diff,filename = paste0('int/sst_diff_ocean_',yrs[1], '-', yrs[5],'.tif'),overwrite=T)
   
}

Reference Point

The layers are rescaled using a single reference point, the 99.99th quantile across all difference rasters.

diffs <- list.files('int',pattern = 'diff',full.names=T)

#get data across all years
vals <- c()

for(i in 1:length(diffs)){
  print(i)
  
  m <- diffs[i]%>%
    raster()%>%
    getValues()
  
  vals <- c(vals,m)
  
}
## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
#get min, max and 99.99th quantile

min_v <- min(vals,na.rm=T)
max_v <- max(vals,na.rm=T)
resc_num  <- quantile(vals,prob=0.9999,na.rm=T) ### 128

rescale <- read.csv("../../supplementary_information/v2016/reference_points_pressures.csv")
rescale$ref_point[rescale$pressure == "Sea Surface Temperature"] <- resc_num 
write.csv(rescale, "../../supplementary_information/v2016/reference_points_pressures.csv", row.names=FALSE)

The minimum value is -133, maximum value is 164 and the reference point is 128.

Rescaling

sprintf('Rescaling')

diffs <- list.files('int',pattern = 'diff',full.names=T)

resc_num <- read.csv("../../supplementary_information/v2016/reference_points_pressures.csv") %>%
  filter(pressure == "Sea Surface Temperature") %>%
  .$ref_point
resc_num <- as.numeric(as.character(resc_num))

foreach(i = 1:length(diffs)) %dopar%{

  r             <- raster(diffs[i])
  
  yrs           <- substr(diffs[i],20,28)
  
  projection(r) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
  
  out = projectRaster(r, crs=mollCRS,over=T)%>%
        calc(.,fun=function(x){ifelse(x>0,ifelse(x>resc_num,1,x/resc_num),0)})%>%
        resample(.,ocean,method='ngb',filename=paste0(dir_M,'/git-annex/globalprep/prs_sst/v2016/output/sst_',yrs,'_1985-1989.tif'),overwrite=T)

}

Results

res <- list.files(file.path(dir_M,'git-annex/globalprep/prs_sst/v2016/output'),full.names=T)

plot(raster(res[23]),col=cols,axes=F,main = 'Sea Surface Temperature Pressure Layer \n OHI 2016')

Citation information

Selig, E.R., K.S. Casey, and J.F. Bruno (2010), New insights into global patterns of ocean temperature anomalies: implications for coral reef health and management, Global Ecology and Biogeography, DOI: 10.1111/j.1466-8238.2009.00522.x.