library(raster)
library(gridExtra)
library(ggpubr)
library(rgdal)
library(RColorBrewer)
library(RStoolbox)
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")
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)")
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)")
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)")
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"
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"
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)
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