Added ozone yield shock functions

.buildlibrary

@@ -1,4 +1,4 @@
-ValidationKey: '12606930'
+ValidationKey: '12813440'
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'
 - 'Warning: namespace ''.*'' is not available and has been replaced'
@@ -8,3 +8,4 @@ AcceptedNotes:
 AutocreateReadme: yes
 allowLinterWarnings: no
 enforceVersionUpdate: no
+skipCoverage: no

CITATION.cff

@@ -2,8 +2,8 @@ cff-version: 1.2.0
 message: If you use this software, please cite it using the metadata from this file.
 type: software
 title: 'mrland: MadRaT land data package'
-version: 0.63.0
-date-released: '2024-10-15'
+version: 0.64.0
+date-released: '2024-10-25'
 abstract: The package provides land related data via the madrat framework.
 authors:
 - family-names: Dietrich
@@ -39,6 +39,8 @@ authors:
   given-names: David
 - family-names: Sauer
   given-names: Pascal
+- family-names: Tommey
+  given-names: Jake
 license: LGPL-3.0
 repository-code: https://github.com/pik-piam/mrland
 doi: 10.5281/zenodo.3822083

DESCRIPTION

@@ -1,8 +1,8 @@
 Type: Package
 Package: mrland
 Title: MadRaT land data package
-Version: 0.63.0
-Date: 2024-10-15
+Version: 0.64.0
+Date: 2024-10-25
 Authors@R: c(
     person("Jan Philipp", "Dietrich", , "[email protected]", role = c("aut", "cre")),
     person("Abhijeet", "Mishra", role = "aut"),
@@ -19,7 +19,8 @@ Authors@R: c(
     person("Stephen", "Wirth", role = "aut"),
     person("Felicitas", "Beier", role = "aut"),
     person("David", "Hoetten", role = "aut"),
-    person("Pascal", "Sauer", role = "aut")
+    person("Pascal", "Sauer", role = "aut"),
+    person("Jake", "Tommey", role = "aut")
   )
 Description: The package provides land related data via the madrat
     framework.

NAMESPACE

@@ -104,6 +104,7 @@ importFrom(raster,aggregate)
 importFrom(raster,brick)
 importFrom(raster,raster)
 importFrom(raster,rasterToPoints)
+importFrom(readxl,read_xlsx)
 importFrom(reshape2,acast)
 importFrom(stats,complete.cases)
 importFrom(stats,median)

R/calcOzoneYieldShock.R

@@ -0,0 +1,89 @@
+#' @title calcOzoneYieldShock
+#' @description calculate Ozone yield shocks
+#' Data from the EAT-Lancet deepdive on Ozone shock effects on crop yields.
+#' @param weighting     use of different weights (totalCrop (default),
+#'                      totalLUspecific, cropSpecific, crop+irrigSpecific,
+#'                      avlCropland, avlCropland+avlPasture)
+#' @param marginal_land  Defines which share of marginal land should be included (see options below) and
+#'                whether suitable land under irrigated conditions ("irrigated"), under rainfed conditions ("rainfed")
+#'                or suitability under rainfed conditions including currently irrigated land (rainfed_and_irrigated)
+#'                should be used. Options combined via ":"
+#'                The different marginal land options are:
+#' \itemize{
+#' \item \code{"all_marginal"}: All marginal land (suitability index between 0-0.33) is included as suitable
+#' \item \code{"q33_marginal"}: The bottom tertile (suitability index below 0.13) of the
+#' marginal land area is excluded.
+#' \item \code{"q50_marginal"}: The bottom  half (suitability index below 0.18) of the
+#' marginal land area is excluded.
+#' \item \code{"q66_marginal"}: The first and second tertile (suitability index below 0.23) of
+#' the marginal land area are excluded.
+#' \item \code{"q75_marginal"}: The first, second and third quartiles (suitability index below 0.25)
+#' of the marginal land are are excluded
+#' \item \code{"no_marginal"}: Areas with a suitability index of 0.33 and lower are excluded.
+#' \item \code{"magpie"}: Returns "all_marginal:rainfed_and_irrigated",
+#'                        "q33_marginal:rainfed_and_irrigated" and
+#'                        "no_marginal:rainfed_and_irrigated" in a magclass object to be used as magpie input.
+#' }
+#' @return magpie object in cellular resolution
+#' @author Jake Tommey
+#' @examples
+#' \dontrun{
+#' calcOutput("OzoneYieldShock")
+#' }
+#' @importFrom mstools toolGetMappingCoord2Country
+
+calcOzoneYieldShock <- function(
+  weighting = "totalCrop",
+  marginal_land = "magpie" #nolint
+) {
+  shocks <- readSource("OzoneYieldShock", convert = "onlycorrect")
+  shocks <- add_columns(shocks, addnm = "y2020", dim = 2, fill = 0)
+  # time interpolate the shock
+  shocks <- magclass::time_interpolate(
+    shocks,
+    interpolated_year = seq(1965, 2150, 5),
+    integrate_interpolated_years = TRUE,
+    extrapolation_type = "constant"
+  )
+
+  mapping <- toolGetMappingCoord2Country()
+  mapping$coordiso <- paste(mapping$coords, mapping$iso, sep = ".")
+  shocks <- toolAggregate(shocks, mapping, from = "iso", to = "coordiso")
+  getSets(shocks) <- c("x.y.iso", "year", "rcp.crop")
+
+  # create a MAgPIE object with the correct shape.
+  cropNames <- toolGetMapping("MAgPIE_LPJmL.csv", type = "sectoral", where = "mappingfolder")$MAgPIE
+  yieldShock <- new.magpie(
+    cells_and_regions = getCells(shocks),
+    years = getYears(shocks),
+    names = c(paste0("rcp2p6.", cropNames), paste0("rcp7p0.", cropNames))
+  )
+
+  yieldShock[, , "tece"] <- shocks[, , "Wheat"]
+  yieldShock[, , "soybean"] <- shocks[, , "soybean"]
+  yieldShock[, , "rice_pro"] <- shocks[, , "Rice"]
+  yieldShock[, , "sugr_beet"] <- shocks[, , "sugarbeet"]
+  yieldShock[, , "maiz"] <- shocks[, , "Maize"]
+  otherCrops <- yieldShock[, , c("tece", "soybean", "rice_pro", "sugr_beet", "maiz"), invert = TRUE]
+  yieldShock[, , getItems(otherCrops, dim = 3)] <- (dimSums(shocks, dim = 3) / length(getItems(shocks, dim = 3)))
+  yieldShock[, , "pasture"] <- 0
+
+  getSets(yieldShock) <- c("x", "y", "iso", "year", "rcp", "crop")
+
+  cropAreaWeights <- calcOutput(
+    "YieldsWeight",
+    weighting = weighting,
+    marginal_land = marginal_land,
+    aggregate = FALSE
+  )
+
+  return(
+    list(
+      x = yieldShock,
+      weight = cropAreaWeights,
+      unit = "%",
+      description = "percentage yield shock due to ozone",
+      isocountries = FALSE
+    )
+  )
+}

R/calcYields.R

@@ -203,82 +203,13 @@ calcYields <- function(source = c(lpjml = "ggcmi_phase3_nchecks_9ca735cb", isimi
 
   }
 
-  # Weight for spatial aggregation
-  if (weighting == "totalCrop") {
-
-    cropAreaWeight <- dimSums(calcOutput("Croparea", sectoral = "kcr", physical = TRUE, irrigation = FALSE,
-                                         cellular = TRUE, cells = cells, aggregate = FALSE,
-                                         years = "y1995", round = 6),
-                              dim = 3) + 10e-10
-
-  } else if (weighting %in% c("totalLUspecific", "cropSpecific", "crop+irrigSpecific")) {
-
-    crop <- calcOutput("Croparea", sectoral = "kcr", physical = TRUE, irrigation = TRUE,
-                       cellular = TRUE, cells = cells, aggregate = FALSE, years = "y1995", round = 6)
-
-    past <- calcOutput("LanduseInitialisation", aggregate = FALSE, cellular = TRUE, nclasses = "seven",
-                       input_magpie = TRUE, cells = cells, years = "y1995", round = 6)[, , "past"]
-
-    if (weighting == "crop+irrigSpecific") {
-
-      cropAreaWeight <- new.magpie(cells_and_regions = getCells(yields),
-                                   years = NULL,
-                                   names = getNames(yields),
-                                   fill = NA)
-      cropAreaWeight[, , findset("kcr")] <- crop + 10e-10
-      cropAreaWeight[, , "pasture"]      <- mbind(setNames(past + 10e-10, "irrigated"),
-                                                  setNames(past + 10e-10, "rainfed"))
-
-    } else if (weighting == "cropSpecific") {
-
-      cropAreaWeight <- new.magpie(cells_and_regions = getCells(yields),
-                                   years = NULL,
-                                   names = getNames(yields, dim = 1),
-                                   fill = NA)
-
-      cropAreaWeight[, , findset("kcr")] <- dimSums(crop, dim = 3.1) + 10e-10
-      cropAreaWeight[, , "pasture"]      <- past + 10e-10
-
-    } else {
-
-      cropAreaWeight <- new.magpie(cells_and_regions = getCells(yields),
-                                   years = NULL,
-                                   names = getNames(yields, dim = 1),
-                                   fill = (dimSums(crop, dim = 3) + 10e-10))
-
-      cropAreaWeight[, , "pasture"] <- past + 10e-10
-
-    }
-
-  } else if (weighting == "avlCropland") {
-
-    cropAreaWeight <- setNames(calcOutput("AvlCropland", marginal_land = marginal_land, cells = cells,
-                                          country_level = FALSE, aggregate = FALSE),
-                               NULL) + 10e-10
-
-  } else if (weighting == "avlCropland+avlPasture") {
-
-    avlCrop <- setNames(calcOutput("AvlCropland", marginal_land = marginal_land, cells = cells,
-                                   country_level = FALSE, aggregate = FALSE), "avlCrop")
-
-    lu1995  <- setYears(calcOutput("LanduseInitialisation", aggregate = FALSE, cellular = TRUE, nclasses = "seven",
-                                   input_magpie = TRUE, cells = cells, years = "y1995", round = 6), NULL)
-
-    cropAreaWeight <- new.magpie(cells_and_regions = getCells(yields),
-                                 years = NULL,
-                                 names = getNames(yields, dim = 1),
-                                 fill = avlCrop)
-
-    cropAreaWeight[, , "pasture"] <- pmax(avlCrop,
-                                          dimSums(lu1995[, , c("primforest", "secdforest", "forestry", "past")],
-                                                  dim = 3)) + 10e-10
-
-  } else {
-
-    stop("Weighting setting is not available.")
-  }
-
-  if (any(is.na(cropAreaWeight))) stop("NAs in weights.")
+  cropAreaWeight <- calcOutput(
+    "YieldsWeight",
+    cells = cells,
+    weighting = weighting,
+    marginal_land = marginal_land,
+    aggregate = FALSE
+  )
 
   # Special case for India case study
   if (indiaYields) {

R/calcYieldsWeight.R

@@ -0,0 +1,135 @@
+#' @title calcYieldsWeight
+#'
+#' @description This function calculates the crop area weightings to use for yields.
+#' @param weighting     use of different weights (totalCrop (default),
+#'                      totalLUspecific, cropSpecific, crop+irrigSpecific,
+#'                      avlCropland, avlCropland+avlPasture)
+#' @param cells         if cellular is TRUE: "magpiecell" for 59199 cells or "lpjcell" for 67420 cells
+#' @param marginal_land  Defines which share of marginal land should be included (see options below) and
+#'                whether suitable land under irrigated conditions ("irrigated"), under rainfed conditions ("rainfed")
+#'                or suitability under rainfed conditions including currently irrigated land (rainfed_and_irrigated)
+#'                should be used. Options combined via ":"
+#'                The different marginal land options are:
+#' \itemize{
+#' \item \code{"all_marginal"}: All marginal land (suitability index between 0-0.33) is included as suitable
+#' \item \code{"q33_marginal"}: The bottom tertile (suitability index below 0.13) of the
+#' marginal land area is excluded.
+#' \item \code{"q50_marginal"}: The bottom  half (suitability index below 0.18) of the
+#' marginal land area is excluded.
+#' \item \code{"q66_marginal"}: The first and second tertile (suitability index below 0.23) of
+#' the marginal land area are excluded.
+#' \item \code{"q75_marginal"}: The first, second and third quartiles (suitability index below 0.25)
+#' of the marginal land are are excluded
+#' \item \code{"no_marginal"}: Areas with a suitability index of 0.33 and lower are excluded.
+#' \item \code{"magpie"}: Returns "all_marginal:rainfed_and_irrigated",
+#'                        "q33_marginal:rainfed_and_irrigated" and
+#'                        "no_marginal:rainfed_and_irrigated" in a magclass object to be used as magpie input.
+#' }
+#' @return magpie object in cellular resolution
+#'
+#' @author Kristine Karstens, Felicitas Beier
+#'
+#' @examples
+#' \dontrun{
+#' calcOutput("YieldsWeight", yields, aggregate = FALSE)
+#' }
+#'
+#' @importFrom magpiesets findset
+#' @importFrom magclass getYears add_columns dimSums time_interpolate
+#' @importFrom madrat toolFillYears toolGetMapping toolTimeAverage
+#' @importFrom mstools toolHarmonize2Baseline
+#' @importFrom mrlandcore toolLPJmLVersion
+#' @importFrom stringr str_split
+#' @importFrom withr local_options
+
+calcYieldsWeight <- function(cells = "lpjcell", weighting = "totalCrop", marginal_land = "magpie") { # nolint
+
+  yieldNames <- toolGetMapping("MAgPIE_LPJmL.csv", type = "sectoral", where = "mappingfolder")$MAgPIE
+  isos <- toolGetMappingCoord2Country()
+  yieldCells <- paste(isos$coords, isos$iso, sep = ".")
+
+
+  # Weight for spatial aggregation
+  if (weighting == "totalCrop") {
+
+    cropAreaWeight <- dimSums(calcOutput("Croparea", sectoral = "kcr", physical = TRUE, irrigation = FALSE,
+                                         cellular = TRUE, cells = cells, aggregate = FALSE,
+                                         years = "y1995", round = 6),
+                              dim = 3) + 10e-10
+
+  } else if (weighting %in% c("totalLUspecific", "cropSpecific", "crop+irrigSpecific")) {
+
+    crop <- calcOutput("Croparea", sectoral = "kcr", physical = TRUE, irrigation = TRUE,
+                       cellular = TRUE, cells = cells, aggregate = FALSE, years = "y1995", round = 6)
+
+    past <- calcOutput("LanduseInitialisation", aggregate = FALSE, cellular = TRUE, nclasses = "seven",
+                       input_magpie = TRUE, cells = cells, years = "y1995", round = 6)[, , "past"]
+
+    if (weighting == "crop+irrigSpecific") {
+
+      cropAreaWeight <- new.magpie(cells_and_regions = yieldCells,
+                                   years = NULL,
+                                   names = c(paste(yieldNames, "irrigated", sep = "."),
+                                             paste(yieldNames, "rainfed", sep = ".")),
+                                   fill = NA)
+      cropAreaWeight[, , findset("kcr")] <- crop + 10e-10
+      cropAreaWeight[, , "pasture"]      <- mbind(setNames(past + 10e-10, "irrigated"),
+                                                  setNames(past + 10e-10, "rainfed"))
+
+    } else if (weighting == "cropSpecific") {
+
+      cropAreaWeight <- new.magpie(cells_and_regions = yieldCells,
+                                   years = NULL,
+                                   names = yieldNames,
+                                   fill = NA)
+
+      cropAreaWeight[, , findset("kcr")] <- dimSums(crop, dim = 3.1) + 10e-10
+      cropAreaWeight[, , "pasture"]      <- past + 10e-10
+
+    } else {
+
+      cropAreaWeight <- new.magpie(cells_and_regions = yieldCells,
+                                   years = NULL,
+                                   names = yieldNames,
+                                   fill = (dimSums(crop, dim = 3) + 10e-10))
+
+      cropAreaWeight[, , "pasture"] <- past + 10e-10
+
+    }
+
+  } else if (weighting == "avlCropland") {
+
+    cropAreaWeight <- setNames(calcOutput("AvlCropland", marginal_land = marginal_land, cells = cells,
+                                          country_level = FALSE, aggregate = FALSE),
+                               NULL) + 10e-10
+
+  } else if (weighting == "avlCropland+avlPasture") {
+
+    avlCrop <- setNames(calcOutput("AvlCropland", marginal_land = marginal_land, cells = cells,
+                                   country_level = FALSE, aggregate = FALSE), "avlCrop")
+
+    lu1995  <- setYears(calcOutput("LanduseInitialisation", aggregate = FALSE, cellular = TRUE, nclasses = "seven",
+                                   input_magpie = TRUE, cells = cells, years = "y1995", round = 6), NULL)
+
+    cropAreaWeight <- new.magpie(cells_and_regions = yieldCells,
+                                 years = NULL,
+                                 names = yieldNames,
+                                 fill = avlCrop)
+
+    cropAreaWeight[, , "pasture"] <- pmax(avlCrop,
+                                          dimSums(lu1995[, , c("primforest", "secdforest", "forestry", "past")],
+                                                  dim = 3)) + 10e-10
+
+  } else {
+
+    stop("Weighting setting is not available.")
+  }
+
+  if (any(is.na(cropAreaWeight))) stop("NAs in weights.")
+
+  return(list(x            = cropAreaWeight,
+              weight       = NULL,
+              unit         = "Mha",
+              description  = "Yields in tons per hectar for different crop types.",
+              isocountries = FALSE))
+}