Required libraries

library(raster)
library(gridExtra)
library(ggpubr)
library(rgdal)
library(RColorBrewer)
library(RStoolbox)

Load raster layers

setwd("N:/Prosjekter/311700_MAREANO/311720_Havbunnskart/c_Metodeprosjekter/xxx_Karbonlagring/OCAR/R")
dir(path = "../input")
##  [1] "land.dbf"                             
##  [2] "land.prj"                             
##  [3] "land.shp"                             
##  [4] "land.shx"                             
##  [5] "NorwegianTrough.cpg"                  
##  [6] "NorwegianTrough.dbf"                  
##  [7] "NorwegianTrough.prj"                  
##  [8] "NorwegianTrough.sbn"                  
##  [9] "NorwegianTrough.sbx"                  
## [10] "NorwegianTrough.shp"                  
## [11] "NorwegianTrough.shp.xml"              
## [12] "NorwegianTrough.shx"                  
## [13] "OCdensity_quantrf_mean.tif"           
## [14] "OCdensity_quantrf_tot.unc.percent.tif"
## [15] "SedRate3_quantrf_mean.tif"            
## [16] "SedRate3_quantrf_tot.unc.percent.tif"
oc <- raster("../input/OCdensity_quantrf_mean.tif")
oc.unc.pc <- raster("../input/OCdensity_quantrf_tot.unc.percent.tif")
sr <- raster("../input/SedRate3_quantrf_mean.tif")
sr.unc.pc <- raster("../input/SedRate3_quantrf_tot.unc.percent.tif")

Calculate organic carbon accumulation rates

OC density in kg m-3

Sedimentation rate in cm yr-1

OCAR in g m-2 yr-1

ocar <- oc * sr * 10
ocar
## class      : RasterLayer 
## dimensions : 2686, 2559, 6873474  (nrow, ncol, ncell)
## resolution : 500, 500  (x, y)
## extent     : 3230952, 4510452, 3023466, 4366466  (xmin, xmax, ymin, ymax)
## crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0.01553161, 66.18131  (min, max)
plot(ocar, main = "OC accumulation rate (g m-2 yr-1)")

Calculate uncertainty (percent)

Error propagation http://ipl.physics.harvard.edu/wp-uploads/2013/03/PS3_Error_Propagation_sp13.pdf

ocar.unc.pc <- sqrt(oc.unc.pc^2 + sr.unc.pc^2)
ocar.unc.pc
## class      : RasterLayer 
## dimensions : 2686, 2559, 6873474  (nrow, ncol, ncell)
## resolution : 500, 500  (x, y)
## extent     : 3230952, 4510452, 3023466, 4366466  (xmin, xmax, ymin, ymax)
## crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 50.04277, 1642.122  (min, max)
plot(ocar.unc.pc, main ="Total uncertainty (% of mean)")

Calculate uncertainty (absolute)

ocar.unc <- ocar * ocar.unc.pc / 100
ocar.unc
## class      : RasterLayer 
## dimensions : 2686, 2559, 6873474  (nrow, ncol, ncell)
## resolution : 500, 500  (x, y)
## extent     : 3230952, 4510452, 3023466, 4366466  (xmin, xmax, ymin, ymax)
## crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 0.1975928, 57.90045  (min, max)
plot(ocar.unc, main ="Total uncertainty (g m-2 yr-1)")

Calculation of the OC accumulation and the associated uncertainty (in Tg/yr)

A <- res(ocar)[1]*res(ocar)[2] #Area of one pixel

oca <- cellStats(ocar, sum)*A/1000000000000
oca_unc <- cellStats(ocar.unc, sum)*A/1000000000000

print(paste0("OC accumulation = ", round(oca, 2), " Tg/yr ± ", round(oca_unc, 2), " Tg/yr"))
## [1] "OC accumulation = 1.43 Tg/yr ± 2.07 Tg/yr"

OC accumulation in the Norwegian Trough

list.files(path = "../input",pattern = ".shp")
## [1] "land.shp"                "NorwegianTrough.shp"    
## [3] "NorwegianTrough.shp.xml"
nt <- readOGR(dsn = "../input", layer = "NorwegianTrough")
## OGR data source with driver: ESRI Shapefile 
## Source: "N:\Prosjekter\311700_MAREANO\311720_Havbunnskart\c_Metodeprosjekter\xxx_Karbonlagring\OCAR\input", layer: "NorwegianTrough"
## with 1 features
## It has 1 fields
ocar_nt <- crop(ocar, nt)
ocar.unc_nt <- crop(ocar.unc, nt)

oca_nt <- cellStats(ocar_nt, sum)*A/1000000000000
oca.unc_nt <- cellStats(ocar.unc_nt, sum)*A/1000000000000

print(paste0("OC accumulation in the Norwegian Trough: ", round(oca_nt, 2), " Tg/yr ± ", round(oca.unc_nt, 2), " Tg/yr"))
## [1] "OC accumulation in the Norwegian Trough: 1.24 Tg/yr ± 1.3 Tg/yr"

Export as geoTiff

writeRaster(ocar, "../output/OCAR.tif", overwrite = TRUE)
writeRaster(ocar.unc, "../output/OCAR_tot.unc.tif", overwrite = TRUE)
writeRaster(ocar.unc.pc, "../output/OCAR_tot.unc.percent.tif", overwrite = TRUE)

Figure for publication

col.pal1 <- brewer.pal(8,"YlGnBu")
col.pal2 <- brewer.pal(8,"YlOrRd")

AoI <- extent(3450000,4450000,3100000,4300000)
land <- readOGR(dsn = "../input", layer = "land")
## OGR data source with driver: ESRI Shapefile 
## Source: "N:\Prosjekter\311700_MAREANO\311720_Havbunnskart\c_Metodeprosjekter\xxx_Karbonlagring\OCAR\input", layer: "land"
## with 8186 features
## It has 1 fields
land <- crop(land, AoI)
ocar <- crop(ocar, AoI)
ocar.unc <- crop(ocar.unc, AoI)
ocar.unc.pc <- crop(ocar.unc.pc, AoI)

p1 <- ggR(ocar, geom_raster = TRUE) + 
  geom_polygon(data = land, aes(x = long, y = lat, group = group)) + 
  scale_fill_gradientn(colours = col.pal1, na.value = NA, name = "") +
  ggtitle(expression(OC~Accumulation~Rate~(g~m^{-2}~yr^{-1}))) +
  theme(axis.title = element_blank()) +
  ggpubr::rotate_y_text()
## Regions defined for each Polygons
p2 <- ggR(ocar.unc, geom_raster = TRUE) + 
  geom_polygon(data = land, aes(x = long, y = lat, group = group)) + 
  scale_fill_gradientn(colours = col.pal2, na.value = NA, name = "") +
  ggtitle(expression(Total~Uncertainty~(g~m^{-2}~yr^{-1}))) +
  theme(axis.title = element_blank()) +
  ggpubr::rotate_y_text()
## Regions defined for each Polygons
p3 <- ggR(ocar.unc.pc, geom_raster = TRUE) + 
  geom_polygon(data = land, aes(x = long, y = lat, group = group)) + 
  scale_fill_gradientn(colours = col.pal2, na.value = NA, name = "") +
  ggtitle(expression(Total~Uncertainty~('%'~of~Mean))) +
  theme(axis.title = element_blank()) +
  ggpubr::rotate_y_text()
## Regions defined for each Polygons
tiff("../figs/ocar_figure_v1.tif", width = 30, height = 10, units = "cm", res = 500, compression = "lzw")
grid.arrange(p1, p2, p3, ncol = 3)
dev.off()
## png 
##   2
tiff("../figs/ocar_figure_v2.tif", width = 20, height = 10, units = "cm", res = 500, compression = "lzw")
grid.arrange(p1, p2, ncol = 2)
dev.off()
## png 
##   2