Update land conversion cost calibration for cropland - FAO as target data set

CHANGELOG.md

@@ -12,6 +12,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **config.cfg** default for `cfg$gms$s29_treecover_max` changed from 0.4 to 1
 - **config.cfg** default for `cfg$gms$s29_fallow_max ` changed from 0.4 to 0
 - **config.cfg** default for `cfg$gms$s35_forest_damage ` changed from 2 to 0
+- **scripts** land conversion cost calibration for cropland - FAO as target data set instead of MAgPIEown
+- **default.cfg** settings for  land conversion cost calibration updated
 
 ### added
 - **scenario_config.csv** added column `NPI-revert`

config/default.cfg

@@ -26,7 +26,7 @@ cfg$input <- c(regional    = "rev4.116_h12_magpie.tgz",
                cellular    = "rev4.116_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz",
                validation  = "rev4.116_h12_validation.tgz",
                additional  = "additional_data_rev4.60.tgz",
-               calibration = "calibration_H12_27Sep24.tgz")
+               calibration = "calibration_H12_FAO30_03Feb25.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,20 +95,24 @@ cfg$recalibrate_landconversion_cost <- "ifneeded" #def "ifneeded"
 # Up to which accuracy shall be recalibrated?
 cfg$calib_accuracy_landconversion_cost <- 0.01         # def = 0.01
 # What is the maximum number of iterations if the precision goal is not reached?
-cfg$calib_maxiter_landconversion_cost <- 40           # def = 40
+cfg$calib_maxiter_landconversion_cost <- 20           # def = 20
 # Restart from existing calibration factors (TRUE or FALSE)
 cfg$restart_landconversion_cost <- FALSE           # def = FALSE
 # Number of lowpass filter iterations applied on calibration factors
 # for time steps 1995-2015
-cfg$lowpass_filter_landconversion_cost <- 1        # def= 1
+cfg$lowpass_filter_landconversion_cost <- 0        # def= 0
 # Set upper limit for cropland calibration factor
-cfg$cost_calib_max_landconversion_cost <- 3         # def= 3
+cfg$cost_calib_max_landconversion_cost <- 2.5         # def= 2.5
 # Set lower limit for cropland calibration factor
-cfg$cost_calib_min_landconversion_cost <- 0.05            # def= 0.05
+cfg$cost_calib_min_landconversion_cost <- 0.2            # def= 0.2
 # Selection type of calibration factors.
 # If FALSE, calibration factors from the last iteration are used.
 # If TRUE, calibration factors from the iteration with the lowest divergence are used.
 cfg$best_calib_landconversion_cost <- FALSE     # def = FALSE
+# Target data set that will be used for cropland calibration at regional level 
+# * "MAgPIEown": same data set as used for initialization of cropland
+# * "FAO": Data from FAOSTAT on cropland area
+cfg$cost_calib_hist_data <- "FAO"            # def = "FAO"
 
 # Settings for NPI/NDC recalculation
 # * (TRUE): NPI/NDC recalculation will be performed

config/projects/scenario_config_fsec.csv

@@ -79,6 +79,5 @@ gms$c73_build_demand;;;;;;;;;;;;;;;;;;;;;;;;50pc;;;;;;;;;;;;;;;;;;;;;;;;;;;
 input['cellular'];rev4.116_FSEC_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;rev4.116_FSEC_0bd54110_cellularmagpie_c200_MRI-ESM2-0-ssp119_lpjml-8e6c5eb1.tgz;rev4.116_FSEC_6819938d_cellularmagpie_c200_MRI-ESM2-0-ssp126_lpjml-8e6c5eb1.tgz;;rev4.116_FSEC_1b5c3817_cellularmagpie_c200_MRI-ESM2-0-ssp245_lpjml-8e6c5eb1.tgz;rev4.116_FSEC_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;rev4.116_FSEC_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz;rev4.116_FSEC_09a63995_cellularmagpie_c200_MRI-ESM2-0-ssp585_lpjml-8e6c5eb1.tgz;;;
 input['regional'];rev4.116_FSEC_magpie.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 input['validation'];rev4.116_FSEC_validation.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-input['additional'];additional_data_rev4.59.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-input['calibration'];calibration_FSEC_19Nov24.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+input['calibration'];calibration_FSEC_FAO_03Feb25.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 magicc_emis_scen;bjoernAR6_C_SSP2-NDC.mif;;;bjoernAR6_C_SSP2-PkBudg900.mif;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;bjoernAR6_C_SSP1-NDC.mif;;;;;;;;;;;;bjoernAR6_C_RemSDP-900-MagSSP1.mif;;

scripts/calibration/landconversion_cost.R

@@ -35,15 +35,27 @@ calibration_run <- function(putfolder, calib_magpie_name, logoption = 3, s_use_g
 
 # get ratio between modelled area and reference area
 
-getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpass_filter = 1) {
+getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpass_filter = 1, histData = "FAO", cost_min = 0.2) {
   require(magclass)
   require(magpie4)
   require(gdx2)
   y <- readGDX(gdx_file,"t")
   magpie <- land(gdx_file)[, y, "crop"]
-  hist <- dimSums(readGDX(gdx_file, "f10_land")[, , "crop"], dim = 1.2)
-  data <- hist[, y, "crop"]
-  shrLost <- (setYears(hist[getRegions(magpie), 2015, ], NULL) - setYears(hist[getRegions(magpie), 1995, ], NULL)) / setYears(hist[getRegions(magpie), 1995, ], NULL)
+  if (histData == "MAgPIEown") {
+    hist <- dimSums(readGDX(gdx_file, "f10_land")[, , "crop"], dim = 1.2)
+    data <- hist[, y, "crop"]
+  } else if (histData == "FAO") {
+    if(file.exists("calib_data.rds")) {
+      val <- readRDS("calib_data.rds")
+    } else {
+      val <- read.report("input/validation.mif", as.list = FALSE)
+      val <- val[getRegions(magpie),getYears(magpie),"historical.FAO_crop_past.Resources|Land Cover|+|Cropland (million ha)"]
+      names(dimnames(val)) <- names(dimnames(magpie))
+      getNames(val) <- "crop"
+      saveRDS(val, file = "calib_data.rds")
+    }
+    data <- val
+  }
   if (nregions(magpie) != nregions(data) | !all(getRegions(magpie) %in% getRegions(data))) {
     stop("Regions in MAgPIE do not agree with regions in reference calibration area data set!")
   }
@@ -52,18 +64,25 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
     out[out == 0] <- 1
     out[is.na(out)] <- 1
     getNames(out) <- NULL
-
-    out[which(out < 0, arr.ind = T)] <- 1
+    out[out < 0] <- 1
+    out <- lowpass(out,i = lowpass_filter)
   } else if (mode == "reward") {
+    shrLostHist <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
+    for (i in 2:length(y)) {
+      shrLostHist[ , y[i], ] <- (setYears(data[, y[i], ], NULL) - setYears(data[, y[i-1], ], NULL)) / setYears(data[, y[i-1], ], NULL)
+    }
+
     out <- magpie / data - 1
     out[is.na(out)] <- 0
+    out[is.infinite(out)] <- 0
     getNames(out) <- NULL
-
-    # only reward if share of cropland lost between 1995 and 2015 exceeds a certain threshold. Otherwise set to 0.
-    out[rownames(which(shrLost > -calib_accuracy, arr.ind = T)),,] <- 0
-    out[which(out < 0, arr.ind = T)] <- 0
+
+    # set reward to 0 if no cropland was lost in historic data set
+    out[shrLostHist < -calib_accuracy] <- 0
+    out[out < 0] <- 0
+    out <- lowpass(out,i = lowpass_filter)
+    out[,1,] <- 0
   }
-  out <- lowpass(out,i = lowpass_filter)
   return(magpiesort(out))
 }
 
@@ -84,14 +103,8 @@ time_series_reward <- function(calib_factor) {
   return(out2)
 }
 
-getHistCrop <- function() {
-  rep <- read.report("input/validation.mif", as.list = FALSE)
-  crop <- collapseNames(rep[, , "historical.FAO_crop_past.Resources|Land Cover|+|Cropland (million ha)"])
-  return(crop)
-}
-
 # Calculate the correction factor and save it
-update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, calib_file, cost_max = 3, cost_min = 0.05, calibration_step = "", n_maxcalib = 40, best_calib = FALSE) {
+update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, calib_file, cost_max = 2.5, cost_min = 0.2, calibration_step = "", n_maxcalib = 20, best_calib = FALSE, histData = "FAO") {
   require(magclass)
   require(magpie4)
   require(gdx2)
@@ -104,6 +117,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
 
   calib_correction_reward <- getCalibFactor(gdx_file, mode = "reward", calib_accuracy = calib_accuracy, lowpass_filter = lowpass_filter)
   calib_divergence_reward <- abs(calib_correction_reward)
+  calib_divergence_reward[calib_correction_reward == 0] <- 0
 
   ### -> in case it is the first step, it forces the initial factors to be equal to 1
   if (file.exists(calib_file)) {
@@ -124,29 +138,31 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
     # use stricter divergence threshold in first 5 calibration_step steps
     cost_acc_reached <- calib_divergence_cost <= ifelse(calibration_step < 6, 0.01, calib_accuracy)
     calib_factor_cost[cost_acc_reached] <- setNames(old_calib[, , "cost"], NULL)[cost_acc_reached]
-
+    calib_correction_cost[cost_acc_reached] <- 1
+
     reward_acc_reached <- calib_divergence_reward <= ifelse(calibration_step < 6, 0.01, calib_accuracy)
     calib_factor_reward[reward_acc_reached] <- setNames(old_calib[, , "reward"], NULL)[reward_acc_reached]
+    calib_correction_reward[reward_acc_reached] <- 0
   }
 
   if (!is.null(cost_max)) {
-    # if reward exists, set cost calibration to cost_max
-    reward_exists <- (calib_factor_reward > 0)
-    calib_factor_cost[reward_exists] <- cost_max
-
+
     # set value for India to cost_max for better convergence
     if ("IND" %in% getRegions(calib_factor_cost)) {
-      calib_factor_cost["IND", , ] <- cost_max
+      calib_factor_cost["IND", -1 , ] <- cost_max
     }
 
-    # set value for CAZ to cost_max for better convergence
-    if ("CAZ" %in% getRegions(calib_factor_cost)) {
-      calib_factor_cost["CAZ", , ] <- cost_max
+    # set value for CAZ to cost_max for better convergence; only needed for LUH2 due to mismatch between LUH2 and FAO historical data for cropland.
+    if (histData == "MAgPIEown") {
+      if ("CAZ" %in% getRegions(calib_factor_cost)) {
+        calib_factor_cost["CAZ", -1 , ] <- cost_max
+      }
     }
 
     above_limit <- (calib_factor_cost >= cost_max)
     calib_factor_cost[above_limit] <- cost_max
     calib_divergence_cost[getRegions(calib_factor_cost), , ][above_limit] <- 0
+
   }
 
   if (!is.null(cost_min)) {
@@ -201,6 +217,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
         " unit: -",
         paste0(" note: Best calibration factor from the run"),
         " origin: scripts/calibration/landconversion_cost.R (path relative to model main directory)",
+        paste(" Calibration settings:",  "calib_accuracy=", calib_accuracy, "lowpass_filter=", lowpass_filter, "cost_max=", cost_max, "cost_min=", cost_min, "n_maxcalib=",n_maxcalib, "best_calib=",best_calib, "histData=",histData),
         paste0(" creation date: ", date())
       )
       write.magpie(round(calib_best_full, 3), calib_file, comment = comment)
@@ -227,6 +244,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
       " unit: -",
       paste0(" note: Calibration step ", calibration_step),
       " origin: scripts/calibration/landconversion_cost.R (path relative to model main directory)",
+      paste(" Calibration settings:",  "calib_accuracy=", calib_accuracy, "lowpass_filter=", lowpass_filter, "cost_max=", cost_max, "cost_min=", cost_min, "n_maxcalib=",n_maxcalib, "best_calib=",best_calib, "histData=",histData),
       paste0(" creation date: ", date())
     )
 
@@ -238,19 +256,20 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
 
 
 
-calibrate_magpie <- function(n_maxcalib = 40,
-                             restart = TRUE,
-                             calib_accuracy = 0.05,
-                             cost_max = 3,
-                             cost_min = 0.05,
+calibrate_magpie <- function(n_maxcalib = 20,
+                             restart = FALSE,
+                             calib_accuracy = 0.01,
+                             cost_max = 2.5,
+                             cost_min = 0.2,
                              calib_magpie_name = "magpie_calib",
                              lowpass_filter = 1,
                              calib_file = "modules/39_landconversion/input/f39_calib.csv",
                              putfolder = "land_conversion_cost_calib_run",
                              data_workspace = NULL,
                              logoption = 3,
                              debug = FALSE,
-                             best_calib = FALSE) {
+                             best_calib = FALSE,
+                             histData = "FAO") {
   require(magclass)
 
   if (!restart) {
@@ -280,6 +299,7 @@ calibrate_magpie <- function(n_maxcalib = 40,
   # delete calib_magpie_gms in the main folder
   unlink(paste0(calib_magpie_name, ".*"))
   unlink("fulldata.gdx")
+  unlink("calib_data.rds")
 
   cat("\nLand conversion cost calibration finished\n")
 }

scripts/start/extra/recalibrateH16.R

@@ -19,7 +19,7 @@ source("scripts/start_functions.R")
 source("config/default.cfg")
 cfg$input['regional'] <- "rev4.116_36f73207_magpie.tgz"
 cfg$input['validation'] <- "rev4.116_36f73207_validation.tgz"
-cfg$input['calibration'] <- "calibration_H16_27Sep24.tgz"
+cfg$input['calibration'] <- "calibration_H16_FAO_03Feb25.tgz"
 cfg$input['cellular'] <- "rev4.116_36f73207_44a213b6_cellularmagpie_c400_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1_clusterweight-ba4466a8.tgz"
 
 cfg$results_folder <- "output/:title:"
@@ -30,4 +30,4 @@ cfg$output <- c("rds_report")
 cfg$force_replace <- TRUE
 cfg$qos <- "priority"
 start_run(cfg,codeCheck=FALSE)
-magpie4::submitCalibration("H16")
+magpie4::submitCalibration("H16_FAO")

scripts/start/extra/recalibrate_FSEC.R

@@ -25,9 +25,8 @@ cfg$results_folder                  <- "output/:title:"
 cfg$recalibrate                     <- TRUE # required when penality_apr22 activated
 cfg$best_calib                      <- TRUE
 cfg$recalibrate_landconversion_cost <- TRUE
-cfg$best_calib_landconversion_cost  <- TRUE
 cfg$output                          <- c("rds_report")
 cfg$force_replace                   <- TRUE
 cfg$qos <- "priority"
 start_run(cfg, codeCheck = FALSE)
-magpie4::submitCalibration("FSEC")
+magpie4::submitCalibration("FSEC_FAO")

scripts/start/extra/recalibrate_default.R

@@ -25,4 +25,4 @@ cfg$output <- c("rds_report")
 cfg$force_replace <- TRUE
 cfg$qos <- "priority"
 start_run(cfg,codeCheck=FALSE)
-magpie4::submitCalibration("H12")
+magpie4::submitCalibration("H12_FAO")

scripts/start/projects/project_ABCDR.R

@@ -33,7 +33,7 @@ cfg$qos <- "standby_dayMax"
 
 cfg$input['regional'] <- "rev4.116_36f73207_magpie.tgz"
 cfg$input['validation'] <- "rev4.116_36f73207_validation.tgz"
-cfg$input['calibration'] <- "calibration_H16_27Sep24.tgz"
+cfg$input['calibration'] <- "calibration_H16_FAO_03Feb25.tgz"
 cfg$input['cellular'] <- "rev4.116_36f73207_bd86374e_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1_clusterweight-ba4466a8.tgz"
 
 ssp <- "SSP2"

scripts/start_functions.R

@@ -476,7 +476,8 @@ start_run <- function(cfg, scenario = NULL, codeCheck = TRUE, lock_model = TRUE,
                      data_workspace = cfg$val_workspace,
                      logoption = 3,
                      debug = cfg$debug,
-                     best_calib = cfg$best_calib_landconversion_cost)
+                     best_calib = cfg$best_calib_landconversion_cost,
+                     histData = cfg$cost_calib_hist_data)
     cat("Land conversion cost calibration factor calculated!\n")
   }