Merge pull request #589 from pfuehrlich-pik/master

magpie release 4.6.10
magpiemodel · Aug 16, 2023 · a9c1230 · a9c1230
2 parents 09a26cf + 49c8b3a
commit a9c1230
Show file tree

Hide file tree

Showing 12 changed files with 220 additions and 131 deletions.
diff --git a/.github/workflows/test-code.yaml b/.github/workflows/test-code.yaml
@@ -18,14 +18,15 @@ jobs:
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
 
       - name: Remove .Rprofile
         run: rm .Rprofile
 
       - uses: r-lib/actions/setup-r@v2
         with:
-          extra-repositories: https://cloud.r-project.org https://rse.pik-potsdam.de/r/packages/
+          use-public-rspm: true
+          extra-repositories: "https://rse.pik-potsdam.de/r/packages"
 
       - name: pre-commit Hook
         run: source(".githooks/pre-commit")

diff --git a/.zenodo.json b/.zenodo.json
@@ -1,6 +1,6 @@
 {
   "title": "MAgPIE - An Open Source land-use modeling framework",
-  "version": "4.6.9",
+  "version": "4.6.10",
   "creators": [
     {
       "name": "Dietrich, Jan Philipp",
@@ -109,5 +109,5 @@
   "license": {
     "id": "AGPL-3.0-or-later"
   },
-  "publication_date": "2023-07-27"
+  "publication_date": "2023-08-16"
 }
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,6 +5,20 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 
+## [4.6.10] - 2023-08-16
+
+### changed
+- **config**  update preprocessing to newest input data v4.88 with new transport costs
+- **GitHub action** the github action is now faster, because it installs binary packages from Posit Package Manager
+- **scenario_config.csv** update preprocessing to newest input data v4.88
+- **scripts** output.R is now faster, because it no longer searches runfolder renvs for full.gms files
+
+### fixed
+- **scripts** check_config does not warn about c_input_gdx_path anymore while running empty model
+- **scripts** fixed erronoeous if clause in output.R
+- **scripts** fixed output/extra/disaggregation_LUH2.R. The script was not working any more because magpie4::protectedArea was changed to return protected area for all land types. Moreover, the script now also works for runs without dynamic forestry (default run) but with a warning message.
+
+
 ## [4.6.9] - 2023-07-27
 
 ### fixed
@@ -737,7 +751,8 @@ This release version is focussed on consistency between the MAgPIE setup and the
 First open source release of the framework. See [MAgPIE 4.0 paper](https://doi.org/10.5194/gmd-12-1299-2019) for more information.
 
 
-[Unreleased]: https://github.com/magpiemodel/magpie/compare/v4.6.9...develop
+[Unreleased]: https://github.com/magpiemodel/magpie/compare/v4.6.10...develop
+[4.6.10]: https://github.com/magpiemodel/magpie/compare/v4.6.9...v4.6.10
 [4.6.9]: https://github.com/magpiemodel/magpie/compare/v4.6.8...v4.6.9
 [4.6.8]: https://github.com/magpiemodel/magpie/compare/v4.6.7...v4.6.8
 [4.6.7]: https://github.com/magpiemodel/magpie/compare/v4.6.6...v4.6.7

diff --git a/CITATION.cff b/CITATION.cff
@@ -125,8 +125,8 @@ authors:
     email: [email protected]
 
 title: MAgPIE - An Open Source land-use modeling framework
-version: 4.6.9
-date-released: 2023-07-27
+version: 4.6.10
+date-released: 2023-08-16
 repository-code: https://github.com/magpiemodel/magpie
 keywords:
   - landuse

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -17,7 +17,7 @@ Imports:
     iamc,
     lucode2 (>= 0.36.0),
     luplot,
-    luscale (>= 2.27.4),
+    luscale (>= 2.27.9),
     lusweave,
     m4fsdp,
     madrat,

diff --git a/README.md b/README.md
@@ -24,8 +24,8 @@ https://www.pik-potsdam.de/research/projects/activities/land-use-modelling/magpi
 A framework description paper has been published in
 Geoscientific Model Development (GMD): https://doi.org/10.5194/gmd-12-1299-2019
 
-The model documentation for version 4.6.9 can be found at
-https://rse.pik-potsdam.de/doc/magpie/4.6.9/
+The model documentation for version 4.6.10 can be found at
+https://rse.pik-potsdam.de/doc/magpie/4.6.10/
 
 A most recent version of the documentation can also be extracted from the
 model source code via the R package goxygen
@@ -188,7 +188,7 @@ [email protected]
 Please contact [email protected]
 
 ## CITATION
-See file CITATION.cff or the [How-to-Cite section](https://rse.pik-potsdam.de/doc/magpie/4.6.9/#how-to-cite) in the model documentation for information how to cite the model.
+See file CITATION.cff or the [How-to-Cite section](https://rse.pik-potsdam.de/doc/magpie/4.6.10/#how-to-cite) in the model documentation for information how to cite the model.
 
 ## AUTHORS
 See list of authors in CITATION.cff

diff --git a/config/default.cfg b/config/default.cfg
@@ -22,11 +22,11 @@ cfg$model <- "main.gms"   #def = "main.gms"
 #### input settings ####
 
 # which input data sets should be used?
-cfg$input <- c(regional    = "rev4.87_h12_magpie.tgz",
-               cellular    = "rev4.87_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz",
-               validation  = "rev4.87_h12_validation.tgz",
+cfg$input <- c(regional    = "rev4.88_h12_magpie.tgz",
+               cellular    = "rev4.88_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz",
+               validation  = "rev4.88_h12_validation.tgz",
                additional  = "additional_data_rev4.43.tgz",
-               calibration = "calibration_H12_per_ton_fao_may22_glo_14Jul23.tgz")
+               calibration = "calibration_H12_per_ton_fao_may22_glo_08Aug23.tgz")
 
 # NOTE: It is recommended to recalibrate the model when changing cellular input data
 #       as well as for any other setting that would affect initial values in the model,
@@ -95,7 +95,7 @@ cfg$restart_landconversion_cost <- FALSE           # def = FALSE
 # calibration factor in each calibration iteration (0-1)
 cfg$damping_factor_landconversion_cost <- 0.96        # def= 0.96
 # Set upper limit for cropland calibration factor
-cfg$cost_calib_max_landconversion_cost <- 3            # def= 3
+cfg$cost_calib_max_landconversion_cost <- 3         # def= 3
 # Set lower limit for cropland calibration factor
 cfg$cost_calib_min_landconversion_cost <- 0.05            # def= 0.05
 # Selection type of calibration factors.

diff --git a/config/scenario_config.csv b/config/scenario_config.csv
@@ -79,7 +79,7 @@ gms$s73_timber_demand_switch;;;;;;;;;;;;;;;;;;;;;;1;1;1;1;1;;;;;1;1;0;;;;;;;;;;;
 gms$s35_forest_damage;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;4;4;4;4;4;;;;;;;
 gms$c32_shock_scenario;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;none;002lin2030;004lin2030;008lin2030;016lin2030;;;;;;;
 gms$c35_shock_scenario;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;none;002lin2030;004lin2030;008lin2030;016lin2030;;;;;;;
-input['cellular'];;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;rev4.87_h12_0bd54110_cellularmagpie_c200_MRI-ESM2-0-ssp119_lpjml-8e6c5eb1.tgz;rev4.87_h12_6819938d_cellularmagpie_c200_MRI-ESM2-0-ssp126_lpjml-8e6c5eb1.tgz;rev4.87_h12_1b5c3817_cellularmagpie_c200_MRI-ESM2-0-ssp245_lpjml-8e6c5eb1.tgz;rev4.87_h12_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;rev4.87_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz;rev4.87_h12_09a63995_cellularmagpie_c200_MRI-ESM2-0-ssp585_lpjml-8e6c5eb1.tgz;
+input['cellular'];;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;rev4.88_h12_0bd54110_cellularmagpie_c200_MRI-ESM2-0-ssp119_lpjml-8e6c5eb1.tgz;rev4.88_h12_6819938d_cellularmagpie_c200_MRI-ESM2-0-ssp126_lpjml-8e6c5eb1.tgz;rev4.88_h12_1b5c3817_cellularmagpie_c200_MRI-ESM2-0-ssp245_lpjml-8e6c5eb1.tgz;rev4.88_h12_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;rev4.88_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz;rev4.88_h12_09a63995_cellularmagpie_c200_MRI-ESM2-0-ssp585_lpjml-8e6c5eb1.tgz;
 gms$c52_land_carbon_sink_rcp;;nocc;nocc_hist;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;RCP19;RCP26;RCP45;RCP60;RCPBU;RCPBU;
 gms$c57_macc_version;;;;;;;;;;;;;;;;;;;;;;;;;;;PBL_2022;PBL_2022;PBL_2022;PBL_2022;;;;;;;;;;;;;;;
 gms$c57_macc_scenario;;;;;;;;;;;;;;;;;;;;;;;;;;;Default;Optimistic;Default;Optimistic;;;;;;;;;;;;;;;

diff --git a/output.R b/output.R
@@ -35,14 +35,16 @@ library(gms)
 runOutputs <- function(comp=NULL, output=NULL, outputdir=NULL, submit=NULL) {
   choose_folder <- function(title="Please choose a folder") {
     # try to use find because it is significantly quicker than list.dirs
-    tmp <- try(system("find ./output -name 'full.gms'", intern=TRUE,  ignore.stderr = TRUE), silent=TRUE)
+    tmp <- try(system("find ./output -path './output/*/renv' -prune -o -name 'full.gms'",
+                      intern=TRUE,  ignore.stderr = TRUE), silent=TRUE)
     if("try-error" %in% class(tmp) || length(tmp)==0) {
       tmp <- base::list.dirs("./output/",recursive=TRUE)
       dirs <- NULL
       for (i in seq_along(tmp)) {
         if (file.exists(file.path(tmp[i],"full.gms"))) dirs <- c(dirs,sub("./output/","",tmp[i]))
       }
     } else {
+      tmp <- grep("/full.gms$", tmp, value = TRUE)
       dirs <- sub("full.gms","",sub("./output/","",tmp, fixed=TRUE), fixed=TRUE)
     }
     dirs <- sort(dirs)
@@ -62,7 +64,7 @@ runOutputs <- function(comp=NULL, output=NULL, outputdir=NULL, submit=NULL) {
     }
     identifier <- tmp
     # PATTERN
-    if(length(identifier==1) && identifier==(length(dirs)+1)){
+    if(length(identifier) == 1 && identifier == length(dirs) + 1) {
       cat("\nInsert the search pattern or the regular expression: ")
       pattern <- gms::getLine()
       id <- grep(pattern=pattern, dirs[-1], perl=TRUE)

diff --git a/scripts/output/extra/disaggregation.R b/scripts/output/extra/disaggregation.R
@@ -28,12 +28,15 @@ map_file <- Sys.glob(file.path(outputdir, "clustermap_*.rds"))
 gdx <- file.path(outputdir, "fulldata.gdx")
 land_hr_file <- file.path(outputdir, "avl_land_full_t_0.5.mz")
 urban_land_hr_file <- file.path(outputdir, "f34_urbanland_0.5.mz")
+wdpa_hr_file <- file.path(outputdir, "wdpa_baseline_0.5.mz")
+consv_prio_hr_file <- file.path(outputdir, "consv_prio_areas_0.5.mz")
+land_consv_hr_out_file <- file.path(outputdir, "cell.conservation_land_0.5.mz")
 land_hr_out_file <- file.path(outputdir, "cell.land_0.5.mz")
 land_hr_share_out_file <- file.path(outputdir, "cell.land_0.5_share.mz")
 croparea_hr_share_out_file <- file.path(outputdir, "cell.croparea_0.5_share.mz")
 land_hr_split_file <- file.path(outputdir, "cell.land_split_0.5.mz")
 land_hr_shr_split_file <- file.path(outputdir, "cell.land_split_0.5_share.mz")
-luh_side_layers        <- file.path(outputdir, "luh2_side_layers_0.5.mz")
+luh_side_layers <- file.path(outputdir, "luh2_side_layers_0.5.mz")
 bii_hr_out_file <- file.path(outputdir, "cell.bii_0.5.mz")
 
 cfg <- gms::loadConfig(file.path(outputdir, "config.yml"))
@@ -56,7 +59,7 @@ if (length(map_file) > 1) {
   write.magpie(x, sub(".mz", ".nc", file), comment = comment, verbose = FALSE)
 }
 
-.dissagcrop <- function(gdx, land_hr, map) {
+.dissagcrop <- function(gdx, land_hr, map_file) {
   area <- croparea(gdx,
     level = "cell", products = "kcr",
     product_aggr = FALSE, water_aggr = FALSE
@@ -69,11 +72,81 @@ if (length(map_file) > 1) {
   crop_shr <- setNames(crop_shr[, getYears(area_shr), "crop"], NULL)
 
   # calculate crop area as share of total cell area
-  area_shr_hr <- madrat::toolAggregate(area_shr, map, to = "cell") * crop_shr
+  area_shr_hr <- madrat::toolAggregate(area_shr, map_file, to = "cell") * crop_shr
   return(area_shr_hr)
 }
 
-.dissagBII <- function(gdx, map, dir) {
+.dissagLandConsv <- function(gdx, cfg, map_file, wdpa_hr_file, consv_prio_hr_file) {
+  land_consv_lr <- readGDX(gdx, "pm_land_conservation", react = "silent")
+  wdpa_hr <- read.magpie(wdpa_hr_file)
+  map <- readRDS(map_file)
+
+  # create full time series
+  land_consv_hr <- new.magpie(map[, "cell"], getYears(land_consv_lr), getNames(wdpa_hr))
+  land_consv_hr[, getYears(land_consv_hr), ] <- wdpa_hr[, nyears(wdpa_hr), ]
+  land_consv_hr[, getYears(wdpa_hr), ] <- wdpa_hr
+
+  if (!all(c(cfg$gms$c22_protect_scenario, cfg$gms$c22_protect_scenario_noselect) %in% "none")) {
+    if (file.exists(consv_prio_hr_file)) {
+      consv_prio_all <- read.magpie(consv_prio_hr_file)
+      consv_prio_hr <- new.magpie(
+        cells_and_regions = map[, "cell"],
+        names = getNames(consv_prio_all, dim = 2), fill = 0
+      )
+      iso <- readGDX(gdx, "iso")
+      consv_iso <- readGDX(gdx, "policy_countries22")
+      consv_iso <- consv_iso[consv_iso %in% getItems(consv_prio_all, dim = 1.1)]
+      consv_select <- cfg$gms$c22_protect_scenario
+      consv_noselect <- cfg$gms$c22_protect_scenario_noselect
+
+      if (consv_noselect != "none") {
+        consv_prio_hr <- collapseDim(consv_prio_all[, , consv_noselect], dim = 3.1)
+      }
+      if (consv_select != "none") {
+        consv_prio_hr[consv_iso, , ] <- collapseDim(consv_prio_all[consv_iso, , consv_select], dim = 3.1)
+      } else if (consv_select == "none") {
+        consv_prio_hr[consv_iso, , ] <- 0
+      }
+      # future conservation only pertains to natveg
+      consv_prio_hr[, , c("crop", "past", "forestry", "urban")] <- 0
+      consv_fader <- readGDX(gdx, "p22_conservation_fader", format = "first_found")
+      consv_prio_hr <- consv_prio_hr * consv_fader[, getYears(land_consv_hr), ]
+
+      # add conservation priority areas
+      land_consv_hr <- (land_consv_hr + consv_prio_hr)
+    } else {
+      warning(paste(
+        "Future land conservation used in MAgPIE run but high resolution",
+        "conservation priority data for disaggregation not found."
+      ))
+    }
+  }
+  # Due to internal constraints and compensation (e.g. NDC forest conservation)
+  # the actual land conservation can sometimes be smaller than the land
+  # conservation in the input data (this can especially happen also if
+  # land restoration is switched off). Therefore a scaling is applied here separately
+  # for grassland and natural vegetation
+  natveg <- c("primforest", "secdforest", "other")
+  consv_sum_lr <- mbind(
+    dimSums(land_consv_lr[, , "past"], 3.2),
+    setNames(dimSums(land_consv_lr[, , natveg], dim = 3), "natveg")
+  )
+  consv_sum_hr_agg <- mbind(
+    toolAggregate(land_consv_hr[, , "past"], map, from = "cell", to = "cluster"),
+    toolAggregate(setNames(dimSums(land_consv_hr[, , natveg], dim = 3), "natveg"),
+      map,
+      from = "cell", to = "cluster"
+    )
+  )
+  consv_scaling <- consv_sum_lr / consv_sum_hr_agg
+  consv_scaling[is.na(consv_scaling) | is.infinite(consv_scaling)] <- 1
+  consv_scaling <- toolAggregate(consv_scaling, map, from = "cluster", to = "cell")
+  land_consv_hr[, , "past"] <- consv_scaling[, , "past"] * land_consv_hr[, , "past"]
+  land_consv_hr[, , natveg] <- consv_scaling[, , "natveg"] * land_consv_hr[, , natveg]
+  return(land_consv_hr)
+}
+
+.dissagBII <- function(gdx, map_file, dir) {
   # Biodiversity intactness indicator (BII) at cluster level
   bii_lr <- BII(gdx,
     file = NULL, level = "cell", mode = "auto", landClass = "all",
@@ -119,7 +192,9 @@ if (length(map_file) > 1) {
 # Prepare data for disaggregation
 # ========================================
 
+# ----------------------
 # Load input data
+# ----------------------
 land_ini_lr <- readGDX(gdx, "f10_land", "f_land", format = "first_found")[, "y1995", ]
 land_lr <- land(gdx, sum = FALSE, level = "cell")
 land_ini_hr <- read.magpie(land_hr_file)[, "y1995", ]
@@ -152,7 +227,9 @@ if (any(land_ini_hr < 0)) {
   land_ini_hr[which(land_ini_hr < 0, arr.ind = T)] <- 0
 }
 
-### Read in hr urban land
+# -----------------------------
+# Read in hr urban land
+# -----------------------------
 if (cfg$gms$urban == "exo_nov21") {
   urban_land_hr <- read.magpie(urban_land_hr_file)
   ssp <- cfg$gms$c34_urban_scenario
@@ -162,37 +239,24 @@ if (cfg$gms$urban == "exo_nov21") {
   urban_land_hr <- "static"
 }
 
-### Get land conservation data
-wdpa_hr <- NULL
-if (file.exists(file.path(outputdir, "wdpa_baseline_0.5.mz"))) {
-  wdpa_hr <- file.path(outputdir, "wdpa_baseline_0.5.mz")
-}
-consv_prio_hr <- NULL
-if (!all(c(cfg$gms$c22_protect_scenario, cfg$gms$c22_protect_scenario_noselect) %in% "none")) {
-  if (file.exists(file.path(outputdir, "consv_prio_areas_0.5.mz"))) {
-    consv_prio_all <- read.magpie(file.path(outputdir, "consv_prio_areas_0.5.mz"))
-    consv_prio_hr <- new.magpie(
-      cells = getCells(consv_prio_all),
-      names = getNames(consv_prio_all, dim = 2), fill = 0
-    )
-    iso <- readGDX(gdx, "iso")
-    consv_iso <- readGDX(gdx, "policy_countries22")
-    consv_iso <- consv_iso[consv_iso %in% getItems(consv_prio_all, dim = 1.1)]
-    consv_select <- cfg$gms$c22_protect_scenario
-    consv_noselect <- cfg$gms$c22_protect_scenario_noselect
-
-    if (consv_noselect != "none") {
-      consv_prio_hr <- collapseDim(consv_prio_all[, , consv_noselect], dim = 3.1)
-    }
-    if (consv_select != "none") {
-      consv_prio_hr[consv_iso, , ] <- collapseDim(consv_prio_all[consv_iso, , consv_select], dim = 3.1)
-    } else if (consv_select == "none") {
-      consv_prio_hr[consv_iso, , ] <- 0
-    }
-  }
+# ----------------------------------------
+# Prepare land conservation data
+# ----------------------------------------
+
+land_consv_hr <- NULL
+if (file.exists(wdpa_hr_file)) {
+  land_consv_hr <- .dissagLandConsv(gdx, cfg, map_file, wdpa_hr_file, consv_prio_hr_file)
+
+  # Write gridded conservation land
+  .writeDisagg(land_consv_hr, land_consv_hr_out_file,
+    comment = "unit: Mha per grid-cell",
+    message = "Write outputs cell.conservation_land"
+  )
 }
 
-### Account for country-specific SNV shares in post-processing
+# -------------------------------------------------------------
+# Account for country-specific SNV shares in post-processing
+# -------------------------------------------------------------
 iso <- readGDX(gdx, "iso")
 snv_pol_iso <- readGDX(gdx, "policy_countries30")
 snv_pol_select <- readGDX(gdx, "s30_snv_shr", "s30_set_aside_shr")
@@ -227,14 +291,13 @@ land_hr <- interpolateAvlCroplandWeighted(
   x = land_lr,
   x_ini_lr = land_ini_lr,
   x_ini_hr = land_ini_hr,
-  avl_cropland_hr = avl_cropland_hr,
   map = map_file,
+  avl_cropland_hr = avl_cropland_hr,
   marginal_land = marginal_land,
-  snv_pol_shr = snv_pol_shr,
-  snv_pol_fader = snv_pol_fader,
   urban_land_hr = urban_land_hr,
-  wdpa_hr = wdpa_hr,
-  consv_prio_hr = consv_prio_hr
+  land_consv_hr = land_consv_hr,
+  snv_pol_shr = snv_pol_shr,
+  snv_pol_fader = snv_pol_fader
 )
 
 # Write output
@@ -375,7 +438,7 @@ if (grepl("grass", cfg$gms$past)) {
 }
 
 # Sort and rename
-land_ini_hr <- land_ini_hr[,,getNames(land_ini_lr)]
+land_ini_hr <- land_ini_hr[, , getNames(land_ini_lr)]
 getSets(land_ini_hr)["d3.1"] <- "land"
 
 # Disaggregate BII values to high resolution
@@ -386,14 +449,13 @@ land_bii_hr <- interpolateAvlCroplandWeighted(
   x = land_lr,
   x_ini_lr = land_ini_lr,
   x_ini_hr = land_ini_hr,
-  avl_cropland_hr = avl_cropland_hr,
   map = map_file,
+  avl_cropland_hr = avl_cropland_hr,
   marginal_land = marginal_land,
+  urban_land_hr = urban_land_hr,
+  land_consv_hr = land_consv_hr,
   snv_pol_shr = snv_pol_shr,
   snv_pol_fader = snv_pol_fader,
-  urban_land_hr = urban_land_hr,
-  wdpa_hr = wdpa_hr,
-  consv_prio_hr = consv_prio_hr,
   unit = "share"
 )