Computing biospheremetrics

biospheremetrics is an R package to calculate two complementary terrestrial biosphere integrity indicators: * human colonization of the biosphere (BioCol) * risk of ecosystem destabilization (EcoRisk)

based on outputs from simulations with a dynamic vegetation model (currently only LPJmL https://github.com/PIK-LPJmL/LPJmL). biospheremetrics utilizes the read functions of the lpjmlkit (https://github.com/PIK-LPJmL/lpjmlkit) package.

Installation lpjmlkit

The dependency package lpjmlkit is currently not available from CRAN, therefore it needs to be installed manually.

Download from github:

git clone [email protected]:PIK-LPJmL/lpjmlkit.git <path_to_lpjmlkit>

and install via devtools:

devtools::install("<path_to_lpjmlkit>")
library("lpjmlkit")

alternatively, you can also just load it from source:

devtools::load_all("<path_to_lpjmlkit>")

Installation biospheremetrics

You can install biospheremetrics by git cloning this repository:

git clone [email protected]:stenzelf/biospheremetrics.git <path_to_biospheremetrics>

and install via devtools:

devtools::install("<path_to_biospheremetrics>")
library("biospheremetrics")

alternatively, you can also just load it from source:

devtools::load_all("<path_to_biospheremetrics>")

Preparations

In order to compute the indicators it is necessary to first run LPJmL, and produce the required outputs. A list of required outputs per metric (e.g. for EcoRisk with nitrogen) can be obtained like this:

biospheremetrics::list_outputs(metric = "ecorisk_nitrogen")
#> $grid
#> $grid$resolution
#> [1] "annual"
#> 
#> $grid$optional
#> [1] FALSE
#> 
#> $grid$file_name
#> [1] "grid"
#> 
#> 
#> $terr_area
#> $terr_area$resolution
#> [1] "annual"
#> 
#> $terr_area$optional
#> [1] FALSE
#> 
#> $terr_area$file_name
#> [1] "terr_area"
#> 
#> 
#> $fpc
#> $fpc$resolution
#> [1] "annual"
#> 
#> $fpc$optional
#> [1] FALSE
#> 
#> $fpc$file_name
#> [1] "fpc"
#> 
#> 
#> $fpc_bft
#> $fpc_bft$resolution
#> [1] "annual"
#> 
#> $fpc_bft$optional
#> [1] FALSE
#> 
#> $fpc_bft$file_name
#> [1] "fpc_bft"
#> 
#> 
#> $cftfrac
#> $cftfrac$resolution
#> [1] "annual"
#> 
#> $cftfrac$optional
#> [1] FALSE
#> 
#> $cftfrac$file_name
#> [1] "cftfrac"
#> 
#> 
#> $firec
#> $firec$resolution
#> [1] "annual"
#> 
#> $firec$optional
#> [1] FALSE
#> 
#> $firec$file_name
#> [1] "firec"
#> 
#> 
#> $npp
#> $npp$resolution
#> [1] "annual"
#> 
#> $npp$optional
#> [1] FALSE
#> 
#> $npp$file_name
#> [1] "npp"
#> 
#> 
#> $gpp
#> $gpp$resolution
#> [1] "annual"
#> 
#> $gpp$optional
#> [1] FALSE
#> 
#> $gpp$file_name
#> [1] "gpp"
#> 
#> 
#> $runoff
#> $runoff$resolution
#> [1] "annual"
#> 
#> $runoff$optional
#> [1] FALSE
#> 
#> $runoff$file_name
#> [1] "runoff"
#> 
#> 
#> $irrig
#> $irrig$resolution
#> [1] "annual"
#> 
#> $irrig$optional
#> [1] FALSE
#> 
#> $irrig$file_name
#> [1] "irrig"
#> 
#> 
#> $transp
#> $transp$resolution
#> [1] "annual"
#> 
#> $transp$optional
#> [1] FALSE
#> 
#> $transp$file_name
#> [1] "transp"
#> 
#> 
#> $vegc
#> $vegc$resolution
#> [1] "annual"
#> 
#> $vegc$optional
#> [1] FALSE
#> 
#> $vegc$file_name
#> [1] "vegc"
#> 
#> 
#> $firef
#> $firef$resolution
#> [1] "annual"
#> 
#> $firef$optional
#> [1] FALSE
#> 
#> $firef$file_name
#> [1] "firef"
#> 
#> 
#> $harvestc
#> $harvestc$resolution
#> [1] "annual"
#> 
#> $harvestc$optional
#> [1] FALSE
#> 
#> $harvestc$file_name
#> [1] "harvestc"
#> 
#> 
#> $timber_harvestc
#> $timber_harvestc$resolution
#> [1] "annual"
#> 
#> $timber_harvestc$optional
#> [1] FALSE
#> 
#> $timber_harvestc$file_name
#> [1] "timber_harvestc"
#> 
#> 
#> $evap
#> $evap$resolution
#> [1] "annual"
#> 
#> $evap$optional
#> [1] FALSE
#> 
#> $evap$file_name
#> [1] "evap"
#> 
#> 
#> $interc
#> $interc$resolution
#> [1] "annual"
#> 
#> $interc$optional
#> [1] FALSE
#> 
#> $interc$file_name
#> [1] "interc"
#> 
#> 
#> $soilc
#> $soilc$resolution
#> [1] "annual"
#> 
#> $soilc$optional
#> [1] FALSE
#> 
#> $soilc$file_name
#> [1] "soilc"
#> 
#> 
#> $litc
#> $litc$resolution
#> [1] "annual"
#> 
#> $litc$optional
#> [1] FALSE
#> 
#> $litc$file_name
#> [1] "litc"
#> 
#> 
#> $soilwater
#> $soilwater$resolution
#> [1] "annual"
#> 
#> $soilwater$optional
#> [1] FALSE
#> 
#> $soilwater$file_name
#> [1] "soilwater"
#> 
#> 
#> $res_storage
#> $res_storage$resolution
#> [1] "annual"
#> 
#> $res_storage$optional
#> [1] FALSE
#> 
#> $res_storage$file_name
#> [1] "res_storage"
#> 
#> 
#> $irrig_stor
#> $irrig_stor$resolution
#> [1] "annual"
#> 
#> $irrig_stor$optional
#> [1] FALSE
#> 
#> $irrig_stor$file_name
#> [1] "irrig_stor"
#> 
#> 
#> $rivervol
#> $rivervol$resolution
#> [1] "annual"
#> 
#> $rivervol$optional
#> [1] FALSE
#> 
#> $rivervol$file_name
#> [1] "rivervol"
#> 
#> 
#> $lakevol
#> $lakevol$resolution
#> [1] "annual"
#> 
#> $lakevol$optional
#> [1] FALSE
#> 
#> $lakevol$file_name
#> [1] "lakevol"
#> 
#> 
#> $prec
#> $prec$resolution
#> [1] "annual"
#> 
#> $prec$optional
#> [1] FALSE
#> 
#> $prec$file_name
#> [1] "prec"
#> 
#> 
#> $vegn
#> $vegn$resolution
#> [1] "annual"
#> 
#> $vegn$optional
#> [1] FALSE
#> 
#> $vegn$file_name
#> [1] "vegn"
#> 
#> 
#> $harvestn
#> $harvestn$resolution
#> [1] "annual"
#> 
#> $harvestn$optional
#> [1] FALSE
#> 
#> $harvestn$file_name
#> [1] "harvestn"
#> 
#> 
#> $firen
#> $firen$resolution
#> [1] "annual"
#> 
#> $firen$optional
#> [1] FALSE
#> 
#> $firen$file_name
#> [1] "firen"
#> 
#> 
#> $soilnh4
#> $soilnh4$resolution
#> [1] "annual"
#> 
#> $soilnh4$optional
#> [1] FALSE
#> 
#> $soilnh4$file_name
#> [1] "soilnh4"
#> 
#> 
#> $soilno3
#> $soilno3$resolution
#> [1] "annual"
#> 
#> $soilno3$optional
#> [1] FALSE
#> 
#> $soilno3$file_name
#> [1] "soilno3"
#> 
#> 
#> $leaching
#> $leaching$resolution
#> [1] "annual"
#> 
#> $leaching$optional
#> [1] FALSE
#> 
#> $leaching$file_name
#> [1] "leaching"
#> 
#> 
#> $n2o_denit
#> $n2o_denit$resolution
#> [1] "annual"
#> 
#> $n2o_denit$optional
#> [1] FALSE
#> 
#> $n2o_denit$file_name
#> [1] "n2o_denit"
#> 
#> 
#> $n2o_nit
#> $n2o_nit$resolution
#> [1] "annual"
#> 
#> $n2o_nit$optional
#> [1] FALSE
#> 
#> $n2o_nit$file_name
#> [1] "n2o_nit"
#> 
#> 
#> $n2_emis
#> $n2_emis$resolution
#> [1] "annual"
#> 
#> $n2_emis$optional
#> [1] FALSE
#> 
#> $n2_emis$file_name
#> [1] "n2_emis"
#> 
#> 
#> $bnf
#> $bnf$resolution
#> [1] "annual"
#> 
#> $bnf$optional
#> [1] FALSE
#> 
#> $bnf$file_name
#> [1] "bnf"
#> 
#> 
#> $ndepos
#> $ndepos$resolution
#> [1] "annual"
#> 
#> $ndepos$optional
#> [1] FALSE
#> 
#> $ndepos$file_name
#> [1] "ndepos"
#> 
#> 
#> $n_volatilization
#> $n_volatilization$resolution
#> [1] "annual"
#> 
#> $n_volatilization$optional
#> [1] FALSE
#> 
#> $n_volatilization$file_name
#> [1] "n_volatilization"
#> 
#> 
#> $nfert_agr
#> $nfert_agr$resolution
#> [1] "annual"
#> 
#> $nfert_agr$optional
#> [1] FALSE
#> 
#> $nfert_agr$file_name
#> [1] "nfert_agr"
#> 
#> 
#> $nmanure_agr
#> $nmanure_agr$resolution
#> [1] "annual"
#> 
#> $nmanure_agr$optional
#> [1] FALSE
#> 
#> $nmanure_agr$file_name
#> [1] "nmanure_agr"
biospheremetrics::list_outputs(metric = "ecorisk_nitrogen")
#> $grid
#> $grid$resolution
#> [1] "annual"
#> 
#> $grid$optional
#> [1] FALSE
#> 
#> $grid$file_name
#> [1] "grid"
#> 
#> 
#> $terr_area
#> $terr_area$resolution
#> [1] "annual"
#> 
#> $terr_area$optional
#> [1] FALSE
#> 
#> $terr_area$file_name
#> [1] "terr_area"
#> 
#> 
#> $fpc
#> $fpc$resolution
#> [1] "annual"
#> 
#> $fpc$optional
#> [1] FALSE
#> 
#> $fpc$file_name
#> [1] "fpc"
#> 
#> 
#> $fpc_bft
#> $fpc_bft$resolution
#> [1] "annual"
#> 
#> $fpc_bft$optional
#> [1] FALSE
#> 
#> $fpc_bft$file_name
#> [1] "fpc_bft"
#> 
#> 
#> $cftfrac
#> $cftfrac$resolution
#> [1] "annual"
#> 
#> $cftfrac$optional
#> [1] FALSE
#> 
#> $cftfrac$file_name
#> [1] "cftfrac"
#> 
#> 
#> $firec
#> $firec$resolution
#> [1] "annual"
#> 
#> $firec$optional
#> [1] FALSE
#> 
#> $firec$file_name
#> [1] "firec"
#> 
#> 
#> $npp
#> $npp$resolution
#> [1] "annual"
#> 
#> $npp$optional
#> [1] FALSE
#> 
#> $npp$file_name
#> [1] "npp"
#> 
#> 
#> $gpp
#> $gpp$resolution
#> [1] "annual"
#> 
#> $gpp$optional
#> [1] FALSE
#> 
#> $gpp$file_name
#> [1] "gpp"
#> 
#> 
#> $runoff
#> $runoff$resolution
#> [1] "annual"
#> 
#> $runoff$optional
#> [1] FALSE
#> 
#> $runoff$file_name
#> [1] "runoff"
#> 
#> 
#> $irrig
#> $irrig$resolution
#> [1] "annual"
#> 
#> $irrig$optional
#> [1] FALSE
#> 
#> $irrig$file_name
#> [1] "irrig"
#> 
#> 
#> $transp
#> $transp$resolution
#> [1] "annual"
#> 
#> $transp$optional
#> [1] FALSE
#> 
#> $transp$file_name
#> [1] "transp"
#> 
#> 
#> $vegc
#> $vegc$resolution
#> [1] "annual"
#> 
#> $vegc$optional
#> [1] FALSE
#> 
#> $vegc$file_name
#> [1] "vegc"
#> 
#> 
#> $firef
#> $firef$resolution
#> [1] "annual"
#> 
#> $firef$optional
#> [1] FALSE
#> 
#> $firef$file_name
#> [1] "firef"
#> 
#> 
#> $harvestc
#> $harvestc$resolution
#> [1] "annual"
#> 
#> $harvestc$optional
#> [1] FALSE
#> 
#> $harvestc$file_name
#> [1] "harvestc"
#> 
#> 
#> $timber_harvestc
#> $timber_harvestc$resolution
#> [1] "annual"
#> 
#> $timber_harvestc$optional
#> [1] FALSE
#> 
#> $timber_harvestc$file_name
#> [1] "timber_harvestc"
#> 
#> 
#> $evap
#> $evap$resolution
#> [1] "annual"
#> 
#> $evap$optional
#> [1] FALSE
#> 
#> $evap$file_name
#> [1] "evap"
#> 
#> 
#> $interc
#> $interc$resolution
#> [1] "annual"
#> 
#> $interc$optional
#> [1] FALSE
#> 
#> $interc$file_name
#> [1] "interc"
#> 
#> 
#> $soilc
#> $soilc$resolution
#> [1] "annual"
#> 
#> $soilc$optional
#> [1] FALSE
#> 
#> $soilc$file_name
#> [1] "soilc"
#> 
#> 
#> $litc
#> $litc$resolution
#> [1] "annual"
#> 
#> $litc$optional
#> [1] FALSE
#> 
#> $litc$file_name
#> [1] "litc"
#> 
#> 
#> $soilwater
#> $soilwater$resolution
#> [1] "annual"
#> 
#> $soilwater$optional
#> [1] FALSE
#> 
#> $soilwater$file_name
#> [1] "soilwater"
#> 
#> 
#> $res_storage
#> $res_storage$resolution
#> [1] "annual"
#> 
#> $res_storage$optional
#> [1] FALSE
#> 
#> $res_storage$file_name
#> [1] "res_storage"
#> 
#> 
#> $irrig_stor
#> $irrig_stor$resolution
#> [1] "annual"
#> 
#> $irrig_stor$optional
#> [1] FALSE
#> 
#> $irrig_stor$file_name
#> [1] "irrig_stor"
#> 
#> 
#> $rivervol
#> $rivervol$resolution
#> [1] "annual"
#> 
#> $rivervol$optional
#> [1] FALSE
#> 
#> $rivervol$file_name
#> [1] "rivervol"
#> 
#> 
#> $lakevol
#> $lakevol$resolution
#> [1] "annual"
#> 
#> $lakevol$optional
#> [1] FALSE
#> 
#> $lakevol$file_name
#> [1] "lakevol"
#> 
#> 
#> $prec
#> $prec$resolution
#> [1] "annual"
#> 
#> $prec$optional
#> [1] FALSE
#> 
#> $prec$file_name
#> [1] "prec"
#> 
#> 
#> $vegn
#> $vegn$resolution
#> [1] "annual"
#> 
#> $vegn$optional
#> [1] FALSE
#> 
#> $vegn$file_name
#> [1] "vegn"
#> 
#> 
#> $harvestn
#> $harvestn$resolution
#> [1] "annual"
#> 
#> $harvestn$optional
#> [1] FALSE
#> 
#> $harvestn$file_name
#> [1] "harvestn"
#> 
#> 
#> $firen
#> $firen$resolution
#> [1] "annual"
#> 
#> $firen$optional
#> [1] FALSE
#> 
#> $firen$file_name
#> [1] "firen"
#> 
#> 
#> $soilnh4
#> $soilnh4$resolution
#> [1] "annual"
#> 
#> $soilnh4$optional
#> [1] FALSE
#> 
#> $soilnh4$file_name
#> [1] "soilnh4"
#> 
#> 
#> $soilno3
#> $soilno3$resolution
#> [1] "annual"
#> 
#> $soilno3$optional
#> [1] FALSE
#> 
#> $soilno3$file_name
#> [1] "soilno3"
#> 
#> 
#> $leaching
#> $leaching$resolution
#> [1] "annual"
#> 
#> $leaching$optional
#> [1] FALSE
#> 
#> $leaching$file_name
#> [1] "leaching"
#> 
#> 
#> $n2o_denit
#> $n2o_denit$resolution
#> [1] "annual"
#> 
#> $n2o_denit$optional
#> [1] FALSE
#> 
#> $n2o_denit$file_name
#> [1] "n2o_denit"
#> 
#> 
#> $n2o_nit
#> $n2o_nit$resolution
#> [1] "annual"
#> 
#> $n2o_nit$optional
#> [1] FALSE
#> 
#> $n2o_nit$file_name
#> [1] "n2o_nit"
#> 
#> 
#> $n2_emis
#> $n2_emis$resolution
#> [1] "annual"
#> 
#> $n2_emis$optional
#> [1] FALSE
#> 
#> $n2_emis$file_name
#> [1] "n2_emis"
#> 
#> 
#> $bnf
#> $bnf$resolution
#> [1] "annual"
#> 
#> $bnf$optional
#> [1] FALSE
#> 
#> $bnf$file_name
#> [1] "bnf"
#> 
#> 
#> $ndepos
#> $ndepos$resolution
#> [1] "annual"
#> 
#> $ndepos$optional
#> [1] FALSE
#> 
#> $ndepos$file_name
#> [1] "ndepos"
#> 
#> 
#> $n_volatilization
#> $n_volatilization$resolution
#> [1] "annual"
#> 
#> $n_volatilization$optional
#> [1] FALSE
#> 
#> $n_volatilization$file_name
#> [1] "n_volatilization"
#> 
#> 
#> $nfert_agr
#> $nfert_agr$resolution
#> [1] "annual"
#> 
#> $nfert_agr$optional
#> [1] FALSE
#> 
#> $nfert_agr$file_name
#> [1] "nfert_agr"
#> 
#> 
#> $nmanure_agr
#> $nmanure_agr$resolution
#> [1] "annual"
#> 
#> $nmanure_agr$optional
#> [1] FALSE
#> 
#> $nmanure_agr$file_name
#> [1] "nmanure_agr"

We have however included some demo data, to check out some functionality without access to LPJmL.

Demo case

The package contains some demo LPJmL simulation data for two cells. Define the folders for run (simulation with climate change and land-use) and pnv (potential natural vegetation - climate change, but no land-use) and a temporary output folder, if you want to test to write to files.

run_folder <- paste0(system.file("extdata","run","lu_1500_2016",package = "biospheremetrics"),"/")
pnv_folder <- paste0(system.file("extdata","run","pnv_1500_2016",package = "biospheremetrics"),"/")
out_folder <- paste0(tempdir(),"/")

BioCol

First let’s have a look at the indicator BioCol and compute it based on the demo data folders:

biocol <- biospheremetrics::calc_biocol(
  path_lu = run_folder,
  path_pnv = pnv_folder,
  gridbased = TRUE,
  start_year = 1500,
  stop_year = 2016,
  reference_npp_time_span = 1510:1539,
  read_saved_data = FALSE,
  save_data = FALSE
)
#> Reading in data from outputs
#> Calculating data

Then let’s plot the progression of the sum of the data over time (for an example with global simulation data see the supplementary manuscript):

biospheremetrics::plot_biocol_ts(
  biocol_data = biocol,
  plot_years = c(1510,2016),
  first_year = 1510,
  min_val = 0,
  max_val = 0.005,
  legendpos = "left",
  highlight_years = NA,
  eps = FALSE
)

EcoRisk

Let’s have a look at the EcoRist indicator next and compute it based on the folders with demo data:

ecorisk <- biospheremetrics::ecorisk_wrapper(
  path_ref = pnv_folder,
  path_scen = run_folder,
  read_saved_data = FALSE,
  nitrogen = TRUE,
  save_data = NULL,
  save_ecorisk = NULL,
  time_span_reference = c(1550:1579),
  time_span_scenario = c(1987:2016)
  )
#> Reading in fpc, fpc_bft, cftfrac
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/vegc.bin.json with unit gC/m2 -> as part of vegetation_carbon_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/vegc.bin.json with unit gC/m2 -> as part of vegetation_carbon_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/soilc.bin.json with unit gC/m2 -> as part of soil_carbon_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/soilc.bin.json with unit gC/m2 -> as part of soil_carbon_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/litc.bin.json with unit gC/m2 -> as part of soil_carbon_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/litc.bin.json with unit gC/m2 -> as part of soil_carbon_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/gpp.bin.json with unit gC/m2/month -> as part of carbon_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/gpp.bin.json with unit gC/m2/month -> as part of carbon_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/gpp.bin.json with unit gC/m2/month -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/gpp.bin.json with unit gC/m2/month -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/npp.bin.json with unit gC/m2/month -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/npp.bin.json with unit gC/m2/month -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/firec.bin.json with unit gC/m2/yr -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/firec.bin.json with unit gC/m2/yr -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/harvestc.bin.json with unit gC/m2/yr -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/harvestc.bin.json with unit gC/m2/yr -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/timber_harvestc.bin.json with unit gC/m2/yr -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/timber_harvestc.bin.json with unit gC/m2/yr -> as part of carbon_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/soilwater.bin.json with unit dm3 -> as part of soil_water_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/soilwater.bin.json with unit dm3 -> as part of soil_water_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/prec.bin.json with unit mm/month -> as part of water_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/prec.bin.json with unit mm/month -> as part of water_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/irrig.bin.json with unit mm/month -> as part of water_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/irrig.bin.json with unit mm/month -> as part of water_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/transp.bin.json with unit mm/month -> as part of water_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/transp.bin.json with unit mm/month -> as part of water_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/evap.bin.json with unit mm/month -> as part of water_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/evap.bin.json with unit mm/month -> as part of water_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/interc.bin.json with unit mm/month -> as part of water_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/interc.bin.json with unit mm/month -> as part of water_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/runoff.bin.json with unit mm/month -> as part of water_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/runoff.bin.json with unit mm/month -> as part of water_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/firef.bin.json with unit yr -> as part of other
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/firef.bin.json with unit yr -> as part of other
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/vegn.bin.json with unit gN/m2 -> as part of vegetation_nitrogen_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/vegn.bin.json with unit gN/m2 -> as part of vegetation_nitrogen_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/soilnh4.bin.json with unit gN/m2 -> as part of soil_mineral_nitrogen_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/soilnh4.bin.json with unit gN/m2 -> as part of soil_mineral_nitrogen_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/soilno3.bin.json with unit gN/m2 -> as part of soil_mineral_nitrogen_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/soilno3.bin.json with unit gN/m2 -> as part of soil_mineral_nitrogen_pool
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/bnf.bin.json with unit gN/m2/month -> as part of nitrogen_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/bnf.bin.json with unit gN/m2/month -> as part of nitrogen_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/nfert_agr.bin.json with unit gN/m2/yr -> as part of nitrogen_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/nfert_agr.bin.json with unit gN/m2/yr -> as part of nitrogen_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/nmanure_agr.bin.json with unit gN/m2/yr -> as part of nitrogen_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/nmanure_agr.bin.json with unit gN/m2/yr -> as part of nitrogen_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/ndepos.bin.json with unit gN/m2/yr -> as part of nitrogen_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/ndepos.bin.json with unit gN/m2/yr -> as part of nitrogen_influx
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/harvestn.bin.json with unit gN/m2/yr -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/harvestn.bin.json with unit gN/m2/yr -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/firen.bin.json with unit gN/m2/yr -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/firen.bin.json with unit gN/m2/yr -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/leaching.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/leaching.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/n2_emis.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/n2_emis.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/n2o_denit.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/n2o_denit.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/n2o_nit.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/n2o_nit.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/lu_1500_2016/n_volatilization.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Reading in /tmp/RtmpDh3ag6/Rinst1189562f6b20/biospheremetrics/extdata/run/pnv_1500_2016/n_volatilization.bin.json with unit gN/m2/month -> as part of nitrogen_outflux
#> Calculating time slice 2002 of 2002-2002

And look at how the aggregated ecorisk_total indicator and the subcomponents vs, lc, gi and eb are evaluated:

ecorisk$ecorisk_total
#>     year
#> cell      2002
#>    0 0.7691418
#>    1 0.9103067
ecorisk$vegetation_structure_change
#>     year
#> cell      2002
#>    0 0.2409668
#>    1 0.6412267
ecorisk$local_change
#>     year
#> cell 2002
#>    0    1
#>    1    1
ecorisk$global_importance
#>     year
#> cell 2002
#>    0    1
#>    1    1
ecorisk$ecosystem_balance
#>     year
#> cell      2002
#>    0 0.8356003
#>    1 1.0000000

Metrics definitions and file names

In the file ./inst/extdata/metric_files.yml the filenames and how they make up certain metric components (mostly EcoRisk) are defined. Usually you should not need to change anything here, start your LPJmL run with the required outputs already with the right names, using the list_outputs() function.

If you need to change sth here, please be careful. To change definitions, you currently need clone the package and modify the ./inst/extdata/metric_files.yml file locally. Afterwards reload/compile the package.

Example use

The following application example calculates the metrics BioCol and EcoRisk from a global run:

run_folder <- "./output/lu_1500_2014/"
pnv_folder <- "./output/pnv_1500_2014/"
out_folder <- "./biospheremetrics/"
lpj_input <- "/path/to/historical/lpjml/input/data/"

################# 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

biocol <- biospheremetrics::calc_biocol(
  path_lu = run_folder,
  path_pnv = pnv_folder,
  gridbased = TRUE,
  start_year = 1980,
  stop_year = 2014,
  reference_npp_time_span = 1510:1539, 
  read_saved_data = FALSE,
  save_data = FALSE,
  npp_threshold = 20,
)

biospheremetrics::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)
biospheremetrics::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,
  read_saved_data = TRUE,
  save_data = FALSE,
  npp_threshold = 20,
  data_file = "</path/to/saved/data/file>" # from script read_in_biocol_data
)

biospheremetrics::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 ################

ecorisk <- biospheremetrics::ecorisk_wrapper(
  path_ref = pnv_folder, 
  path_scen = run_folder, 
  read_saved_data = FALSE,
  nitrogen = TRUE,
  weighting = "equal",
  save_data = NULL,
  save_ecorisk = NULL,
  time_span_reference = c(1550:1579),
  time_span_scenario = c(1985:2014),
  dimensions_only_local = FALSE
)

# plot ecorisk
biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "ecorisk_total",
  file = paste0(out_folder,"EcoRisk/ecorisk.png"),
  title = "ecorisk"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "vegetation_structure_change",
  file = paste0(out_folder, "EcoRisk/vs.png"),
  title = "vegetation structure change"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "local_change",
  file = paste0(out_folder, "EcoRisk/lc.png"),
  title = "local change"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "global_importance",
  file = paste0(out_folder, "EcoRisk/gi.png"),
  title = "global importance"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "ecosystem_balance",
  file = paste0(out_folder, "EcoRisk/eb.png"),
  title = "ecosystem balance")

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "carbon_stocks",
  file = paste0(out_folder, "EcoRisk/cs.png"),
  title = "carbon_stocks"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "carbon_fluxes",
  file = paste0(out_folder, "EcoRisk/cf.png"),
  title = "carbon_fluxes"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "water_stocks",
  file = paste0(out_folder, "EcoRisk/ws.png"),
  title = " water_stocks"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "water_fluxes",
  file = paste0(out_folder, "EcoRisk/wf.png"),
  title = " water_fluxes"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "nitrogen_stocks",
  file = paste0(out_folder, "EcoRisk/ns.png"),
  title = " nitrogen_stocks"
)

biospheremetrics::plot_ecorisk_map(
  ecorisk_object = ecorisk,
  plot_dimension = "nitrogen_fluxes",
  file = paste0(out_folder, "EcoRisk/nf.png"),
  title = " nitrogen_fluxes"
)

################# ecorisk biomes ################

biome_classes <- biospheremetrics::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 = paste0(lpj_input,"/GSWP3-W5E5/tas_gswp3-w5e5_1901-2016.clm"),
    elevation = paste0(lpj_input,"/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 <- biospheremetrics::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 = paste0(lpj_input,"/GSWP3-W5E5/tas_gswp3-w5e5_1901-2016.clm"),
    elevation = paste0(lpj_input,"/input_VERSION2/elevation.bin")
  ),
  time_span_reference = as.character(1901:1910), 
  savanna_proxy = list(pft_lai = 6),
  montane_arctic_proxy = list(elevation = 1000) 
)

biospheremetrics::plot_biomes_mercator(biome_ids = biome_classes$biome_id,
            file = paste0(out_folder,"EcoRisk/biomes_2005-2014.png"),
            order_legend = seq_len(19)
            )
biospheremetrics::plot_biomes_mercator(biome_ids = biome_classes_pi$biome_id,
            file = paste0(out_folder,"EcoRisk/biomes_1901-1910.png"),
            order_legend = seq_len(19)
            )

# compute median ecorisk values for biomes/large worldregions
ecorisk_disaggregated_full <- biospheremetrics::disaggregate_into_biomes(
  data = ecorisk,
  biome_class = biome_classes,
  type = "quantile",
  classes = "allbiomes"
)
ecorisk_disaggregated_full[is.na(ecorisk_disaggregated_full)] <- 0

ecorisk_disaggregated_4regions <- biospheremetrics::disaggregate_into_biomes(
  data = ecorisk,
  biome_class = biome_classes,
  type = "quantile",
  classes = "4biomes"
)

biospheremetrics::plot_ecorisk_radial_panel(
  data = ecorisk_disaggregated_full[-c(17,18,19),,], 
  biome_names = get_biome_names(1)[-c(17,18,19)],
  file = paste0(out_folder,"EcoRisk/EcoRisk_panel_1564_vs_2002.png"),
  use_quantile = TRUE,
  eps = TRUE
)

biospheremetrics::plot_ecorisk_radial_panel(
  data = ecorisk_disaggregated_4regions[,,], 
  biome_names = c("tropics","temperate","boreal","arctic"),
  file = paste0(out_folder,"EcoRisk/EcoRisk_4regions_1564_vs_2002.png"),
  use_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("/path/to/ecorisk_overtime_gamma.RData") # containing ecorisk object

ecorisk_overtime_allbiomes <- biospheremetrics::disaggregate_into_biomes(
  data = ecorisk,
  biome_class = biome_classes,
  type = "quantile",
  classes = "allbiomes"
)

biospheremetrics::plot_ecorisk_over_time_panel(
  data = ecorisk_overtime_allbiomes,
  timerange = c(1916,2003),
  biome_names = c("tropic","temperate","boreal","arctic"),
  file = paste0(out_folder,"overtime_panel.png"),
  eps = TRUE
)

ecorisk_overtime_biome16 <- biospheremetrics::disaggregate_into_biomes(
  data = ecorisk,
  biome_class = biome_classes,
  type = "quantile",
  classes = "allbiomes"
)

biospheremetrics::plot_ecorisk_over_time_panel(
  data = ecorisk_overtime_biome16[-c(3,17,18),,,],
  timerange = c(1916,2003),
  biome_names = get_biome_names(1)[-c(3,17,18)], 
  file = paste0(out_folder,"overtime_panel_16.png"),
  eps = TRUE
)