# Script encoding: UTF-8 # R version: 4.1.1 # Author: Emmanuel Aramendia # This script builds Figures 7, 8, and 9 from the main paper # Libraries --------------------------------------------------------------- library(dplyr) library(tidyr) library(stringr) library(IEATools) library(Recca) library(ggplot2) library(matsindf) library(matsbyname) library(forcats) # Defining list of countries and products --------------------------------- # Selecting which products we want to plot in the paper list_of_products <- c("Coke oven coke", "Electricity", "Wind", "Natural gas", "Gas/diesel oil excl. biofuels", "Motor gasoline", "Crude oil", "Coking coal", "Heat") # Selecting the countries for which we want to plot these products in the paper country_list <- c("Brazil", "China", "EU27", "India", "Mexico", "Russian Federation", "United States", "Turkey") # Figuring out embodied primary energy by product through energy multipliers -------- IEA_io_mats_gma %>% glimpse() # Building kp vector, then multiplying... calc_energy_multipliers_df <- IEA_io_mats_gma %>% dplyr::mutate( kp = kvec_from_template_byname( a = L_pxp %>% select_rows_byname(retain_pattern = paste0("\\[from Resources\\]")), colname = "Product" ), energy_multiplier = matsbyname::matrixproduct_byname( matsbyname::transpose_byname(kp), L_pxp ) %>% matsbyname::transpose_byname(), energy_multiplier_by_product = matsbyname::matrixproduct_byname( matsbyname::hatize_byname(kp), L_pxp ) ) %>% glimpse() # Pulling energy multipliers: energy_multipliers_df <- calc_energy_multipliers_df %>% dplyr::select(Country, Method, Energy.type, Last.stage, Year, energy_multiplier) %>% tidyr::pivot_longer(cols = energy_multiplier, values_to = "Energy_multiplier", names_to = "matnames") %>% matsindf::expand_to_tidy(matvals = "Energy_multiplier", rownames = "Products") %>% dplyr::select(-matnames, -colnames, -rowtypes, -coltypes) %>% dplyr::mutate( Country = stringr::str_extract(Products, "\\{.*\\}"), Country = stringr::str_remove(Country, "\\{"), Country = stringr::str_remove(Country, "\\}"), Products = stringr::str_remove(Products, "\\{.*\\}_") ) %>% dplyr::filter( ! stringr::str_detect(Products, "Resources") ) %>% glimpse() # Filtering products filtered_energy_multipliers_df <- energy_multipliers_df %>% dplyr::filter( Products %in% list_of_products ) %>% dplyr::mutate( `Primary to final energy efficiency` = 1/Energy_multiplier ) %>% dplyr::mutate( Products = dplyr::case_when( Products == "Gas/diesel oil excl. biofuels" ~ "Gas/diesel oil", TRUE ~ Products ) ) # Try to get plots of embodied primary energy by type of energy. # Oil, nat gas, renewables, nuclear. energy_multiplier_by_product_df <- calc_energy_multipliers_df %>% dplyr::select(Country, Method, Energy.type, Last.stage, Year, energy_multiplier_by_product) %>% tidyr::pivot_longer(cols = energy_multiplier_by_product, names_to = "matnames", values_to = "matvals") %>% dplyr::filter(Year == c(2000, 2017)) %>% matsindf::expand_to_tidy() %>% dplyr::select(-matnames, -rowtypes, -coltypes) %>% dplyr::rename( Embodied_Resource = rownames, Embodying_Product = colnames, Embodied_Energy = matvals ) %>% dplyr::mutate( Country = stringr::str_extract(Embodying_Product, "\\{.*\\}"), Country = stringr::str_remove(Country, "\\{"), Country = stringr::str_remove(Country, "\\}"), Embodied_Resource = stringr::str_remove(Embodied_Resource, "\\{.*\\}_"), Embodied_Resource = stringr::str_remove(Embodied_Resource, " \\[from Resources\\]"), Embodying_Product = stringr::str_remove(Embodying_Product, "\\{.*\\}_") ) %>% dplyr::mutate( Embodied_Resource = dplyr::case_when( Embodied_Resource == "Anthracite" ~ "Coal products", Embodied_Resource == "Biodiesels" ~ "Bioenergy", Embodied_Resource == "Biogases" ~ "Bioenergy", Embodied_Resource == "Biogasoline" ~ "Bioenergy", Embodied_Resource == "Coking coal" ~ "Coal products", Embodied_Resource == "Crude oil" ~ "Oil products", Embodied_Resource == "Geothermal" ~ "Renewable energy", Embodied_Resource == "Hydro" ~ "Renewable energy", Embodied_Resource == "Industrial waste" ~ "Other", Embodied_Resource == "Lignite" ~ "Coal products", Embodied_Resource == "Natural gas" ~ "Natural gas", Embodied_Resource == "Natural gas liquids" ~ "Oil products", Embodied_Resource == "Other bituminous coal" ~ "Coal products", Embodied_Resource == "Primary solid biofuels" ~ "Bioenergy", Embodied_Resource == "Solar photovoltaics" ~ "Renewable energy", Embodied_Resource == "Solar thermal" ~ "Renewable energy", Embodied_Resource == "Sub-bituminous coal" ~ "Coal products", Embodied_Resource == "Wind" ~ "Renewable energy", Embodied_Resource == "Heat" ~ "Other", Embodied_Resource == "Nuclear" ~ "Nuclear", Embodied_Resource == "Other liquid biofuels" ~ "Bioenergy", Embodied_Resource == "Municipal waste (non-renewable)" ~ "Other", Embodied_Resource == "Municipal waste (renewable)" ~ "Other", Embodied_Resource == "Other hydrocarbons" ~ "Oil products", Embodied_Resource == "Tide, wave and ocean" ~ "Renewable energy", Embodied_Resource == "Additives/blending components" ~ "Oil products", Embodied_Resource == "Oil shale and oil sands" ~ "Oil products", Embodied_Resource == "Other sources" ~ "Other", Embodied_Resource == "Peat" ~ "Coal products", TRUE ~ Embodied_Resource ) ) %>% dplyr::group_by( Country, Method, Energy.type, Last.stage, Year, Embodied_Resource, Embodying_Product ) %>% dplyr::summarise( Embodied_Energy = sum(Embodied_Energy) ) %>% print() # Defining factor levels order country_order <- c("China", "United States", "EU27", "Russian Federation") # Plot embodied energy by energy type, for different energy products # Figure 8 in paper x11() energy_multiplier_by_product_df %>% dplyr::filter(Country %in% c("China", "EU27", "United States", "Russian Federation")) %>% dplyr::mutate( Country = dplyr::case_when( Country == "Russian Federation" ~ "Russia", TRUE ~ Country ) ) %>% dplyr::filter(Year %in% c(2000, 2017)) %>% dplyr::filter( Embodying_Product %in% c("Electricity", "Heat") ) %>% dplyr::mutate( Year = as.character(Year) ) %>% ggplot(aes(x = Year, y = Embodied_Energy)) + geom_bar(stat = "identity", position = "stack", aes(fill = Embodied_Resource), width = 0.5) + ylab("Embodied primary energy per unit of energy product (w/o unit)") + labs(fill = "Energy source") + scale_fill_viridis_d(option = "turbo", begin = 0.02, end = 0.98) + facet_grid(rows = vars(Country), cols = vars(Embodying_Product)) + theme_bw() + theme(strip.text = element_text(size = 11), legend.title = element_text(size = 12), legend.text = element_text(size = 11), legend.position = "right", axis.title = element_text(size = 12), axis.text.y = element_text(size = 12), axis.text.x = element_text(size = 10)) ggsave("embodied_primary_energy_by_type.pdf", width = 11, height = 10) # Figuring out embodied greenhouse gases emissions ------------------------ # First, defining the emission vectors (emissions by each (product, country), but defined only by product (country doesn't matter)) # (i) For the combustion process define_emission_combustion_vector_df <- calc_energy_multipliers_df %>% tidyr::pivot_longer(cols = -c("Country", "Method", "Energy.type", "Last.stage", "Year"), names_to = "matnames", values_to = "matvalues") %>% dplyr::filter(matnames == "kp") %>% matsindf::expand_to_tidy(matnames = "matnames", matvals = "matvalues", rowtypes = "rowtypes", coltypes = "coltypes") %>% dplyr::mutate( Product_name = stringr::str_remove(rownames, "\\{.*\\}_"), Product_name = stringr::str_remove(Product_name, " \\[from Resources\\]"), matnames = "e_c", matvalues = dplyr::case_when( Product_name == "Anthracite" ~ 26.8, Product_name == "Coking coal" ~ 25.8, Product_name == "Crude oil" ~ 20, Product_name == "Lignite" ~ 27.6, Product_name == "Natural gas" ~ 15.3, Product_name == "Natural gas liquids" ~ 17.5, Product_name == "Other bituminous coal" ~ 25.8, Product_name == "Sub-bituminous coal" ~ 26.2, Product_name == "Other hydrocarbons" ~ 20, Product_name == "Oil shale and oil sands" ~ 29.1, Product_name == "Peat" ~ 28.9, TRUE ~ 0 ) ) %>% print() # Pushing back to matrix format the e_c vector back_to_matrix_format_combustion <- define_emission_combustion_vector_df %>% dplyr::select(-Product_name) %>% dplyr::group_by(Country, Method, Energy.type, Last.stage, Year) %>% matsindf::collapse_to_matrices( matvals = "matvalues" ) %>% dplyr::rename( e_c = matvalues ) %>% glimpse() # (ii) For the extraction process define_emission_extraction_vector_df <- calc_energy_multipliers_df %>% tidyr::pivot_longer(cols = -c("Country", "Method", "Energy.type", "Last.stage", "Year"), names_to = "matnames", values_to = "matvalues") %>% dplyr::filter(matnames == "kp") %>% matsindf::expand_to_tidy(matnames = "matnames", matvals = "matvalues", rowtypes = "rowtypes", coltypes = "coltypes") %>% dplyr::mutate( Product_name = stringr::str_remove(rownames, "\\{.*\\}_"), Product_name = stringr::str_remove(Product_name, " \\[from Resources\\]"), matnames = "e_e", matvalues = dplyr::case_when( Product_name == "Anthracite" ~ 0, Product_name == "Coking coal" ~ 0, Product_name == "Crude oil" ~ 1.7, Product_name == "Lignite" ~ 0, Product_name == "Natural gas" ~ 2.3, Product_name == "Natural gas liquids" ~ 2.7, Product_name == "Other bituminous coal" ~ 0, Product_name == "Sub-bituminous coal" ~ 0, Product_name == "Other hydrocarbons" ~ 1.7, Product_name == "Oil shale and oil sands" ~ 2.5, Product_name == "Peat" ~ 0, TRUE ~ 0 ) ) %>% print() # Pushing back to matrix format the e_e vector back_to_matrix_format_extraction <- define_emission_extraction_vector_df %>% dplyr::select(-Product_name) %>% dplyr::group_by(Country, Method, Energy.type, Last.stage, Year) %>% matsindf::collapse_to_matrices( matvals = "matvalues" ) %>% dplyr::rename( e_e = matvalues ) %>% glimpse() # Binding co2 extension vector, and conducting calculations: calc_co2_multiplier <- calc_energy_multipliers_df %>% dplyr::left_join(back_to_matrix_format_combustion, by = c("Country", "Method", "Energy.type", "Last.stage", "Year")) %>% dplyr::left_join(back_to_matrix_format_extraction, by = c("Country", "Method", "Energy.type", "Last.stage", "Year")) %>% dplyr::mutate( co2_c = matsbyname::transpose_byname(e_c) %>% matsbyname::matrixproduct_byname(matsbyname::hatize_byname(kp)) %>% matsbyname::matrixproduct_byname(L_pxp) %>% matsbyname::transpose_byname(), co2_e = matsbyname::transpose_byname(e_e) %>% matsbyname::matrixproduct_byname(matsbyname::hatize_byname(kp)) %>% matsbyname::matrixproduct_byname(L_pxp) %>% matsbyname::transpose_byname(), ) %>% dplyr::mutate( Z_eiou = matsbyname::matrixproduct_byname(U_EIOU, matsbyname::hatinv_byname(g)), co2_c_eiou = matsbyname::matrixproduct_byname( matsbyname::transpose_byname(co2_c), matsbyname::matrixproduct_byname(Z_eiou, L_ixp) ) %>% matsbyname::transpose_byname(), co2_c_downstream = matsbyname::difference_byname(co2_c, co2_c_eiou) ) %>% glimpse() # Pulling CO2 multipliers co2_multipliers_df <- calc_co2_multiplier %>% dplyr::select(Country, Method, Energy.type, Last.stage, Year, co2_c, co2_e, co2_c_eiou, co2_c_downstream) %>% tidyr::pivot_longer(cols = c("co2_c", "co2_e", "co2_c_eiou", "co2_c_downstream"), values_to = "CO2_multiplier", names_to = "matnames") %>% matsindf::expand_to_tidy(matvals = "CO2_multiplier", rownames = "Products") %>% dplyr::select(-colnames, -rowtypes, -coltypes) %>% dplyr::mutate( Country = stringr::str_extract(Products, "\\{.*\\}"), Country = stringr::str_remove(Country, "\\{"), Country = stringr::str_remove(Country, "\\}"), Products = stringr::str_remove(Products, "\\{.*\\}_") ) %>% dplyr::filter( ! stringr::str_detect(Products, "Resources") ) %>% dplyr::rename( Process = matnames ) %>% glimpse() # Filtering chosen products: filtered_co2_multipliers_df <- co2_multipliers_df %>% dplyr::filter( Products %in% list_of_products ) %>% dplyr::filter(Process %in% c("co2_e", "co2_c_eiou", "co2_c_downstream")) %>% dplyr::mutate( Products = dplyr::case_when( Products == "Gas/diesel oil excl. biofuels" ~ "Gas/diesel oil", TRUE ~ Products ) ) %>% dplyr::filter(Products != "Wind") %>% glimpse() # Plotting CO2 multipliers - Figure 7 in paper x11() filtered_co2_multipliers_df %>% dplyr::mutate( Process = dplyr::case_when( Process == "co2_e" ~ "Methane flaring and leakages", Process == "co2_c_downstream" ~ "Downstream energy use", Process == "co2_c_eiou" ~ "Energy use by the energy industry", ) ) %>% dplyr::filter(Year %in% c(2000, 2017)) %>% dplyr::filter(Country %in% c("China", "EU27", "United States", "Russian Federation")) %>% dplyr::filter(Products %in% c(list_of_products, "Gas/diesel oil")) %>% dplyr::mutate( Products = dplyr::case_when( Products == "Coke oven coke" ~ "Coke", Products == "Gas/diesel oil" ~ "Gas oil", TRUE ~ Products ), Year = as.character(Year), Country = dplyr::case_when( Country == "Russian Federation" ~ "Russia", TRUE ~ Country ) ) %>% ggplot(aes(x = Year, y = CO2_multiplier)) + geom_bar(aes(fill = factor(Process, levels = c("Methane flaring and leakages", "Energy use by the energy industry", "Downstream energy use"))), position = "stack", stat = "identity", width = 0.8) + scale_fill_viridis_d(begin = 0.05, end = 0.5, option = "turbo", direction = -1) + ylab("Energy-related CO2 emissions (kgCO2e/GJ)") + labs( fill = "Process" ) + facet_grid(rows = vars(Country), cols = vars(Products)) + theme_bw() + theme( legend.position = "bottom", legend.text = element_text(size = 11), legend.title = element_text(size = 12), axis.title = element_text(size = 12), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 10.5), strip.text = element_text(size = 11) ) ggsave("co2_multipliers.pdf", width = 11.5, height = 10) # Now, calculate CO2 footprint by sector ---------------------------------- # Creating empty DF product_demand_by_sector_region <- tibble::tibble( Country = character(), Year = numeric(), Y_tau_times_L_pxp = list() ) # Filling in data frame for (region in country_list){ print(region) # Observations to add for region iteration to_add <- IEA_io_mats_gma %>% dplyr::mutate( Y_tau = matsbyname::select_cols_byname( Y, retain_pattern = paste0("\\{", region, "\\}_") ), Y_tau_times_L_pxp = matsbyname::matrixproduct_byname(L_pxp, Y_tau) ) %>% dplyr::mutate( Country = region ) %>% dplyr::select(Country, Year, Y_tau_times_L_pxp) # Adding rows to DF product_demand_by_sector_region <- dplyr::bind_rows( product_demand_by_sector_region, to_add ) } # Calculating CO2 footprint by sector co2_footprint_by_sector <- product_demand_by_sector_region %>% dplyr::mutate( matnames = "Y_tau_times_Lpxp" ) %>% dplyr::rename( matvals = Y_tau_times_L_pxp ) %>% matsindf::expand_to_tidy() %>% dplyr::select(-rowtypes, -coltypes) %>% dplyr::mutate( colnames = stringr::str_remove(colnames, "\\{.*\\}_"), rownames = stringr::str_remove(rownames, "\\{.*\\}_") ) %>% dplyr::rename( Sector = colnames, Product = rownames ) %>% dplyr::filter(stringr::str_detect(Product, "\\[from Resources\\]")) %>% dplyr::mutate( Product = stringr::str_remove(Product, " \\[from Resources\\]") ) %>% dplyr::left_join( define_emission_extraction_vector_df %>% dplyr::select(matvalues, Product_name) %>% dplyr::rename(extraction_intensities = matvalues, Product = Product_name) %>% dplyr::distinct(), by = "Product" ) %>% dplyr::left_join( define_emission_combustion_vector_df %>% dplyr::select(matvalues, Product_name) %>% dplyr::rename(combustion_intensities = matvalues, Product = Product_name) %>% dplyr::distinct(), by = "Product" ) %>% dplyr::mutate( Extraction_emissions = (matvals * extraction_intensities) * 41.868 * 1e3 / 1e9, # Conversion to (1) GJ, and (2) MtCO2e Combustion_emissions = (matvals * combustion_intensities) * 41.868 * 1e3 / 1e9, # Conversion to (1) GJ, and (2) MtCO2e ) %>% glimpse() # List sectors to look at sector_list <- c("Road", "Rail", "Residential", "Chemical and petrochemical", "Iron and steel", "Mining and quarrying", "Non-ferrous metals", "Non-metallic minerals") # Plotting emissions per sector - Figure 9 in paper x11() co2_footprint_by_sector %>% tidyr::pivot_longer( cols = c("Extraction_emissions", "Combustion_emissions"), names_to = "Scope", values_to = "Emissions" ) %>% dplyr::mutate( Product_Group = dplyr::case_when( Product == "Anthracite" ~ "Coal products", Product == "Biodiesels" ~ "Bioenergy", Product == "Biogases" ~ "Bioenergy", Product == "Biogasoline" ~ "Bioenergy", Product == "Coking coal" ~ "Coal products", Product == "Crude oil" ~ "Oil products", Product == "Geothermal" ~ "Renewable energy", Product == "Hydro" ~ "Renewable energy", Product == "Industrial waste" ~ "Other", Product == "Lignite" ~ "Coal products", Product == "Natural gas" ~ "Natural gas", Product == "Natural gas liquids" ~ "Oil products", Product == "Other bituminous coal" ~ "Coal products", Product == "Primary solid biofuels" ~ "Bioenergy", Product == "Solar photovoltaics" ~ "Renewable energy", Product == "Solar thermal" ~ "Renewable energy", Product == "Sub-bituminous coal" ~ "Coal products", Product == "Wind" ~ "Renewable energy", Product == "Heat" ~ "Other", Product == "Nuclear" ~ "Nuclear", Product == "Other liquid biofuels" ~ "Bioenergy", Product == "Municipal waste (non-renewable)" ~ "Other", Product == "Municipal waste (renewable)" ~ "Other", Product == "Other hydrocarbons" ~ "Oil products", Product == "Tide, wave and ocean" ~ "Renewable energy", Product == "Additives/blending components" ~ "Oil products", Product == "Oil shale and oil sands" ~ "Oil products", Product == "Other sources" ~ "Other", Product == "Peat" ~ "Coal products", TRUE ~ Product ) ) %>% dplyr::group_by(Country, Year, Sector, Product_Group) %>% dplyr::summarise(Emissions = sum(Emissions)) %>% dplyr::filter(Emissions != 0) %>% dplyr::filter(Country %in% c("EU27", "China", "United States", "India")) %>% dplyr::filter(Sector %in% c("Road", "Iron and steel", "Chemical and petrochemical")) %>% ggplot(aes(x = Year, y = Emissions)) + geom_area(aes(fill = Product_Group)) + scale_fill_viridis_d(option = "turbo", begin = 0.05, end = 0.3) + ylab("Greenhouse gas emissions (MtCO2e)") + theme_bw() + labs(fill = "Energy source") + theme( legend.position = "bottom", legend.text = element_text(size = 11), legend.title = element_text(size = 12), axis.title = element_text(size = 12), axis.text.y = element_text(size = 11), axis.text.x = element_text(size = 10.5), strip.text = element_text(size = 11) ) + facet_grid( cols = vars(Sector), rows = vars(Country), scales = "free_y" ) ggsave("emissions_by_sector.pdf")