ohi logo
OHI Science | Citation policy

Summary

This script takes the raw netCDF data and does the following:

  1. Calculates the annual mean for all years from 1958 - 2016
  2. Rescales each annual raster layer from 0 to 1 based on a biological threshold (Ω <= 1) and the proportional change compared to a historical mean
  3. Interpolates the data to gap-fill for cells where there is no data
  4. Resamples the rescaled raster layer to 1km^2 cell resolution
  5. Mask the resampled data to select only those cells within the ocean

Updates from 2016 assessment

A significant update to the process but not the results is the fact that we reran the OA pressure layers for all years 1958 - 2020 after receiving additional years of data from Ivan Lima on a regular grid. Previously the data was provided on an irregular grid requiring reprojection and additional interpolation. We don’t expect data changes to previous years but all years were rerun for consistency. Previously interpolation was done using Arcpy, but all necessary interpolation is done in R using Inverse Distance Weigthing.


Data Source

Reference: Feely et al.(2009)

Downloaded: July 19, 2017

Description: Aragonite Saturation State \(\Omega_{arg}\)

Native data resolution: 1 degree cells

Time range: 1880-1889 and 1958-2100, monthly data provided for each year. Future years are based on model projections for RCP 8.5. Previous years are hindcast/historical data.

Format: NetCDF

Notes about the data:

This data was shared with us by Ivan Lima from Woods Hole Institue for Oceanography in July 2017. The data came as NetCDFs with a resolution of about 1 degree. The data values are monthly average surface Ω aragonite saturation state.


Methods

Setup

#set options for all chunks in code
knitr::opts_chunk$set(warning=FALSE, message=FALSE,fig.width=6, fig.height=6)

#libraries
library(raster)
library(ncdf4)
library(maps)
library(parallel)
library(foreach)
library(doParallel)
library(RColorBrewer)

source("~/github/ohiprep/src/R/common.R")

#define paths for the raw data and OA folder held on git-annex on our NCEAS server, Mazu

raw_dir    = file.path(dir_M, 'git-annex/globalprep/_raw_data')
oagit_dir  = file.path(dir_M, 'git-annex/globalprep/prs_oa')

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

Load raw data

We received new data from Ivan Lima at WHOI from 1958 - 2016. This is provided in a single NetCDF file and on a regular 1x1 degree grid.

# read in 1958-2016 data

raw <- stack(file.path(raw_dir, 'WHOI/cesm_co2sys_1958-2016.1x1d.nc'), varname="OARG") #grab aragonite variable using varname
                    
plot(raw[[1]]) #we'll have to rotate the rasters

plot(rotate(raw[[1]]))
maps::map('world', col='gray95', fill=T, border='gray80', add=T)

Historical Mean

The historical mean Ω aragonite saturation state from 1880 - 1889 was calculated for OHI 2015. The same raster was used for OHI 2016 and is again used here for comparison.

hist <- raster(file.path(dir_M,'git-annex/globalprep/prs_oa/v2015/working/global_oa_1880_1889_arag_mean_moll.tif'))

plot(hist,main='Mean Ωaragonite 1880-1889', col=rev(cols), box=F,axes=F)

Create annual means

Since this data is now on a regular grid, but the historical data from 1880-1889 was transformed previously from an irregular grid, there will be issues matching up the resolution and extent so we will need to use this historical mean raster that we created for OHI 2015 and on, to use as a template for resampling.

# define mollweide projection
mollCRS <- CRS('+proj=moll')

#register parallel cores

registerDoParallel(10)

foreach (i = seq(1,708,by=12)) %dopar% {

  j = i+11
  
  yr = substr(names(raw[[i]]),2,5)
  
  yr_mean = raw[[i:j]]%>%
            calc(fun=function(x){mean(x,na.rm=T)})%>%
            rotate()%>%
            projectRaster(crs = mollCRS,over=T)%>%
            resample(hist,method = 'ngb')%>%
          writeRaster(.,filename = paste0(file.path(oagit_dir),'/v2017/int/annual_avg_moll/global_arag_avg_moll_',yr,'.tif'),overwrite=T)
 
}
#Here is the output for 2016

r_2016 <- raster(file.path(oagit_dir,'v2017/int/annual_avg_moll/global_arag_avg_moll_2016.tif'))

plot(r_2016,box=F,axes=F,main="Mean Ωaragonite saturation state 2016", col=rev(cols))

Rescale from 0 to 1

This pressure layer is rescaled so that all values lie between 0 and 1 using both a historical reference period and a biological reference point. All cells with values less than one, indicating an undersaturated state, are set equal to the highest stressor level, 1. For all other cells, rescaling the aragonite staturation state value to between 0 and 1 relies upon the change in saturation relative to the reference period.

Deviation from aragonite saturation state is determined for each year in the study period using this equation:

\[\Delta \Omega_{year} = \frac{(\Omega_{base} - \Omega_{year})}{(\Omega_{base} - 1)}\]

Note that the current value is subtracted from the baseline; this way, a reduction in \(\Omega\) becomes a positive pressure value. It is then normalized by the current mean state; so a decrease in \(\Omega\) while the current state is high indicates less pressure than the same decrease when the current state is near 1.

\(\Delta \Omega_{year}\) is then modified to account for increases in aragonite saturation state (pressure = 0) and arag sat state less than 1 (pressure = 1).

The oaRescale function rescales each of the annual rasters. If the current value is less than or equal to 1, it is set to 1, otherwise the value is calculated from the above equation.

#for each layer, all values <=1 are assigned a 1, otherwise old-new/(old-1)

oaRescale <- function(file){
  
  yr   = substr(file, nchar(file)-7, nchar(file)-4)  #get year of file
  mean = raster(file)                              #get seasonal mean aragonite raster for given year
  diff = (hist-mean)/(hist-1)
  mean[mean<=1] <- 1                                 #all values at or less than 1 are given a value of 1
  mean[mean>1] <- diff[mean>1]                     #all cells with values greater than 1 are swapped out with their amount of change scaled to how close to 1 
  mean[mean<0] <- 0                                  #all values less than 0 (indicating a decrease in acidity) are capped at 0

    writeRaster(mean, filename=paste0(oagit_dir, '/v2017/int/annual_avg_moll_rescaled/oa_rescaled_', yr, sep=""), format='GTiff', overwrite=TRUE)
}

files = list.files(file.path(oagit_dir, 'v2017/int/annual_avg_moll'), full.names=TRUE)

mclapply(files, oaRescale, mc.cores = 16)
r <- raster(file.path(oagit_dir,'v2017/int/annual_avg_moll_rescaled/oa_rescaled_2016.tif'))
plot(r,col=cols,box=F,axes=F, main = 'Rescaled Ωaragonite layer for 2016')

Interpolate

Since there are oceanic cells with no information in the raw data, we need to fill in these gaps. We do this by interpolating across the globe using the data we have with an Inverse Distance Weighting (IDW) function. Previously this was done using arcpy from ArcGIS.

library(gstat)
#register parallel cores

registerDoParallel(24)

files <- list.files(file.path(oagit_dir,'v2017/int/annual_avg_moll_rescaled'), full.names=TRUE)

foreach(file = files) %dopar%{

r  <- raster(file) #oa raster
yr <- substr(file, nchar(file)-7, nchar(file)-4)
xy <- data.frame(xyFromCell(r, 1:ncell(r)))                         #get xy coords into dataframe
v  <- getValues(r)                                                  #get cell values 
tmpdf <- cbind(xy, v)%>%filter(!is.na(v))                           #create dataframe of x,y, and values. remove NAs (throws error since these are cells we are interpolating over)
mg <- gstat(id = "v", formula = v~1, locations = ~x+y, data=tmpdf,
            nmax=7, set=list(idp = 2)) #define model. power function = 2, this is default for idw models
z <- interpolate(r, mg, progress='text')                            #interpolate across NA cells

writeRaster(z,filename  = paste0(file.path(oagit_dir),"/v2017/int/annual_avg_moll_rescaled_int/oa_resc_int_",yr,".tif"),overwrite=T)
}
plot(raster(file.path(oagit_dir,'v2017/int/annual_avg_moll_rescaled_int/oa_resc_int_2016.tif')), col=cols, box=FALSE, axes=FALSE, main='Rescaled and Interpolated Ωaragonite layer for 2016')

Resample & Land Mask

All pressure layers need to be resampled to 1km2 cell resolution. We have a template ocean raster with cells at this resolution that we use to resample all pressure layers. You won’t see any difference between the plot above and this one since we are using the nearest neighbor method when resampling which maintains the original cell value for each of resampled cell.

files <- list.files(file.path(oagit_dir,'v2017/int/annual_avg_moll_rescaled_int'),full.names=T)

foreach(file = files)%dopar% {

  yr <- substr(file,nchar(file)-7,nchar(file)-4)
  
  raster(file)%>%
  resample(ocean,method = 'ngb')%>%
  mask(ocean,filename = paste0(oagit_dir,'/v2017/output/oa_prs_layer_',yr,'.tif'),overwrite=T)
  
}

Final Pressure Layer

plot(raster(file.path(oagit_dir, 'v2017/output/oa_prs_layer_2016.tif')), box=F, axes=F, col=cols, main='Final Ocean Acidification \nPressure Layer 2016')

Gap-filled cells

We want to create a raster layer that shows all cells that were gap-filled. Since they were the same cells interpolated across all years, we only need to create one raster.

#Rescaled data before interpolation 
pre_int = raster(file.path(oagit_dir, 'v2017/int/annual_avg_moll_rescaled/oa_rescaled_2016.tif'))%>%
                resample(ocean, progress='text', method = 'ngb')

#after interpolation,
r_int = raster(file.path(oagit_dir, 'v2017/output/oa_prs_layer_2016.tif'))
    
#interpolated (or gap-filled) cells    
interp_cells = mask(r_int, pre_int, inverse=TRUE, filename = file.path(oagit_dir, 'v2017/output/oa_interpolated_cells.tif'))
plot(raster(file.path(oagit_dir, 'v2017/output/oa_interpolated_cells.tif')), col=cols, box=F, axes=F, main='Interpolated cells')

Citation information

Woods Hole Oceanographic Institution. 2014 update to data originally published in: Feely, R.A., S.C. Doney, and S.R. Cooley. 2009. Ocean acidification: Present conditions and future changes in a high-CO2 world. Oceanography 22(4):36–47