-
Notifications
You must be signed in to change notification settings - Fork 17
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #49 from jdrtommey/ozone_final
Added ozone yield shock functions
- Loading branch information
Showing
14 changed files
with
443 additions
and
86 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 | ||
| ) | ||
| ) | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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)) | ||
| } |
Oops, something went wrong.