Plastics production data from PlasticsEurope, TextileExchange and IRSG

NAMESPACE

@@ -1,6 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
 export(MplotMulti)
+export(calcPlGaoCabrera2025)
+export(calcPlProduction)
+export(calcPlSectorPolymerSplit)
+export(calcPlSyntheticFibre)
+export(calcPlSyntheticRubber)
 export(getSources_mrmfa)
 import(GDPuc)
 import(madrat)
@@ -27,16 +32,22 @@ importFrom(dplyr,ungroup)
 importFrom(magclass,"getComment<-")
 importFrom(magclass,"getItems<-")
 importFrom(magclass,"getNames<-")
+importFrom(magclass,"getSets<-")
 importFrom(magclass,add_dimension)
 importFrom(magclass,as.magpie)
+importFrom(magclass,collapseDim)
 importFrom(magclass,collapseNames)
 importFrom(magclass,dimSums)
+importFrom(magclass,getComment)
 importFrom(magclass,getItems)
 importFrom(magclass,getNames)
 importFrom(magclass,getYears)
+importFrom(magclass,magpiesort)
 importFrom(magclass,mbind)
 importFrom(magclass,mselect)
+importFrom(magclass,new.magpie)
 importFrom(magclass,setNames)
+importFrom(magclass,setYears)
 importFrom(magclass,where)
 importFrom(magpiesets,findset)
 importFrom(magrittr,"%>%")

R/calcCeTrade.R

@@ -69,6 +69,7 @@ calcCeTrade <- function(subtype, category, HS = "92", target_years = NULL) {
   x <- as.magpie(df, temporal = 1, spatial = 2, datacol = 4)
   x <- toolCountryFill(x, fill = NA, verbosity = 2)
   x <- replace_non_finite(x, replace = 0)
+  x <- collapseDim(x)
 
   # ----------------------------------------------------------------------------
   # Aggregation to avoid intra-regional trade

R/calcPlBACI.R

@@ -64,45 +64,40 @@ calcPlBACI <- function(subtype, category, HS) {
       rename("sector" = "Target")
 
     # map UNEP-NGP polymers to polymers used in REMIND-MFA
-    if (category == "Waste"){
+    if (category == "Waste") {
       polymer_map <- toolGetMapping("polymermappingUNEP_NGP_waste.csv", type = "sectoral", where = "mrmfa")
     } else {
       polymer_map <- toolGetMapping("polymermappingUNEP_NGP.csv", type = "sectoral", where = "mrmfa")
     }
-    # use polymer use by sector as weights (summarize over all Regions,
-    # as polymer share by sector is constant over all Regions in OECD data);
+    # use polymer use by sector as weights (summarize over all Regions and Years);
     # use total polymer use over all sectors as weights for "General" sector
-    use <- calcOutput("PlOECD", subtype = "Use_2019_region", aggregate = TRUE) %>% as.data.frame()
-    use_by_sector <- use %>%
-      rename(sector = "Data2", polymer = "Data1") %>%
-      filter(.data$polymer != "Total") %>%
-      group_by(.data$sector, .data$polymer) %>%
-      summarize(value = sum(.data$Value)) %>%
-      mutate(sector = case_when(.data$sector == "Total" ~ "General", .default = .data$sector))
-    split <- merge(polymer_map, use_by_sector, by.y = "polymer", by.x = "Target") %>%
+    use_sector <- dimSums(calcOutput("PlGaoCabrera2025", aggregate = TRUE), dim = c("region", "year"))
+    use_total <- dimSums(use_sector, dim = "sector") %>% addDim(dim = 3.3, dimName = "sector", item = "General")
+    use_by_sector <- mbind(use_sector, use_total) %>% as.data.frame(rev = 3)
+    split <- merge(polymer_map, use_by_sector, by.y = "polymer", by.x = "Target_polymer") %>%
       group_by(.data$sector, .data$Source) %>%
       mutate(
-        total = sum(.data$value, na.rm = TRUE),
-        weight = .data$value / .data$total
+        total = sum(.data$.value, na.rm = TRUE),
+        weight = .data$.value / .data$total
       ) %>%
-      select(-"total", -"value")
+      select(-"total", -".value")
     df <- left_join(df, split, by = c("polymer" = "Source", "sector")) %>%
       mutate(value = .data$value * .data$weight)
-    # some weights are NaN for polymers that are not used in a specific sector according to OECD,
+    # some weights are NaN for polymers that are not used in a specific sector according to Gao & Cabrera,
     # throw a warning if these combination appear in the dataset
     nan <- df %>% filter(is.na(.data$weight))
     if (nrow(nan) > 0) {
       warning(paste(
         "The following sector-polymer combinations cannot be mapped from the BACI data
-        as they do not exist in the OECD dataset used for weighting:\n",
+        as they do not exist in the GaoCabrera dataset used for weighting:\n",
         paste(utils::capture.output(print(nan)), collapse = "\n")
       ))
     }
 
     df <- df %>%
       select(-"polymer", -"weight") %>%
-      rename("polymer" = "Target") %>%
-      group_by(.data$t, .data$exporter, .data$importer, .data$polymer, .data$sector) %>%
+      rename("polymer" = "Target_polymer") %>%
+      group_by(.data$t, .data$exporter, .data$importer, .data$type, .data$polymer, .data$sector) %>%
       summarize(value = sum(.data$value)) %>%
       ungroup()
   }
@@ -114,6 +109,9 @@ calcPlBACI <- function(subtype, category, HS) {
   df$importer[df$importer == "SCG"] <- "SRB"
   df$exporter[df$exporter == "ANT"] <- "CUW"
   df$importer[df$importer == "ANT"] <- "CUW"
+  # ZA1 is the Southern African Customs Union, assigned to South Africa
+  df$exporter[df$exporter == "ZA1"] <- "ZAF"
+  df$importer[df$importer == "ZA1"] <- "ZAF"
 
   # this mapgie object should actually contain two spatial dimensions "importer" and
   # "exporter", but since madrat does not support regional aggregation for two spatial

R/calcPlEoL_shares.R

@@ -52,7 +52,7 @@ calcPlEoL_shares <- function(subtype) {
   # - US data (1960-2018, fill gaps by linear interpolation)
   # - calculate shares
   # ---------------------------------------------------------------------------
-  eu <- readSource("PlasticsEurope", subtype = "PlasticEoL_EU", convert = FALSE) %>%
+  eu <- readSource("PlasticsEurope", subtype = "EoL_EU", convert = FALSE) %>%
     as.data.frame()
   cn_raw <- readSource("China_PlasticEoL", convert = FALSE)
   cn <- time_interpolate(cn_raw, interpolated_year = 1978:2021, integrate_interpolated_years = TRUE) %>%
@@ -125,15 +125,24 @@ calcPlEoL_shares <- function(subtype) {
     1 - (x_backcast[, 1950:2000, "Recycled"] + x_backcast[, 1950:2000, "Incinerated"])
   )
 
+  # ---------------------------------------------------------------------------
+  # Forecast EoL rates to 2022 (keep constant)
+  # ---------------------------------------------------------------------------
+  x_forecast <- time_interpolate(x_backcast,
+    interpolated_year = 2019:2022,
+    integrate_interpolated_years = TRUE,
+    extrapolation_type = "constant"
+  )
+
   # ---------------------------------------------------------------------------
   # Select data based on subtype
   # ---------------------------------------------------------------------------
   x <- switch(subtype,
-    "Recycled" = mselect(x_backcast, Data1 = "Recycled"),
-    "Landfilled" = mselect(x_backcast, Data1 = "Landfilled"),
-    "Incinerated" = mselect(x_backcast, Data1 = "Incinerated"),
-    "Collected" = mselect(x_backcast, Data1 = "Collected"),
-    "All" = x_backcast,
+    "Recycled" = mselect(x_forecast, Data1 = "Recycled"),
+    "Landfilled" = mselect(x_forecast, Data1 = "Landfilled"),
+    "Incinerated" = mselect(x_forecast, Data1 = "Incinerated"),
+    "Collected" = mselect(x_forecast, Data1 = "Collected"),
+    "All" = x_forecast,
     stop("Unsupported subtype: ", subtype)
   )
   if (subtype != "All") {
@@ -151,6 +160,7 @@ calcPlEoL_shares <- function(subtype) {
   weight[, , ] <- 0
   weight[, seq(2000, 2019, 1)] <- plOECD[, seq(2000, 2019, ), "Total"]
   weight[, seq(1950, 1999, 1)] <- plOECD[, 2000, "Total"]
+  weight[, seq(2020, 2022, 1)] <- plOECD[, 2019, "Total"]
 
   # ---------------------------------------------------------------------------
   # Return results

R/calcPlGaoCabrera.R

@@ -12,7 +12,7 @@
 #' @seealso \code{\link{readGaoCabrera2025}}, \code{\link{calcPlPottinger}}
 #' @importFrom magclass dimSums getNames<-
 calcPlGaoCabrera <- function() {
-  x <- readSource("GaoCabrera2025")
+  x <- readSource("GaoCabrera2025", subtype = "consumption")
 
   # sum all 14 polymers into a single total and convert kt -> Mt
   x <- dimSums(x, dim = 3) / 1000

R/calcPlGaoCabrera2025.R

@@ -0,0 +1,78 @@
+#' Calculate absolute consumption per polymer and sector from Gao & Cabrera-Serrenho (2025)
+#'
+#' @description
+#' Combines the apparent polymer consumption of Gao & Cabrera-Serrenho (2025)
+#' (\code{readSource("GaoCabrera2025", "consumption")})
+#' with their polymer-specific end-use sector distribution
+#' (\code{readSource("GaoCabrera2025", "sector_shares")}) to obtain, for every
+#' country and year, the consumption of each polymer in each of the 8 end-use
+#' sectors.
+#'
+#' @return A list in \code{\link[madrat]{calcOutput}} format with the absolute
+#'   consumption by country, year, type, polymer and sector (1978-2021).
+#'   Years 2020-2021 rely on incomplete Gao trade data.
+#' @author Leonie Schweiger
+#' @seealso \code{\link{readGaoCabrera2025}}, \code{\link{calcPlProduction}}
+#' @examples
+#' \dontrun{
+#' a <- calcOutput("PlGaoCabrera2025")
+#' }
+#' @importFrom magclass collapseDim
+#' @export
+calcPlGaoCabrera2025 <- function() {
+  # ---------------------------------------------------------------------------
+  # Inputs: apparent polymer consumption (ISO3 x year x polymer, kt) and the
+  # polymer-specific end-use sector distribution (sector x polymer, unitless).
+  # ---------------------------------------------------------------------------
+  cons <- readSource("GaoCabrera2025", subtype = "consumption")
+  # sector_shares carries both sector and polymer as name subdims (sector.polymer)
+  # with only dummy spatial/temporal dims; collapse those so it broadcasts cleanly.
+  shares <- collapseDim(readSource("GaoCabrera2025", subtype = "sector_shares", convert = FALSE))
+
+  # ---------------------------------------------------------------------------
+  # Outer product cons x shares -> absolute consumption per polymer and sector.
+  # `cons` has name dim polymer; `shares` has name dims sector.polymer. The
+  # product matches on the shared polymer subdim and broadcasts cons over sector
+  # (and the single shares value over region/year).
+  # ---------------------------------------------------------------------------
+  absolute <- cons * shares / 1000 # (region, year, polymer.sector), kt -> Mt
+
+  # Merge LLDPE into LDPE (report both polyethylenes together as LDPE)
+  polyMap <- data.frame(from = getItems(absolute, dim = "polymer"))
+  polyMap$to <- ifelse(polyMap$from == "LLDPE", "LDPE", polyMap$from)
+  absolute <- toolAggregate(absolute, rel = polyMap, dim = 3.1, from = "from", to = "to")
+  # TODO Other fibre category only includes acrylic fibres; however, for remind-mfa,
+  # this includes also PP fibre and Elastane (PU fibre) which need to be subtracted from PP and PUR
+
+  # ---------------------------------------------------------------------------
+  # Backcast data to 1950 using Geyer et al. (2023) global polymer consumption data
+  # ---------------------------------------------------------------------------
+  geyer <- readSource("Geyer", subtype = "Prod_1950-2015", convert = FALSE)
+  absoluteBackcasted <- toolBackcastByReference(absolute, geyer)
+
+  # ---------------------------------------------------------------------------
+  # Add the type subdim (Fibre / Rubber / Plastics) by splitting on the polymer
+  # groups.
+  # ---------------------------------------------------------------------------
+  fibrePolymers <- c("PET fibre", "Polyamide fibre", "Other fibre")
+  rubberPolymers <- "Rubbers"
+  plasticsPolymers <- setdiff(getItems(absoluteBackcasted, dim = "polymer"), c(fibrePolymers, rubberPolymers))
+
+  tagType <- function(polys, typeName) {
+    add_dimension(mselect(absoluteBackcasted, polymer = polys), dim = 3.1, add = "type", nm = typeName)
+  }
+  absTyped <- mbind(
+    tagType(plasticsPolymers, "Plastics"),
+    tagType(fibrePolymers, "Fibre"),
+    tagType(rubberPolymers, "Rubber")
+  )
+
+  return(list(
+    x = absTyped,
+    weight = NULL,
+    unit = "Mt",
+    description = paste(
+      "Apparent consumption per polymer and end-use sector, from Gao & Cabrera-Serrenho (2025)."
+    )
+  ))
+}

R/calcPlHVCinput.R

@@ -58,26 +58,16 @@ calcPlHVCinput <- function(subtype) {
 
   # merge with HVC input and calculate total HVC input by element
   HVC_input <- merge(HVC_input, carbon_content, by.x = "input", by.y = "Data1") %>%
-    mutate(input_by_element = .data$tonnes * .data$Value, hvc_type = case_when(.data$input == "butadiene" ~ "C4", .default = "HVC w/o C4")) %>%
+    mutate(
+      input_by_element = .data$tonnes * .data$Value,
+      hvc_type = case_when(.data$input == "butadiene" ~ "C4", .default = "HVC w/o C4")) %>%
     group_by(.data$polymer, .data$element, .data$hvc_type) %>%
     summarize(HVC_input = sum(.data$input_by_element))
   HVC_input_magclass <- as.magpie(HVC_input)
 
   # ---------------------------------------------------------------------------
-  # map polymers to REMIND-MFA polymers using production weights for ABS&SAN mapping to ABS_ASA_SAN
-  # ---------------------------------------------------------------------------
-  polymer_map <- toolGetMapping("polymermappingLeviCullen.csv", type = "sectoral", where = "mrmfa")
-  production_volumes <- readSource("LeviCullen", subtype = "Production", convert = FALSE)
-  weights <- mselect(production_volumes, stage = "total production", type = "plastics", to = polymerization$polymer, collapseNames = TRUE) %>% collapseDim()
-  HVC_input_mapped <- HVC_input_magclass %>%
-    toolAggregate(rel = polymer_map, dim = 3.1, from = "Source", to = "Target", weight = weights, wdim = 3.1)
-  polymerization_mapped <- polymerization_magclass %>%
-    toolAggregate(rel = polymer_map, dim = 3.1, from = "Source", to = "Target", weight = weights, wdim = 3.1)
-
-  # ---------------------------------------------------------------------------
-  # fill data for Fibres, Bioplastics and Other with assumptions:
-  # for Fibres, assume same input as for PET, since 85% of global Fibre production is polyester
-  # for Bioplastics and Other, regress HVC inputs by carbon content and assume typical polymerization yield
+  # fill data for Other thermosets with assumption:
+  # regress HVC inputs by carbon content and assume typical polymerization yield
   # according to Levi and Cullen 2018 (total input = 1.02)
   # ---------------------------------------------------------------------------
 
@@ -86,27 +76,20 @@ calcPlHVCinput <- function(subtype) {
     mselect(Element = "C") %>%
     collapseDim() # keep only C content, result: (polymer)
 
-  # --- Fibres: copy HVC input and polymerization yield from PET ---
-  fibres_HVC <- HVC_input_mapped[, , "PET"]
-  getItems(fibres_HVC, dim = 3.1) <- "Fibres"
-  fibres_polym <- polymerization_mapped[, , "PET"]
-  getItems(fibres_polym, dim = 3.1) <- "Fibres"
-
-  # --- Bioplastics & Others: linear regression of HVC input on carbon content ---
   # build a dataframe of known polymer carbon contents and their HVC inputs per (element, HVC type)
-  known_polymers <- getItems(HVC_input_mapped, dim = 3.1)
+  known_polymers <- getItems(HVC_input_magclass, dim = 3.1)
   cc_known <- as.data.frame(mselect(carbon_content_pl, Polymer = known_polymers)) %>%
     select("Data1", "Value") %>%
     rename(polymer = "Data1", carbon = "Value")
 
-  hvc_known <- as.data.frame(HVC_input_mapped) %>%
+  hvc_known <- as.data.frame(HVC_input_magclass) %>%
     select(-"Cell", -"Region", -"Year") %>%
     rename(polymer = "Data1", element = "Data2", hvc_type = "Data3", hvc_input = "Value")
 
   reg_data <- merge(hvc_known, cc_known, by = "polymer")
 
-  # for each element, fit lm on HVC w/o C4 only and predict; C4 input is assumed zero
-  missing_polymers <- c("BioPlastics", "Others")
+  # fit lm on HVC w/o C4 only and predict; C4 input is assumed zero
+  missing_polymers <- c("Other thermosets")
   cc_missing <- as.data.frame(mselect(carbon_content_pl, Polymer = missing_polymers)) %>%
     select("Data1", "Value") %>%
     rename(polymer = "Data1", carbon = "Value")
@@ -116,8 +99,8 @@ calcPlHVCinput <- function(subtype) {
     group_by(.data$element) %>%
     tidyr::nest() %>%
     mutate(
-      model = purrr::map(.data$data, ~ lm(hvc_input ~ carbon, data = .x)),
-      pred  = purrr::map(.data$model, ~ predict(.x, newdata = cc_missing))
+      model = purrr::map(.data$data, ~ stats::lm(hvc_input ~ carbon, data = .x)),
+      pred  = purrr::map(.data$model, function(m) stats::predict(m, newdata = cc_missing))
     ) %>%
     tidyr::unnest("pred") %>%
     mutate(
@@ -130,16 +113,28 @@ calcPlHVCinput <- function(subtype) {
   missing_HVC <- add_columns(missing_HVC, addnm = "C4", dim = "hvc_type", fill = 0)
   getItems(missing_HVC, dim = 1) <- NULL
 
-  # polymerization yield of 98% for Bioplastics and Others
+  # polymerization yield of 98%
   missing_polym <- new.magpie(
     names = missing_polymers,
     sets = c("region", "year", "polymer"), fill = 1.02
   )
   getItems(missing_polym, dim = 1) <- NULL
 
+  # ---------------------------------------------------------------------------
+  # map polymers to REMIND-MFA polymers using production weights
+  # ---------------------------------------------------------------------------
+  polymer_map <- toolGetMapping("polymermappingLeviCullen.csv", type = "sectoral", where = "mrmfa")
+  production_volumes <- readSource("LeviCullen", subtype = "Production", convert = FALSE)
+  weights <- mselect(production_volumes, stage = "total production", type = "plastics", to = polymerization$polymer,
+                     collapseNames = TRUE) %>% collapseDim()
+  HVC_input_mapped <- HVC_input_magclass %>%
+    toolAggregate(rel = polymer_map, dim = 3.1, from = "Source", to = "Target", weight = weights, wdim = 3.1)
+  polymerization_mapped <- polymerization_magclass %>%
+    toolAggregate(rel = polymer_map, dim = 3.1, from = "Source", to = "Target", weight = weights, wdim = 3.1)
+
   # --- combine all polymers ---
-  HVC_input_full <- mbind(HVC_input_mapped, fibres_HVC, missing_HVC)
-  polymerization_full <- mbind(polymerization_mapped, fibres_polym, missing_polym)
+  HVC_input_full <- mbind(HVC_input_mapped, missing_HVC)
+  polymerization_full <- mbind(polymerization_mapped, missing_polym)
 
   # get HVC input ratio
   HVC_input_ratio <- HVC_input_full / polymerization_full
@@ -164,7 +159,7 @@ calcPlHVCinput <- function(subtype) {
     subtype,
     " of plastics by polymer. ",
     "Data derived from Levi and Cullen 2018 https://doi.org/10.1021/acs.est.7b04573 for available polymers; ",
-    "for PUR, Fibres, Others and Bioplastics based on assumptions regarding their composition."
+    "for 'Other thermosets' based on assumptions regarding their composition."
   )
   output <- list(
     x = x,

R/calcPlLifetime.R

@@ -14,8 +14,6 @@ calcPlLifetime <- function(subtype) {
     rel = sector_map, dim = 3,
     from = "Source", to = "Target"
   )
-  # remove category "Industrial Machinery"
-  data <- data[, , getItems(data, 3) != "Industrial Machinery"]
   description <- paste(
     subtype,
     " of plastic goods by use sector. ",