Skip to content
Snippets Groups Projects
Commit 0ea57ed4 authored by Fabian Stenzel's avatar Fabian Stenzel
Browse files

updated package renamed a few functions and added get_biome_names

parent c94df3d4
No related branches found
No related tags found
1 merge request!4full review
Showing
with 1528 additions and 722 deletions
^.*\.Rproj$
^\.Rproj\.user$
.Rproj.user
.Rhistory
.RData
.Ruserdata
Package: biospheremetrics
Type: Package
Title: Biosphere integrity metrics for LPJmL
Version: 0.1.1
Version: 0.1.2
Author: Fabian Stenzel, Johanna Braun, Jannes Breier
Maintainer: Fabian Stenzel <stenzel@pik-potsdam.de>
Description: Functions to compute Biosphere integrity metrics M-COL and M-ECO based on output from LPJmL
License: GPL-3
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.2.2
RoxygenNote: 7.2.3
Depends:
lpjmliotools
# Generated by roxygen2: do not edit by hand
export(average_nyear_window)
export(calcDeltaV)
export(calcGamma)
export(calcMCOL)
export(calcMECO)
export(classify_biomes)
export(disaggregateGammaIntoBiomes)
export(evaluateGamma)
export(gammaCrossTable)
export(disaggregateMECOintoBiomes)
export(get_biome_names)
export(list_needed_outputs)
export(mecoCrossTable)
export(plotBiomeAverages)
export(plotBiomeAveragesToScreen)
export(plotBiomes)
export(plotBiomesAverage)
export(plotBiomesAveragesToScreen)
export(plotBiomesToScreen)
export(plotGammaCrossTable)
export(plotGammaCrossTableToScreen)
export(plotGammaMap)
export(plotGammaMapToScreen)
export(plotGammaRadial)
export(plotGammaRadial4)
export(plotGammaRadialToScreen)
export(plotMCOL)
export(plotMCOLmap)
export(plotMCOLovertime)
export(plotMECOcrossTable)
export(plotMECOcrossTableToScreen)
export(plotMECOmap)
export(plotMECOmapToScreen)
export(plotMECOradial)
export(plotMECOradialPanel)
export(plotMECOradialToScreen)
export(plot_biomes)
export(readGammaData)
export(readMCOLdata)
export(readMECOData)
export(replaceRefDataWithAverageRefBiomeCell)
importFrom(magrittr,"%>%")
......@@ -78,9 +78,9 @@ plotMCOLmap <- function(data, file, title, legendtitle, zeroThreshold = 0.1, eps
#' }
#' @export
plotMCOLovertime <- function(mcolData, file, firstyr, plotyrs, highlightyrs = 2000, minVal = 0,
maxVal = 100, legendpos = "topleft", ext = FALSE, eps = FALSE){
maxVal = 100, legendpos = "topleft", ext = FALSE, eps = FALSE, ref = "pi"){
lastyr = firstyr + length(mcolData$npp_act_overtime) - 1
colz = c("slateblue","gold","green3","red3","darkorange","black")
colz = c("slateblue","gold","green3","darkorange","black","red3")
if (eps) {
file = strsplit(file,".",fixed=TRUE)[[1]]
file = paste(c(file[1:(length(file) - 1)],"eps"),collapse=".")
......@@ -90,29 +90,33 @@ plotMCOLovertime <- function(mcolData, file, firstyr, plotyrs, highlightyrs = 20
png(file, width=3.5, height = 3, units = "in", res = 300, pointsize = 6,type="cairo")
}
par(bty="o",oma=c(0,0,0,0),mar=c(4,5,1,3))
plot(x=seq(firstyr,lastyr,1),y=mcolData$npp_pot_overtime,ylab="GtC/yr",xlab="Year",xlim=plotyrs,
ylim=c(minVal, maxVal),type = "l",col=colz[1],xaxs="i",yaxs="i")
plot(NA,ylab="GtC/yr",xlab="Year",xlim=plotyrs,
ylim=c(minVal, maxVal),xaxs="i",yaxs="i")
grid()
lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_pot_overtime,type = "l",col=colz[1])
lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_act_overtime,type = "l",col=colz[2])
lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_eco_overtime,type = "l",col=colz[3])
lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_luc_overtime,type = "l",col=colz[4])
lines(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime,type = "l",col=colz[5])
par(bty="n",oma=c(0,0,0,0),mar=c(4,5,1,3), new = T)
if (ref == "pi") {
plot(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime_perc_piref,ylab="",xlab="",xlim=plotyrs,
ylim=c(10, 30),type = "l",col=colz[4],xaxs="i", yaxs="i", axes = F)
ylim=c(10, 30),type = "l",col=colz[6],xaxs="i", yaxs="i", axes = F)
} else if (ref == "act") {
plot(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime,ylab="",xlab="",xlim=plotyrs,
ylim=c(10, 30),type = "l",col=colz[4],xaxs="i", yaxs="i", axes = F)
ylim=c(10, 30),type = "l",col=colz[6],xaxs="i", yaxs="i", axes = F)
}else stop(paste0("Unknown value for parameter ref: ",ref," - Aborting."))
axis(side = 4, col = colz[4],col.axis = colz[4])
mtext(text = "%", col=colz[4], side = 4,line = 2)
axis(side = 4, col = colz[6],col.axis = colz[6])
mtext(text = "%", col=colz[6], side = 4,line = 2)
if (!is.null(highlightyrs)){
for (y in highlightyrs){
lines(x=c(y,y),y=c(minVal,maxVal),col="grey40")
}
}
legend(legendpos,legend = c("NPPpot (PNV)","NPPact (landuse)","NPPeco","M-COL [% NPPpi]"),col=colz[1:4] ,lty=1,cex=1)
legend(legendpos,legend = c("NPPpot (PNV)","NPPact (landuse)","NPPeco","NPPluc","M-COLabs","M-COL [% NPPpi]"),col=colz ,lty=1,cex=1)
dev.off()
}
#' Calculate the ecosystem change metric gamma between 2 simulations/timesteps
......@@ -510,13 +514,15 @@ calcMCOL <- function(inFol_lu, inFol_pnv, startyr, stopyr, gridbased = T, npp_th
timber_harvest_overtime <- colSums(timber*cellarea)/10^15 # from gC/m2 to GtC
fire_overtime <- colSums(fire*cellarea)/10^15 # from gC/m2 to GtC
mcol_overtime <- harvest_cft_overtime + rharvest_cft_overtime +
harvest_grasslands_overtime + harvest_bioenergy_overtime +
timber_harvest_overtime + fire_overtime + npp_luc_overtime
mcol_overtime_perc_piref <- mcol_overtime/mean(npp_pot_overtime[1:10])*100
mcol <- harvest_cft + rharvest_cft + harvest_grasslands + harvest_bioenergy + timber + fire + ynpp_potential - ynpp
mcol[ynpp_potential<npp_threshold] <- 0 # set to 0 below lower threshold of NPP
mcol_luc <- ynpp_potential - ynpp
mcol_luc_piref <- rep(rowMeans(ynpp_potential[,1:10]),times = length(ynpp[1,])) - ynpp # always compare to pi_ref
mcol_harvest <- harvest_cft + rharvest_cft + harvest_grasslands + harvest_bioenergy + timber + fire
mcol <- mcol_harvest + mcol_luc
mcol[abs(ynpp_potential)<npp_threshold] <- 0 # set to 0 below lower threshold of NPP
mcol_perc <- mcol/ynpp_potential*100 #actual NPPpot as ref
mcol_perc_piref <- mcol/rowMeans(ynpp_potential[,1:10])*100 # NPPpi as ref
......@@ -529,7 +535,8 @@ calcMCOL <- function(inFol_lu, inFol_pnv, startyr, stopyr, gridbased = T, npp_th
harvest_cft_overtime = harvest_cft_overtime, npp_luc_overtime = npp_luc_overtime,
rharvest_cft_overtime = rharvest_cft_overtime, fire_overtime = fire_overtime,
timber_harvest_overtime = timber_harvest_overtime, harvest_cft = harvest_cft,
grassland_scaling_factor_cellwise = grassland_scaling_factor_cellwise))
grassland_scaling_factor_cellwise = grassland_scaling_factor_cellwise,
mcol_harvest = mcol_harvest, mcol_luc = mcol_luc, mcol_luc_piref = mcol_luc_piref))
}else{
return(list(mcol_overtime = mcol_overtime, mcol = mcol, mcol_perc = mcol_perc,
mcol_overtime_perc_piref = mcol_overtime_perc_piref,
......@@ -538,9 +545,9 @@ calcMCOL <- function(inFol_lu, inFol_pnv, startyr, stopyr, gridbased = T, npp_th
npp_eco_overtime = npp_eco_overtime, #npp_bioenergy_overtime = npp_bioenergy_overtime,
harvest_cft_overtime = harvest_cft_overtime, npp_luc_overtime = npp_luc_overtime,
rharvest_cft_overtime = rharvest_cft_overtime, fire_overtime = fire_overtime,
timber_harvest_overtime = timber_harvest_overtime, harvest_cft = harvest_cft))
timber_harvest_overtime = timber_harvest_overtime, harvest_cft = harvest_cft,
mcol_harvest = mcol_harvest, mcol_luc = mcol_luc, mcol_luc_piref = mcol_luc_piref))
}
}
#' Plot absolute MCOL, overtime, maps, and npp into given folder
......@@ -575,19 +582,34 @@ plotMCOL <- function(mcolData, outFol, plotyears, minVal, maxVal, legendpos,
print(paste0("Plotting MCOL figures"))
dir.create(file.path(outFol),showWarnings = F)
lpjmliotools::plotGlobal(data = rowMeans(mcolData$mcol[,(mapindex-mapyear_buffer):(mapindex+mapyear_buffer)]),
file = paste0(outFol,"MCOL",mapyear,"_absolute.png"),
title = paste0("MCOL in ",mapyear), min = 0, max = 1000,
legendtitle = "GtC", legYes = T, onlyPos = F, eps = eps, type = "lin")
file = paste0(outFol,"MCOL",mapyear,"_absolute.png"), type = "exp",
title = paste0("MCOL_abs in ",mapyear), pow2min = 0, pow2max = 12,
legendtitle = "GtC", legYes = T, onlyPos = F, eps = eps)
lpjmliotools::plotGlobal(data = rowMeans(mcolData$mcol_luc[,(mapindex-mapyear_buffer):(mapindex+mapyear_buffer)]),
file = paste0(outFol,"MCOL",mapyear,"_luc.png"), type = "exp",
title = paste0("MCOL_luc in ",mapyear), pow2min = 0, pow2max = 12,
legendtitle = "GtC", legYes = T, onlyPos = F, eps = eps)
lpjmliotools::plotGlobal(data = rowMeans(mcolData$mcol_luc_piref[,(mapindex-mapyear_buffer):(mapindex+mapyear_buffer)]),
file = paste0(outFol,"MCOL",mapyear,"_luc_piref.png"), type = "exp",
title = paste0("MCOL_luc piref in ",mapyear), pow2min = 0, pow2max = 12,
legendtitle = "GtC", legYes = T, onlyPos = F, eps = eps)
lpjmliotools::plotGlobal(data = rowMeans(mcolData$mcol_harvest[,(mapindex-mapyear_buffer):(mapindex+mapyear_buffer)]),
file = paste0(outFol,"MCOL",mapyear,"_harv.png"), type = "exp",
title = paste0("MCOL_harv in ",mapyear), pow2min = 0, pow2max = 12,
legendtitle = "GtC", legYes = T, onlyPos = F, eps = eps)
plotMCOLovertime(mcolData = mcolData, file = paste0(outFol,"MCOL_overtime_LPJmL_",plotyears[1],"-",plotyears[2],".png"),
firstyr = startyr, plotyrs = plotyears, minVal = minVal, ref = "pi",
legendpos = legendpos, maxVal = maxVal, eps = eps, highlightyrs = highlightyear)
plotMCOLmap(data = rowMeans(mcolData$mcol_perc[,(mapindex-mapyear_buffer):(mapindex+mapyear_buffer)]),
file = paste0(outFol,"MCOL",mapyear,"_LPJmL.png"),legendtitle = "% of NPPpot", eps = eps,
title = paste0(mapindex-mapyear_buffer, " - ",mapindex+mapyear_buffer) )
title = paste0("MCOL_perc ",mapyear-mapyear_buffer, " - ",mapyear+mapyear_buffer) )
plotMCOLmap(data = rowMeans(mcolData$mcol_perc_piref[,(mapindex-mapyear_buffer):(mapindex+mapyear_buffer)]),
file = paste0(outFol,"MCOL_piref_",mapyear,"_LPJmL.png"),
title = paste0(mapindex-mapyear_buffer, " - ",mapindex+mapyear_buffer),legendtitle = "% of NPPpi", eps = eps)
lpjmliotools::plotGlobalMan(data = rowMeans(mcolData$ynpp[,(mapindex-mapyear_buffer):(mapindex+mapyear_buffer)]),file = paste0(outFol,"NPP",mapyear,"_LPJmL.png"), brks = seq(0,1000,100),
title = paste0("MCOL_perc ",mapyear-mapyear_buffer, " - ",mapyear+mapyear_buffer),legendtitle = "% of NPPpi", eps = eps)
lpjmliotools::plotGlobalMan(data = rowMeans(mcolData$ynpp[,(mapindex-mapyear_buffer):(mapindex+mapyear_buffer)]),
file = paste0(outFol,"NPP",mapyear,"_LPJmL.png"), brks = seq(0,1000,100),
palette = c("orangered4","orangered","orange","gold","greenyellow","limegreen","green4","darkcyan","darkslategrey","navy"),
title = paste0("NPP average ",mapindex-mapyear_buffer, "-",mapindex+mapyear_buffer),legendtitle = "gC/m2",legYes = T)
title = paste0("NPP average ",mapyear-mapyear_buffer, "-",mapyear+mapyear_buffer),
legendtitle = "gC/m2",legYes = T)
} # end of plotMCOL
This diff is collapsed.
#' Calculate averages (mean) for defined window sizes
#'
#' Define window sizes (nyear_window) to be used to calculate averages (mean)
#' for each window (`dim(x)[3] / nyear_window`). Instead of discrete windows,
#' also moving averages can be computed as well as years inbetween interpolated.
#'
#' @param x LPJmL output array with `dim(x)=c(cell, month, year)`
#'
#' @param nyear_window integer, if supplied it defines the years for each window
#' to be averaged over in `dim(x)[3]`. If `nyear_window == 1` values are used
#' directly (instead of calculating an average). nyear_window has to be smaller
#' than `dim(x)[3]` and `dim(x)[3]` is ideally a multipe of nyear_window.
#' Defaults to `NULL`
#'
#' @param moving_average logical. If `TRUE` moving average is computed. start
#' and end are interpolated using spline interpolation.
#'
#' @param interpolate logical. If `TRUE` and nyear_window is defined (with
#' `moving_average == FALSE` years are interpolated (spline) to return array
#' with same dimensions as `x` (mainly`dim(x)[3]` -> year).
#'
#' @param nyear_reference integer, if supplied (default NULL), it defines a
#' time_span for ideally reference runs to be used as a baseline. E.g.
#' `nyear_reference = 30` to be used for preindustrial climate reference.
#'
#' @return array with same amount of cells and months as x. 3rd dimension is
#' defined by nyear_window, basically `dim(x)[3]/nyear_window` or equal to
#' dim(x)[3] if `moving_average == TRUE` or `interpolate == TRUE`
#'
#' @md
#' @importFrom magrittr %>%
#' @export
average_nyear_window <- function(x,
nyear_window = NULL,
moving_average = FALSE,
interpolate = FALSE,
nyear_reference = NULL) {
third_dim <- names(dim(x))[
!names(dim(x)) %in% c("cell", "year")
] %>% {
if (rlang::is_empty(.)) NULL else .
}
if (length(third_dim) > 1) {
stop(paste0("x has to have dimensions \"cell\", \"year\" and can have ",
"one third dimension (e.g. \"month\", \"year\""))
}
# check validity of x dimensions
if (!all(names(dim(x)) %in% c("cell", third_dim, "year"))) {
stop(paste0("x has to have dimensions \"cell\"",
ifelse(is.null(third_dim),
"",
paste0(", \"", third_dim, "\"")),
" and \"year\"."))
}
# moving average function - spline interpolation to fill NAs at start/end
moving_average_fun <- function(x, n) {
stats::filter(x, rep(1 / n, n), sides = 2) %>%
zoo::na.spline()
}
# utility function to interpolate inbetween averaging windows via spline
interpolate_spline <- function(x, y, nyear_window) {
rep(NA, dim(y)["year"]) %>%
`[<-`(seq(round(nyear_window / 2), dim(y)["year"], nyear_window),
value = x) %>%
zoo::na.spline()
}
# if nyear_window is supplied not all years are used for averaging
if (!is.null(nyear_window)) {
if (!is.null(nyear_reference)) {
orig_x <- x
x <- lpjmlKit::subset_array(x, list(year = 1:nyear_reference))
}
# only valid for nyear_window < years of x (dim(x)[3])
if (nyear_window > 1 & nyear_window <= dim(x)["year"]) {
# check if multiple (can also be left out)
# if (dim(x)[3] %% nyear_window == 0) {
if (moving_average) {
y <- aperm(apply(x,
c("cell", third_dim),
moving_average_fun,
nyear_window),
c("cell", third_dim, ""))
} else {
# calculate mean for discret windows/bins with size of nyear_window
y <- array(x,
# set correct dimensions (with names)
dim = c(dim(x)[c("cell", third_dim)],
year = nyear_window,
window = dim(x)[["year"]] / nyear_window),
# set correct dimensions names with nyear and windows
dimnames = append(dimnames(x)[c("cell", third_dim)],
list(year = seq_len(nyear_window),
window = dimnames(x)[["year"]][
seq(round(nyear_window / 2),
dim(x)[["year"]],
nyear_window)
]))) %>%
apply(c("cell", third_dim, "window"), mean)
if (interpolate) {
y <- aperm(apply(y,
c("cell", third_dim),
interpolate_spline,
x,
nyear_window),
c("cell", third_dim, ""))
}
}
# }
} else if (nyear_window == 1) {
y <- x
} else {
stop(paste0("Amount of nyear_window (", nyear_window, ") not supported."))
}
# recycle nyear_reference subset for original x (years)
if (!is.null(nyear_reference)) {
# multiple factor
nmultiple <- round(dim(orig_x)[["year"]] / nyear_reference)
replace_multiple_id <- nmultiple * dim(y)[[3]]
# if average window is returned as year dimension
if (!moving_average & !interpolate) {
z <- array(NA,
dim = c(dim(y)[c("cell", third_dim)],
window = replace_multiple_id),
dimnames = append(dimnames(y)[c("cell", third_dim)],
list(window = rep(dimnames(y)[[3]],
nmultiple))))
# return as original year dimension
} else {
# years vector also for non multiples (subset only partly recycled)
years <- rep(NA, dim(orig_x)[["year"]]) %>%
`[<-`(, value = dimnames(x)[["year"]]) %>%
suppressWarnings()
z <- array(NA,
dim = dim(orig_x),
dimnames = append(dimnames(y)[c("cell", third_dim)],
list(year = years)))
}
# recycle subset y for rest of original (x) years in z
z[, , seq_len(replace_multiple_id)] <- y
# check if not multiple - then only partly recylce array
if ((dim(z)[3] - replace_multiple_id) > 0) {
z[, , (replace_multiple_id + 1):dim(z)[3]] <- (
y[, , seq_len(dim(z)[3] - replace_multiple_id)]
)
}
return(z)
} else {
# rename dimnames of array
if (all(dim(y) == dim(x))) {
dim(y) <- dim(x)
dimnames(y) <- dimnames(x)
}
}
} else {
y <- apply(x, c("cell", third_dim), mean)
if (is.null(dim(y))) {
y <- array(
y,
dim = c(cell = dim(x)[["cell"]], 1),
dimnames = list(cell = dimnames(x)[["cell"]], 1)
)
}
}
return(y)
}
This diff is collapsed.
......@@ -30,12 +30,12 @@ plot_biomes <- function(biome_data,
lpjml_extent <- c(-180, 180, -60, 85)
bounding_box <- system.file("extdata", "ne_110m_wgs84_bounding_box.shp",
package = "pbfunctions") %>%
package = "biospheremetrics") %>%
rgdal::readOGR(layer = "ne_110m_wgs84_bounding_box", verbose = FALSE) %>%
{ if(to_robinson) sp::spTransform(., sp::CRS("+proj=robin")) else . } # nolint
countries <- system.file("extdata", "ne_110m_admin_0_countries.shp",
package = "pbfunctions") %>%
package = "biospheremetrics") %>%
rgdal::readOGR(layer = "ne_110m_admin_0_countries", verbose = FALSE) %>%
crop(., lpjml_extent) %>%
{ if(to_robinson) sp::spTransform(., CRS("+proj=robin")) else . } # nolint
......@@ -47,7 +47,7 @@ plot_biomes <- function(biome_data,
"#FFFFD4", "white", "#dad4d4")
biome_mapping <- system.file("extdata", "biomes.csv",
package = "pbfunctions") %>%
package = "biospheremetrics") %>%
readr::read_csv2()
names(biome_cols) <- biome_mapping$short_name
......
id;name;short_name;abbreviation
1;Tropical Rainforest;Tropical Rain.;TrRF
2;Tropical Seasonal & Deciduous Forest;Tropical Decid. Forest;TrDF
3;Temperate Broadleaved Evergreen Forest;Temp. Broad. Ever. Forest;TeBE
4;Temperate Broadleaved Deciduous Forest;Temp. Broad. Decid. Forest;TeBD
5;Mixed Forest;Mixed Forest;MF
6;Temperate Coniferous Forest;Temp. Conif. Forest;TeCF
7;Boreal Evergreen Forest;Bor. Ever. Forest;BoEF
8;Boreal Deciduous Forest;Bor. Decid. Forest;BoDF
9;Warm Woody Savanna, Woodland & Shrubland;Warm Woodland;WaWo
10;Warm Savanna & Open Shrubland;Warm Savanna;WaSa
11;Warm Grassland;Warm Grassland;WaGr
12;Temperate Woody Savanna, Woodland & Shrubland;Temp. Woodland;TeWo
13;Temperate Savanna & Open Shrubland;Temp. Savanna;TeSa
14;Temperate Grassland;Temp. Grassland;TeGr
15;Montane Grassland;Montane Grassland;MoGr
16;Arctic Tundra;Arctic Tundra;ArTu
17;Desert;Desert;Des
18;Rocks and Ice;Rocks and Ice;RoIc
19;Water;Water;Wat
This diff is collapsed.
UTF-8
\ No newline at end of file
File added
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]
\ No newline at end of file
File added
File added
File added
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]
\ No newline at end of file
File added
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment