From 147dcf2f0070e0610fd046557fdeb9520e39ea72 Mon Sep 17 00:00:00 2001 From: David Hoetten Date: Mon, 27 May 2024 14:07:31 +0200 Subject: [PATCH 1/7] added output report with corresponding r markdown --- scripts/output/projects/MAgPIE_EU_report.R | 3 + .../projects/MAgPIE_EU_report_notebook.Rmd | 249 ++++++++++++++++++ 2 files changed, 252 insertions(+) create mode 100644 scripts/output/projects/MAgPIE_EU_report.R create mode 100644 scripts/output/projects/MAgPIE_EU_report_notebook.Rmd diff --git a/scripts/output/projects/MAgPIE_EU_report.R b/scripts/output/projects/MAgPIE_EU_report.R new file mode 100644 index 000000000..dae269932 --- /dev/null +++ b/scripts/output/projects/MAgPIE_EU_report.R @@ -0,0 +1,3 @@ +readArgs("outputdir") +rmarkdown::render("./scripts/output/projects/MAgPIE_EU_report_notebook.Rmd", + output_dir = paste0(outputdir, "/validation")) diff --git a/scripts/output/projects/MAgPIE_EU_report_notebook.Rmd b/scripts/output/projects/MAgPIE_EU_report_notebook.Rmd new file mode 100644 index 000000000..3af76a13c --- /dev/null +++ b/scripts/output/projects/MAgPIE_EU_report_notebook.Rmd @@ -0,0 +1,249 @@ +--- +title: "Compare MAgPIE EU with LUH Data" +output: + html_document: + df_print: paged + code_folding: hide +--- + + +This output report validates spatial pattern of the EU regions against LUH2v2 data. + +# Setup and read data +```{r, setup, include=FALSE} +knitr::opts_knit$set(root.dir = dirname(dirname(dirname(getwd())))) +``` + +```{r load-libraries, echo = T, results = 'hide', message=FALSE, warning=FALSE} +# setup working dir and packages, read outputdir +readArgs("outputdir") +library(madrat) +library(mrcommons) +library(magpie4) + +library(ggplot2) +library(plotly) +library(gridExtra) +library(patchwork) + +library(dplyr) +library(tidyr) +library(stringr) +``` + +```{r} +print(paste0("Script started for output directory: ", outputdir)) +``` + +```{r read-data, echo = T, results = 'hide', message = FALSE, warning = FALSE} +compareYears <- c(1995, 2000, 2005, 2010) +# ----- Read model output and luh2v2 data, process and bind them to a validation object ----- +#### load and process reference data to match with model data +luh2v2 <- calcOutput("LanduseInitialisation", aggregate = FALSE, cellular = TRUE, + nclasses = "seven", cells = "lpjmlcell") +luh2v2 <- luh2v2[, compareYears, ] +# add subdimension for reference data +getNames(luh2v2) <- paste0(getNames(luh2v2), ".ref") + +#### load and process model data to match with reference data +path2landoutput <- file.path(outputdir, "cell.land_0.5.mz") +cellLand <- read.magpie(path2landoutput) +cellLand <- cellLand[, compareYears, ] +# add subdimension for model data +getNames(cellLand) <- paste0(getNames(cellLand), ".mod") +load(file.path(outputdir, "spatial_header.rda")) + +##### bind model and reference data +validationObj <- mbind(cellLand, luh2v2) + +##### extract EU countries +mapping <- toolGetMapping("regionmappingH12.csv") +countriesEU <- mapping$CountryCode[mapping$RegionCode == "EUR"] +``` + + +# Plot data + +```{r, load-plot-functions, echo = T, results = 'hide'} +# ----- Create a plotting functions ----- +convertToWideDataframe <- function(magclassObj) { + as.data.frame(magclassObj) %>% + pivot_wider(names_from = "Data2", values_from = "Value") +} + +subsetCountryCluster <- function(magclassObj, countries) { + magclassObj[intersect(countries, getItems(magclassObj, split=TRUE)[[1]][[1]]), , ] +} + +# combines a scatter plot with maps such that facets are aligned +plotCombined <- function(validationObj) { + # create the plots + p1 <- plotScatter(validationObj) + p2 <- plotMaps(validationObj) + + # combine and align the plots using patchwork + combinedPlot <- p1 / p2 + + # print the combined plot + print(combinedPlot) +} + +# create scatter plot +plotScatter <- function(validationObj) { + validationObj <- toolCoord2Isocell(validationObj, "lpjcell") + validationObj <- subsetCountryCluster(validationObj, countriesEU) + + ### Create quality measure dataframe with RMSE, MAE a.o. + qualityMeasuresDf <- NULL # list of dataframes, one for each year + numericYears <- as.numeric(gsub("y", "", getItems(validationObj, 2))) + for (year in numericYears) { + # get a vector with the quality measures + qualityMeasures <- luplot::qualityMeasure( + validationObj[, year, "mod"], validationObj[, year, "ref"], + measures = c("Willmott", "Willmott refined", "Nash Sutcliffe", "RMSE", "MAE") + ) + + # Create text for the year facet + text <- "" + for (i in seq_along(qualityMeasures)) { + text <- paste0(text, stringr::str_trunc(names(qualityMeasures)[i], 12, ellipsis = "."), + " : ", qualityMeasures[i], "\n") + } + + # Create df for that year + qualityMeasuresDf[[year]] <- data.frame( + Year = year, + text = text + ) + } + qualityMeasBinded <- dplyr::bind_rows(qualityMeasuresDf) # bind all year dataframes + + # Create the plot with facets + validationDF <- convertToWideDataframe(validationObj) + + p <- ggplot(validationDF, + aes(x = mod, y = ref, reg = Region, cell = Cell, relativeErr = (mod - ref) / ref)) + + geom_point(size = 0.5) + + geom_abline(color = "#663939", size = 1.5) + + facet_wrap(~ Year, scales = "free", ncol = 4) + # Create facets based on 'Year' + labs(x = "MAgPIE output", y = "luh2v2") + + theme(panel.background = element_rect(fill = "gray55"), + panel.grid.major = element_line(color = "gray62"), + panel.grid.minor = element_line(color = "grey58")) + + # Add quality measure text + p <- p + geom_text(x = 0, y = Inf, aes(label = text), color = "#ffb5b5", # bright color contrasting dark background + hjust = 0, vjust = 1.1, nudge_x = 10, size = 2.6, family = "sans", fontface = "bold", + data = qualityMeasBinded, inherit.aes = FALSE) + + return(p) +} + +# create maps +plotMaps <- function(validationObj) { + # create a ggplot object using luplot + relErr <- (validationObj[, , paste0("mod")] - validationObj[, , paste0("ref")]) / + (validationObj[, , paste0("ref")]) + + p <- luplot::plotmap2(relErr, + legend_range = c(-2, 2), legendname = "relative\n diff \n to \n LUH2v2", ncol = 4, + midcol = "#ffffff", lowcol = "blue", highcol = "red", midpoint = 0, + title = "" + ) + + # adjust the plot + p <- p + coord_sf(xlim = c(-10, 40), ylim = c(35, 70)) + facet_wrap(~ Year, ncol = 4) + + theme(aspect.ratio = 1, legend.title = element_text(size = 8)) + + p <- p + guides(fill = guide_colorbar(barheight = 5, barwidth = 0.2)) + return(p) +} + +# create and interactive scatterplot for eu countries of a year +plotInteractiveSub <- function(validationObj, year) { + p <- plotScatter(validationObj[, year, ]) + + p <- p + ggtitle(paste0("EU Countries ", year)) + p <- plotly::ggplotly(p) + p +} +``` + +## crop +```{r crop,message=FALSE, warning=FALSE, out.width="100%"} +type <- "crop" +plotCombined(validationObj[, , type]) +``` + +```{r crop2, message=FALSE, warning=FALSE, fig.align="center"} +plotInteractiveSub(validationObj[, , type], 2010) +``` + +## pasture +```{r past,message=FALSE, warning=FALSE, out.width="100%"} +type <- "past" +plotCombined(validationObj[, , type]) +``` + +```{r past2, message=FALSE, warning=FALSE, fig.align="center"} +plotInteractiveSub(validationObj[, , type], 2010) +``` + +## forestry +```{r forestry,message=FALSE, warning=FALSE, out.width="100%"} +type <- "forestry" +plotCombined(validationObj[, , type]) +``` + +```{r forestry2, message=FALSE, warning=FALSE, fig.align="center"} +plotInteractiveSub(validationObj[, , type], 2010) +``` + +## primforest +```{r primforest,message=FALSE, warning=FALSE, out.width="100%"} +type <- "primforest" +plotCombined(validationObj[, , type]) +``` + +```{r primforest2, message=FALSE, warning=FALSE, fig.align="center"} +plotInteractiveSub(validationObj[, , type], 2010) +``` + +## secdforest +```{r secdforest,message=FALSE, warning=FALSE, out.width="100%"} +type <- "secdforest" +plotCombined(validationObj[, , type]) +``` + +```{r secdforest2, message=FALSE, warning=FALSE, fig.align="center"} +plotInteractiveSub(validationObj[, , type], 2010) +``` + +## urban + +```{r urban,message=FALSE, warning=FALSE, out.width="100%"} +type <- "urban" +plotCombined(validationObj[, , type]) +``` + +```{r urban2, message=FALSE, warning=FALSE, fig.align="center"} +plotInteractiveSub(validationObj[, , type], 2010) +``` + +## other +```{r other, message=FALSE, warning=FALSE, out.width="100%"} +type <- "other" +plotCombined(validationObj[, , type]) +``` + +```{r other2, message=FALSE, warning=FALSE, fig.align="center"} +plotInteractiveSub(validationObj[, , type], 2010) +``` From eec7b46dd0e8e8c6bcd7386e7ea6310a8a748faf Mon Sep 17 00:00:00 2001 From: David Hoetten Date: Mon, 27 May 2024 14:57:33 +0200 Subject: [PATCH 2/7] renamed report --- scripts/output/projects/{MAgPIE_EU_report.R => EU_report.R} | 2 +- .../projects/{MAgPIE_EU_report_notebook.Rmd => EU_report.Rmd} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename scripts/output/projects/{MAgPIE_EU_report.R => EU_report.R} (53%) rename scripts/output/projects/{MAgPIE_EU_report_notebook.Rmd => EU_report.Rmd} (100%) diff --git a/scripts/output/projects/MAgPIE_EU_report.R b/scripts/output/projects/EU_report.R similarity index 53% rename from scripts/output/projects/MAgPIE_EU_report.R rename to scripts/output/projects/EU_report.R index dae269932..d8b2126a2 100644 --- a/scripts/output/projects/MAgPIE_EU_report.R +++ b/scripts/output/projects/EU_report.R @@ -1,3 +1,3 @@ readArgs("outputdir") -rmarkdown::render("./scripts/output/projects/MAgPIE_EU_report_notebook.Rmd", +rmarkdown::render("./scripts/output/projects/EU_report.Rmd", output_dir = paste0(outputdir, "/validation")) diff --git a/scripts/output/projects/MAgPIE_EU_report_notebook.Rmd b/scripts/output/projects/EU_report.Rmd similarity index 100% rename from scripts/output/projects/MAgPIE_EU_report_notebook.Rmd rename to scripts/output/projects/EU_report.Rmd From 8b0fe3bda3a0abb1f5550ee5e7a7bf6a73d09da8 Mon Sep 17 00:00:00 2001 From: David Hoetten Date: Mon, 27 May 2024 14:59:14 +0200 Subject: [PATCH 3/7] added changelog entry --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b558901db..6690fd6cd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). - **35_natveg** `vm_land(j2,"forestry")` included in NPI/NDC constraint `q35_min_forest` ### added +- **scripts** added Magpie EU report `EU_report.R` that uses `EU_report.Rmd` - **default.cfg** added cropland growth constraint `cfg$gms$s30_annual_max_growth` - **default.cfg** added technical cost for missing BII increase `cfg$gms$s44_cost_bii_missing` - **default.cfg** added settings for new price-driven bioenergy realization `1st2ndgen_priced_feb24`: `cfg$gms$s60_2ndgen_bioenergy_dem_min_post_fix`, `cfg$gms$c60_bioenergy_subsidy_fix_SSP2`, `s60_bioenergy_gj_price_1st`, From 48917e7c6319cbbe93465c6197648be6fb84dc38 Mon Sep 17 00:00:00 2001 From: DavidhoPIK <101278418+DavidhoPIK@users.noreply.github.com> Date: Mon, 27 May 2024 15:05:10 +0200 Subject: [PATCH 4/7] Update CHANGELOG.md --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6690fd6cd..457320063 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,7 +21,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/). - **35_natveg** `vm_land(j2,"forestry")` included in NPI/NDC constraint `q35_min_forest` ### added -- **scripts** added Magpie EU report `EU_report.R` that uses `EU_report.Rmd` +- **scripts** added output report `EU_report.R` that uses `EU_report.Rmd` - **default.cfg** added cropland growth constraint `cfg$gms$s30_annual_max_growth` - **default.cfg** added technical cost for missing BII increase `cfg$gms$s44_cost_bii_missing` - **default.cfg** added settings for new price-driven bioenergy realization `1st2ndgen_priced_feb24`: `cfg$gms$s60_2ndgen_bioenergy_dem_min_post_fix`, `cfg$gms$c60_bioenergy_subsidy_fix_SSP2`, `s60_bioenergy_gj_price_1st`, From 94ab66ac512ac0ce100db996e5f9413033ccb4b1 Mon Sep 17 00:00:00 2001 From: DavidhoPIK <101278418+DavidhoPIK@users.noreply.github.com> Date: Tue, 28 May 2024 10:42:45 +0200 Subject: [PATCH 5/7] Minor update of comment in EU_report.Rmd --- scripts/output/projects/EU_report.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/output/projects/EU_report.Rmd b/scripts/output/projects/EU_report.Rmd index 3af76a13c..88afd2d1a 100644 --- a/scripts/output/projects/EU_report.Rmd +++ b/scripts/output/projects/EU_report.Rmd @@ -6,7 +6,7 @@ output: code_folding: hide --- + This output report validates spatial pattern of the EU regions against LUH2v2 data. -# Setup and read data +# Setup data and plot functions ```{r, setup, include=FALSE} knitr::opts_knit$set(root.dir = dirname(dirname(dirname(getwd())))) ``` @@ -45,15 +47,15 @@ print(paste0("Script started for output directory: ", outputdir)) ```{r read-data, echo = T, results = 'hide', message = FALSE, warning = FALSE} compareYears <- c(1995, 2000, 2005, 2010) -# ----- Read model output and luh2v2 data, process and bind them to a validation object ----- -#### load and process reference data to match with model data +# ----- Read and process model output and luh2v2 data; then bind data to a validation object ----- +# load and process reference data to match with model data luh2v2 <- calcOutput("LanduseInitialisation", aggregate = FALSE, cellular = TRUE, nclasses = "seven", cells = "lpjmlcell") luh2v2 <- luh2v2[, compareYears, ] # add subdimension for reference data getNames(luh2v2) <- paste0(getNames(luh2v2), ".ref") -#### load and process model data to match with reference data +# load and process model data to match with reference data path2landoutput <- file.path(outputdir, "cell.land_0.5.mz") cellLand <- read.magpie(path2landoutput) cellLand <- cellLand[, compareYears, ] @@ -61,29 +63,28 @@ cellLand <- cellLand[, compareYears, ] getNames(cellLand) <- paste0(getNames(cellLand), ".mod") load(file.path(outputdir, "spatial_header.rda")) -##### bind model and reference data +# bind model and reference data validationObj <- mbind(cellLand, luh2v2) -##### extract EU countries +# get EU countries mapping <- toolGetMapping("regionmappingH12.csv") countriesEU <- mapping$CountryCode[mapping$RegionCode == "EUR"] ``` - -# Plot data - ```{r, load-plot-functions, echo = T, results = 'hide'} # ----- Create a plotting functions ----- +# helper function to convert magclass obj to wide dataframe convertToWideDataframe <- function(magclassObj) { as.data.frame(magclassObj) %>% pivot_wider(names_from = "Data2", values_from = "Value") } +# helper function to subset all cluster of a country subsetCountryCluster <- function(magclassObj, countries) { magclassObj[intersect(countries, getItems(magclassObj, split=TRUE)[[1]][[1]]), , ] } -# combines a scatter plot with maps such that facets are aligned +# function combines a scatter plot with maps such that facets are aligned plotCombined <- function(validationObj) { # create the plots p1 <- plotScatter(validationObj) @@ -96,12 +97,12 @@ plotCombined <- function(validationObj) { print(combinedPlot) } -# create scatter plot +# create scatter plot for eu countries faceted by year plotScatter <- function(validationObj) { validationObj <- toolCoord2Isocell(validationObj, "lpjcell") validationObj <- subsetCountryCluster(validationObj, countriesEU) - ### Create quality measure dataframe with RMSE, MAE a.o. + # create quality measure dataframe with RMSE, MAE a.o. qualityMeasuresDf <- NULL # list of dataframes, one for each year numericYears <- as.numeric(gsub("y", "", getItems(validationObj, 2))) for (year in numericYears) { @@ -111,14 +112,14 @@ plotScatter <- function(validationObj) { measures = c("Willmott", "Willmott refined", "Nash Sutcliffe", "RMSE", "MAE") ) - # Create text for the year facet + # create text for the year facet text <- "" for (i in seq_along(qualityMeasures)) { text <- paste0(text, stringr::str_trunc(names(qualityMeasures)[i], 12, ellipsis = "."), " : ", qualityMeasures[i], "\n") } - # Create df for that year + # create dataframe for that year qualityMeasuresDf[[year]] <- data.frame( Year = year, text = text @@ -126,7 +127,7 @@ plotScatter <- function(validationObj) { } qualityMeasBinded <- dplyr::bind_rows(qualityMeasuresDf) # bind all year dataframes - # Create the plot with facets + # create the plot with facets validationDF <- convertToWideDataframe(validationObj) p <- ggplot(validationDF, @@ -139,7 +140,7 @@ plotScatter <- function(validationObj) { panel.grid.major = element_line(color = "gray62"), panel.grid.minor = element_line(color = "grey58")) - # Add quality measure text + # add quality measure text p <- p + geom_text(x = 0, y = Inf, aes(label = text), color = "#ffb5b5", # bright color contrasting dark background hjust = 0, vjust = 1.1, nudge_x = 10, size = 2.6, family = "sans", fontface = "bold", data = qualityMeasBinded, inherit.aes = FALSE) @@ -147,7 +148,7 @@ plotScatter <- function(validationObj) { return(p) } -# create maps +# create EU maps faceted by year plotMaps <- function(validationObj) { # create a ggplot object using luplot relErr <- (validationObj[, , paste0("mod")] - validationObj[, , paste0("ref")]) / @@ -167,16 +168,17 @@ plotMaps <- function(validationObj) { return(p) } -# create and interactive scatterplot for eu countries of a year +# create an interactive scatterplot for eu countries of a year plotInteractiveSub <- function(validationObj, year) { p <- plotScatter(validationObj[, year, ]) - p <- p + ggtitle(paste0("EU Countries ", year)) p <- plotly::ggplotly(p) - p + return(p) } ``` +# Plot data + ## crop ```{r crop,message=FALSE, warning=FALSE, out.width="100%"} type <- "crop" From 773ace41ed54694d5c6e22a6d557a2463bcf9e5f Mon Sep 17 00:00:00 2001 From: David Hoetten Date: Fri, 21 Jun 2024 10:50:40 +0200 Subject: [PATCH 7/7] removed creation of validation subfolder --- scripts/output/projects/EU_report.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/output/projects/EU_report.R b/scripts/output/projects/EU_report.R index d8b2126a2..0eec8d917 100644 --- a/scripts/output/projects/EU_report.R +++ b/scripts/output/projects/EU_report.R @@ -1,3 +1,3 @@ readArgs("outputdir") rmarkdown::render("./scripts/output/projects/EU_report.Rmd", - output_dir = paste0(outputdir, "/validation")) + output_dir = outputdir)