OHI Science | Citation policy
This script takes the raw netCDF data and does the following:
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.
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.
#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("../../../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 schemeWe 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)
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)
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))
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')
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')
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.
ocean <- raster(file.path(dir_M, "git-annex/globalprep/spatial/v2017/ocean.tif"))
files <- list.files(file.path(oagit_dir,'v2017/int/annual_avg_moll_rescaled_int'),
full.names = TRUE)
foreach(file = files)%dopar% {
#for(file in files){ # file=files[1] # for some reason cores weren't all working for me. Did it the long way!
yr <- substr(file,nchar(file)-7,nchar(file)-4)
cat(yr)
raster(file)%>%
resample(ocean, method = 'ngb')%>%
mask(ocean, filename = paste0(oagit_dir, '/v2017/output/oa_prs_layer_',yr,'.tif'), overwrite=TRUE)
}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')
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')
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