diff --git a/CHANGELOG.md b/CHANGELOG.md index a854989ca..31e584de0 100644 --- a/CHANGELOG.md +++ b/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` diff --git a/DESCRIPTION b/DESCRIPTION index f927eaf72..7e2b5847f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -14,7 +14,7 @@ Imports: gms (>= 0.24.0), here, iamc, - lucode2 (>= 0.36.0), + lucode2 (>= 0.51.1), luplot, luscale (>= 2.27.9), lusweave, diff --git a/README.md b/README.md index 55d364db8..0c3ec04f1 100644 --- a/README.md +++ b/README.md @@ -33,7 +33,7 @@ model source code via the R package goxygen package and run the main function (goxygen) in the main folder of the model. The resulting documentation can be found in the folder "doc". -Please find a set of tutorials here https://magpiemodel.github.io/tutorials/. +Please find a set of tutorials here https://magpiemodel.github.io/tutorials. This guide will give you a brief technical introduction in how to install, run and use the model and how to analyse the model output. diff --git a/config/default.cfg b/config/default.cfg index e0987238b..3c255d640 100644 --- a/config/default.cfg +++ b/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 diff --git a/config/projects/scenario_config_fsec.csv b/config/projects/scenario_config_fsec.csv index c830f7fed..636aa5fb0 100644 --- a/config/projects/scenario_config_fsec.csv +++ b/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;; diff --git a/main.gms b/main.gms index 6f9a944e8..cc705b7bc 100644 --- a/main.gms +++ b/main.gms @@ -160,13 +160,13 @@ $title magpie * md5sum: NA * Repository: https://rse.pik-potsdam.de/data/magpie/public * -* Used data set: additional_data_rev4.59.tgz +* Used data set: additional_data_rev4.60.tgz * md5sum: NA * Repository: https://rse.pik-potsdam.de/data/magpie/public * -* Used data set: calibration_H12_27Sep24.tgz -* md5sum: NA -* Repository: https://rse.pik-potsdam.de/data/magpie/public +* Used data set: calibration_H12_FAO30_03Feb25.tgz +* md5sum: aba0b877f383fefabc558d79180fc43f +* Repository: /Users/flo/Development/input_data/ * * Low resolution: c200 * High resolution: 0.5 @@ -179,11 +179,11 @@ $title magpie * * Regionscode: 62eff8f7 * -* Regions data revision: 4.114 +* Regions data revision: 4.116 * * lpj2magpie settings: * * LPJmL data: MRI-ESM2-0:ssp370 -* * Revision: 4.114 +* * Revision: 4.116 * * aggregation settings: * * Input resolution: 0.5 @@ -195,7 +195,7 @@ $title magpie * * Call: withCallingHandlers(expr, message = messageHandler, warning = warningHandler, error = errorHandler) * * -* Last modification (input data): Sun Oct 27 00:37:36 2024 +* Last modification (input data): Thu Feb 6 14:08:31 2025 * *###################### R SECTION END (VERSION INFO) ########################### diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R index 12ab1bb0e..bf3356d86 100644 --- a/scripts/calibration/landconversion_cost.R +++ b/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,11 +256,11 @@ 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", @@ -250,7 +268,8 @@ calibrate_magpie <- function(n_maxcalib = 40, 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") } diff --git a/scripts/start/extra/recalibrateH16.R b/scripts/start/extra/recalibrateH16.R index 4cc1c33ff..b4c1071b6 100644 --- a/scripts/start/extra/recalibrateH16.R +++ b/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") diff --git a/scripts/start/extra/recalibrate_FSEC.R b/scripts/start/extra/recalibrate_FSEC.R index f651a5f03..b4fd187ed 100644 --- a/scripts/start/extra/recalibrate_FSEC.R +++ b/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") diff --git a/scripts/start/extra/recalibrate_default.R b/scripts/start/extra/recalibrate_default.R index 78de9855e..b792deb4e 100644 --- a/scripts/start/extra/recalibrate_default.R +++ b/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") diff --git a/scripts/start/projects/project_ABCDR.R b/scripts/start/projects/project_ABCDR.R index c37eb4194..dfbd722ca 100644 --- a/scripts/start/projects/project_ABCDR.R +++ b/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" diff --git a/scripts/start_functions.R b/scripts/start_functions.R index 24ae800f6..cfc129755 100644 --- a/scripts/start_functions.R +++ b/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") }