biospheremetrics
The goal of biospheremetrics is to provide functions to calculate and plot the biosphere integrity metrics M-ECO and M-ECO in an R package based on outputs of LPJmL. biospheremetrics utilizes the read functions of the lpjmlkit package.
Installation
You can install biospheremetrics by git cloning this repository:
git clone https://gitlab.pik-potsdam.de/stenzel/biospheremetrics.git <path_to_biospheremetrics>and install via devtools:
devtools::install("<path_to_biospheremetrics>")
library("biospheremetrics")alternatively, you can also load it from source:
devtools::load_all("/p/projects/open/Fabian/LPJbox/biospheremetrics_paper/")Example
The ./scripts folder contains scripts to be used on the PIK cluster to
compute longer timeseries with higher RAM demand.
Example
The following application example calculates the metrics BioCol and EcoRisk:
library(devtools)
library(lpjmlkit)
library(sf)
library(terra)
devtools::load_all("/p/projects/open/Fabian/LPJbox/biospheremetrics_paper/")
run_folder <- "/p/projects/open/Fabian/runs/metrics_202306/output/lu_1500_2014/"
pnv_folder <- "/p/projects/open/Fabian/runs/metrics_202306/output/pnv_1500_2014/"
out_folder <- "/p/projects/open/Fabian/Metrics/"
lpj_input <- "/p/projects/lpjml/input/historical/"
# read grid
grid <- lpjmlkit::read_io(paste0(run_folder,"grid.bin.json"))
# calculate cell area
lat <- grid[, , 2]
lon <- grid[, , 1]
cellarea <- lpjmlkit::calc_cellarea(grid)
################# calculate BioCol ################
# 16GB of RAM are enough to calculate BioCol for a smaller analysis window (~40 years)
# for longer spans (500 years) - use separate script ("read_in_BioCol_data.R")
# and submit as cluster job using "sbatch R_read_in_BioCol_data.sh" - analysis for "biocol overtime" below
vars_biocol <- data.frame(
row.names = c("grid", "npp", "pft_npp", "pft_harvest", "pft_rharvest",
"firec", "timber_harvest", "cftfrac", "fpc"),
outname = c("grid.bin.json", "mnpp.bin.json", "pft_npp.bin.json",
"pft_harvest.pft.bin.json","pft_rharvest.pft.bin.json",
"firec.bin.json","timber_harvestc.bin.json","cftfrac.bin.json",
"fpc.bin.json")
)
biocol <- calc_biocol(
path_lu = run_folder,
path_pnv = pnv_folder,
gridbased = TRUE,
start_year = 1980,
stop_year = 2014,
reference_npp_time_span = 1510:1539,
reference_npp_file = "/p/projects/open/Fabian/runs/metrics_202306/output/pnv_1500_2014/mnpp.bin.json",
read_saved_data = FALSE,
save_data = TRUE,
npp_threshold = 20,
data_file = "/p/projects/open/Fabian/Metrics/BioCol_202306.RData",
external_fire = FALSE,
external_wood_harvest = TRUE,
external_fire_file = "/p/projects/open/Fabian/LPJbox/human_ignition_fraction.RData",
external_wood_harvest_file = "/p/projects/open/LanduseData/LUH2_v2h/wood_harvest_biomass_sum_1500-2014_67420.RData",
varnames = vars_biocol,
grass_scaling = FALSE,
include_fire = FALSE
)
plot_biocol(
biocol_data = biocol,
path_write = paste0(out_folder,"BioCol/"),
plotyears = c(1980,2014),
min_val = 0,
max_val = 90,
legendpos = "left",
start_year = 1980,
mapyear = 2000,
highlightyear = 2000,
eps = FALSE
)
############## analyse and plot biocol overtime #################
# first submit `R_read_BioCol_data.sh` to cluster via slurm to read in and process the input files, a lot of memory is required for this
# then here only read the preprocessed data file (read_saved_data = TRUE)
biocol_overtime <- calc_biocol(
path_lu = run_folder,
path_pnv = pnv_folder,
gridbased = TRUE,
start_year = 1500,
stop_year = 2014,
reference_npp_time_span = 1550:1579,
reference_npp_file = "/p/projects/open/Fabian/runs/metrics_202306/output/pnv_1500_2014/mnpp.bin.json",
read_saved_data = TRUE,
save_data = FALSE,
npp_threshold = 20,
data_file = "/p/projects/open/Fabian/Metrics/data/BioCol_202306_overtime.RData",
external_fire = FALSE,
external_wood_harvest = TRUE,
external_fire_file = "/p/projects/open/Fabian/LPJbox/human_ignition_fraction.RData",
external_wood_harvest_file = "/p/projects/open/LanduseData/LUH2_v2h/wood_harvest_biomass_sum_1500-2014_67420.RData",
varnames = vars_biocol,
grass_scaling = FALSE,
include_fire = FALSE
)
plot_biocol(
biocol_data = biocol_overtime,
path_write = paste0(out_folder,"BioCol/"),
plotyears = c(1550,2014),
min_val = 0,
max_val = 90,
legendpos = list(x=1550,y=23),
start_year = 1500,
mapyear = 2000,
highlightyear = 2000,
eps = FALSE
)
################# compute EcoRisk ################
vars_ecorisk <- data.frame(
row.names = c("grid","fpc", "fpc_bft", "cftfrac", "firec", "npp", "runoff",
"transp", "vegc", "firef", "rh", "harvestc", "rharvestc",
"pft_harvestc", "pft_rharvestc", "evap", "interc", "discharge",
"soilc", "litc", "swc", "vegn", "soilnh4", "soilno3",
"leaching", "n2o_denit", "n2o_nit", "n2_emis", "bnf",
"n_volatilization"),
outname = c("grid.bin.json", "fpc.bin.json", "fpc_bft.bin.json",
"cftfrac.bin.json", "firec.bin.json", "mnpp.bin.json",
"mrunoff.bin.json", "mtransp.bin.json", "vegc.bin.json",
"firef.bin.json", "mrh.bin.json", "flux_harvest.bin.json",
"flux_rharvest.bin.json", "pft_harvest.pft.bin.json",
"pft_rharvest.pft.bin.json", "mevap.bin.json",
"minterc.bin.json", "mdischarge.bin.json", "soilc.bin.json",
"litc.bin.json", "mswc.bin.json", "vegn.bin.json",
"soilnh4.bin.json", "soilno3.bin.json", "mleaching.bin.json",
"mn2o_denit.bin.json", "mn2o_nit.bin.json", "mn2_emis.bin.json",
"mbnf.bin.json", "mn_volatilization.bin.json")
)
ecorisk <- ecorisk_wrapper(
path_ref = pnv_folder,
path_scen = run_folder,
read_saved_data = TRUE,
nitrogen = TRUE,
varnames = vars_ecorisk,
weighting = "equal",
save_data = "/p/projects/open/Fabian/Metrics/data/ecorisk_202306_data.RData",
save_ecorisk = "/p/projects/open/Fabian/Metrics/data/ecorisk_202306_gamma.RData",
time_span_reference = c(1550:1579),
time_span_scenario = c(1985:2014),
dimensions_only_local = FALSE
)
# plot ecorisk
plot_ecorisk_map(
ecorisk$ecorisk_total,
file = paste0(out_folder,"EcoRisk/ecorisk.png"),
title="ecorisk"
)
plot_ecorisk_map(
ecorisk$vegetation_structure_change,
file = paste0(out_folder, "EcoRisk/vs.png"),
title = "vegetation structure change"
)
plot_ecorisk_map(
ecorisk$local_change,
file = paste0(out_folder, "EcoRisk/lc.png"),
title = "local change"
)
plot_ecorisk_map(
ecorisk$global_importance,
file = paste0(out_folder, "EcoRisk/gi.png"),
title = "global importance"
)
plot_ecorisk_map(
ecorisk$ecosystem_balance,
file = paste0(out_folder, "EcoRisk/eb.png"),
title = "ecosystem balance")
plot_ecorisk_map(
ecorisk$carbon_stocks,
file = paste0(out_folder, "EcoRisk/cs.png"),
title = "carbon_stocks"
)
plot_ecorisk_map(
ecorisk$carbon_fluxes,
file = paste0(out_folder, "EcoRisk/cf.png"),
title = "carbon_fluxes"
)
plot_ecorisk_map(
ecorisk$water_stocks,
file = paste0(out_folder, "EcoRisk/ws.png"),
title = " water_stocks"
)
plot_ecorisk_map(
ecorisk$water_fluxes,
file = paste0(out_folder, "EcoRisk/wf.png"),
title = " water_fluxes"
)
plot_ecorisk_map(
ecorisk$nitrogen_stocks,
file = paste0(out_folder, "EcoRisk/ns.png"),
title = " nitrogen_stocks"
)
plot_ecorisk_map(
ecorisk$nitrogen_fluxes,
file = paste0(out_folder, "EcoRisk/nf.png"),
title = " nitrogen_fluxes"
)
################# ecorisk biomes ################
biome_classes <- classify_biomes(
path_reference = pnv_folder,
files_reference = list(
grid = paste0(pnv_folder,"grid.bin.json"),
fpc = paste0(pnv_folder,"fpc.bin.json"),
vegc = paste0(pnv_folder,"vegc.bin.json"),
pft_lai = paste0(pnv_folder,"pft_lai.bin.json"),
temp = "/p/projects/lpjml/input/historical/GSWP3-W5E5/tas_gswp3-w5e5_1901-2016.clm",
elevation = "/p/projects/lpjml/input/historical/input_VERSION2/elevation.bin"
),
time_span_reference = as.character(1985:2014),
savanna_proxy = list(pft_lai = 6),
montane_arctic_proxy = list(elevation = 1000)
)
biome_classes_pi <- classify_biomes(
path_reference = pnv_folder,
files_reference = list(
grid = paste0(pnv_folder,"grid.bin.json"),
fpc = paste0(pnv_folder,"fpc.bin.json"),
vegc = paste0(pnv_folder,"vegc.bin.json"),
pft_lai = paste0(pnv_folder,"pft_lai.bin.json"),
temp = "/p/projects/lpjml/input/historical/GSWP3-W5E5/tas_gswp3-w5e5_1901-2016.clm",
elevation = "/p/projects/lpjml/input/historical/input_VERSION2/elevation.bin"
),
time_span_reference = as.character(1901:1910),
savanna_proxy = list(pft_lai = 6),
montane_arctic_proxy = list(elevation = 1000)
)
plot_biomes(biome_data = biome_classes,
display_area = TRUE,
cellarea = cellarea,
file_name = paste0(out_folder,"EcoRisk/biomes_2005-2014.png"),
order_legend = 1:19,
to_robinson = FALSE)
plot_biomes(biome_data=biome_classes_pi,
display_area = TRUE,
cellarea = cellarea,
file_name = paste0(out_folder,"EcoRisk/biomes_1901-1910.png"),
order_legend = 1:19,
to_robinson = FALSE)
# compute median ecorisk values for biomes/large worldregions
ecorisk_disaggregated_full <- disaggregate_into_biomes(
data = ecorisk,
biome_class = biome_classes,
type = "quantile",
classes = "allbiomes"
)
ecorisk_disaggregated_full[is.na(meco_disaggregated_full)] <- 0
ecorisk_disaggregated_4regions <- disaggregate_into_biomes(
data = ecorisk,
biome_class = biome_classes,
type = "quantile",
classes = "4biomes"
)
plot_ecorisk_radial_panel(
data = ecorisk_disaggregated_full[-c(17,18,19),,],
biomeNames = get_biome_names(1)[-c(17,18,19)],
file = paste0(out_folder,"EcoRisk/EcoRisk_panel_1564_vs_2002.png"),
quantile = TRUE,
eps = TRUE
)
plot_ecorisk_radial_panel(
data = ecorisk_disaggregated_4regions[,,],
biomeNames = c("tropics","temperate","boreal","arctic"),
file = paste0(out_folder,"EcoRisk/EcoRisk_4regions_1564_vs_2002.png"),
quantile = TRUE,
eps = TRUE
)
################# ecorisk overtime ################
# first use the script `R_calc_ecorisk_overtime.sh` to read in and process the data
# on the PIK cluster this takes about a day for 100 years and 80GB of memory
load("/p/projects/open/Fabian/Metrics/data/ecorisk_202306_overtime_gamma.RData")
ecorisk_overtime_allbiomes <- disaggregate_into_biomes(
ecorisk = ecorisk,
biome_class = biome_classes,
type = "quantile",
classes = "allbiomes"
)
plot_ecorisk_over_time_panel(
data = ecorisk_overtime_allbiomes,
timerange = c(1916,2003),
biomeNames = c("tropic","temperate","boreal","arctic"),
file = paste0(out_folder,"overtime_panel.png"),
eps=TRUE
)
ecorisk_overtime_biome16 <- disaggregate_into_biomes(
data = ecorisk,
biome_class = biome_classes,
type = "quantile",
classes = "allbiomes"
)
plot_ecorisk_over_time_panel(
data = ecorisk_overtime_biome16[-c(3,17,18),,,],
timerange = c(1916,2003),
biomeNames = get_biome_names(1)[-c(3,17,18)],
file = paste0(out_folder,"overtime_panel_16.png"),
eps=TRUE
)
################# compare to average PI biome cell #################
intra_biome_distrib_PI <- calculate_within_biome_diffs(
biome_classes = biome_classes_pi,
intra_biome_distrib_file = "/p/projects/open/Fabian/Metrics/data/ecorisk_PNV_intra_biome_distrib_file_202306.RData",
dataFile_base = "/p/projects/open/Fabian/Metrics/data/ecorisk_202306_data.RData",
create = TRUE, plotting = TRUE, res = 0.02, vars_ecorisk = vars_ecorisk,
plot_folder = out_folder, time_span_reference = as.character(1891:1920))
plot_biome_internal_distribution(
data = intra_biome_distrib_PI[,"ecorisk_total",],
file = paste0(out_folder,"EcoRisk_newCol/distribution_PI_within_biome_differences.png"),
biomes_abbrv = get_biome_names(1),
scale = 4,
eps=TRUE,
palette = paletteNew
)
################# cross table average biomes today #################
dataFile_base = "/p/projects/open/Fabian/Metrics/data/ecorisk_202306_data.RData"
data_file = "/p/projects/open/Fabian/Metrics/data/ecorisk_202306_crosstable_data.RData"
ecoriskFile = "/p/projects/open/Fabian/Metrics/data/ecorisk_202306_crosstable_gamma.RData"
ecorisk_cross_table(dataFileIn = dataFile_base,
dataFileOut = data_file,
biome_classes_in = biome_classes) #pickCells = pickcells)
nbiomes <- length(biome_classes$biome_names)
ecorisk_crosstable_today <- ecorisk_wrapper(
path_ref = NULL,
path_scen = NULL,
read_saved_data = TRUE,
save_data = data_file,
save_ecorisk = ecoriskFile,
varnames = vars_ecorisk,
time_span_reference = as.character(1985:2014),
time_span_scenario = as.character(1985:2014)
#ncells = nbiomes^2
)
# if written previously, load crosstable data
if (FALSE) {
load(ecoriskFile)
ecorisk_crosstable_today <- ecorisk
}
crosstable <- ecorisk_crosstable_today$ecorisk_total
dim(crosstable) <- c(nbiomes,nbiomes)
colnames(crosstable) <- get_biome_names(1)
rownames(crosstable) <- get_biome_names(2)
plot_ecorisk_cross_table(
data = crosstable[-c(3,8,18,19),-c(3,8,18,19)],
file = paste0(out_folder,"/EcoRisk_newCol/crosstable_today.png"),
lmar=12,
palette = paletteNew
)