figures.Rmd

title: "figures"
output:
  bookdown::word_document2:
    fig_caption: yes

knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  echo = FALSE,
  comment = "#>",
  fig.path = "../figures/",
  dpi = 300
)

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse,
               here,
               wesanderson,
               patchwork,
               readxl)

pal <- wes_palette("Cavalcanti1", 5, type = "discrete")
extrafont::loadfonts()

# load data wrangling functions
source(here("analysis", "R", "wrangler_functions.R"))

## load result data for EU deciles
eu_q_count = 10

# summary countries aggregated by country quintiles and eu ntile
dat_country_summary_by_cquint_and_euntile = get_country_summary_by_cquint_and_euntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank")
pdat_country_summary_by_cquint_and_euntile =
  pivot_results_longer_adorn(dat_country_summary_by_cquint_and_euntile, cols_ex)

# summary of countries by EU quantile without sectoral resolution
dat_country_summary_by_eu_ntile = get_country_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "eu_q_rank")
pdat_country_summary_by_eu_ntile = 
  pivot_results_longer_adorn(dat_country_summary_by_eu_ntile, cols_ex)

# summary of countries by country quintile with aggregate sectoral resolution
dat_sector_summary_by_country_quintile = get_sector_summary_by_country_quintile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "iso2", "quint", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_country_quintile =
  pivot_results_longer_adorn(dat_sector_summary_by_country_quintile, cols_ex)

# summary of eu ntile with aggregate sectoral resolution
dat_sector_summary_by_eu_ntile = get_sector_summary_by_eu_ntile(eu_q_count)
# pivot to long format for plotting and attach readable indicator names
cols_ex = c("year", "eu_q_rank", "sector_agg_id")
pdat_sector_summary_by_eu_ntile =
  pivot_results_longer_adorn(dat_sector_summary_by_eu_ntile, cols_ex)

Fig. 1

p1 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_fd_me") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Expenditure (trn€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_energy_use_tj") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy footprint (EJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p3 = pdat_country_summary_by_eu_ntile %>%
  filter(year == 2015, 
         indicator == "total_co2eq_kg") %>%
  group_by(eu_q_rank) %>%
  summarise(value = sum(value)*0.000000001,
            eu_ntile_name = first(eu_ntile_name)) %>%
  ggplot(aes(x=eu_ntile_name, y=value)) +
    geom_col(position = position_dodge(), fill=pal[1]) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon footprint (MtCO2eq)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 


p_top = p1 + p2 + p3

p1 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_co2eq_kg)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_co2eq_kg = weighted.mean(pe_co2eq_kg,total_fd_me)),
               aes(y=pe_co2eq_kg, yend=pe_co2eq_kg, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per expenditure (kgCO2eq/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

p2 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  ggplot(aes(x=factor(eu_q_rank), y=pe_energy_use_mj)) +
  geom_violin(aes(weight=total_fd_me), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat_country_summary_by_eu_ntile %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(pe_energy_use_mj = weighted.mean(pe_energy_use_mj,total_fd_me)),
               aes(y=pe_energy_use_mj, yend=pe_energy_use_mj, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Energy intensity per expenditure (MJ/€)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 

dat3 = dat_country_summary_by_eu_ntile %>%
  filter(year == 2015) %>%
  mutate(intensity_e_c = total_co2eq_kg*0.001/total_energy_use_tj)

p3 = dat3 %>%
  ggplot(aes(x=factor(eu_q_rank), y=intensity_e_c)) +
  geom_violin(aes(weight=total_energy_use_tj), fill=pal[1], color=pal[1], alpha=0.5) +
  geom_point( alpha=0.3) +
  geom_segment(data=dat3 %>%
                 filter(year == 2015) %>%
                 group_by(eu_q_rank) %>%
                 summarise(intensity_e_c = weighted.mean(intensity_e_c,total_energy_use_tj)),
               aes(y=intensity_e_c, yend=intensity_e_c, x=eu_q_rank-0.3, xend=eu_q_rank+0.3), size=1.5) +
    #theme_ipsum() +
    theme_minimal() +
    theme(text=element_text(family="Liberation Sans Narrow")) +
    labs(x="", y="Carbon intensity per energy (gCO2eq/TJ)") +
    theme(axis.text.x = element_text(angle = 90)) +
    scale_x_discrete(labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) 


p_bottom = p2 + p1 + p3

#p3 & theme(plot.margin =  margin(1, 1, 1, 1, "mm"))

a = p_top / p_bottom + plot_annotation(tag_levels = 'a') +
  plot_layout(guides = 'collect')  & 
  theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "cm"),
        legend.position = 'bottom',
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

ggsave(here("analysis", "figures", "figure1.pdf"), device=cairo_pdf)

knitr::include_graphics(here::here("analysis", "figures", "figure1.pdf"))

Fig.2

pdat_basket = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), by="sector_agg_id") %>% 
  group_by(eu_q_rank) %>%
  mutate(value = total_fd_me/sum(total_fd_me)*100) %>%
  select(five_sectors, eu_q_rank, 
         value) %>%
  mutate(indicator = "Expenditure share (%)")

pdat_int_co2eq = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), by="sector_agg_id") %>% 
  mutate(value = (total_co2eq_kg)/(total_fd_me*1000000)) %>%
  select(five_sectors, eu_q_rank, 
         value) %>%
  mutate(indicator = "Carbon intensity (kgCO2eq/€)")

pdat_int_energy = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), by="sector_agg_id") %>% 
  mutate(value = (total_energy_use_tj)/(total_fd_me)) %>%
  select(five_sectors, eu_q_rank,
         value) %>%
  mutate(indicator = "Energy intensity (MJ/€)")

library(viridis)
pdat = pdat_int_co2eq %>%
  bind_rows(pdat_int_energy) %>%
  bind_rows(pdat_basket) 


p1 = pdat_basket %>%
  mutate(Decile = if_else(eu_q_rank<10, paste0("D0", eu_q_rank), paste0("D", eu_q_rank))) %>% 
ggplot(aes(x=value, y=five_sectors, color=Decile)) +
  #geom_point(shape=18, alpha=0.75) +
  geom_text(label="I", alpha=0.75, fontface="bold") +
  scale_color_viridis(option = "A", end = 0.8,
                      direction = -1, 
                      discrete = T,
                      name="European\ndecile") +
  theme_minimal() +
  labs(x=unique(pdat_basket$indicator),y="")

p2 = pdat_int_energy %>%
  mutate(Decile = if_else(eu_q_rank<10, paste0("D0", eu_q_rank), paste0("D", eu_q_rank))) %>% 
ggplot(aes(x=value, y=five_sectors, color=Decile)) +
  #geom_point(shape=18, alpha=0.75) +
  geom_text(label="I", alpha=0.75, fontface="bold") +
  scale_color_viridis(option = "A", end = 0.8,
                      direction = -1, 
                      discrete = T,
                      name="European\ndecile") +
  theme_minimal() +
  labs(x=unique(pdat_int_energy$indicator),y="")

p3 = pdat_int_co2eq %>%
  mutate(Decile = if_else(eu_q_rank<10, paste0("D0", eu_q_rank), paste0("D", eu_q_rank))) %>% 
ggplot(aes(x=value, y=five_sectors, color=Decile)) +
  #geom_point(shape=18, alpha=0.75) +
  geom_text(label="I", alpha=0.75, fontface="bold") +
  scale_color_viridis(option = "A", end = 0.8,
                      direction = -1, 
                      discrete = T,
                      name="European\ndecile") +
  theme_minimal() +
  labs(x=unique(pdat_int_co2eq$indicator),y="")

#plot_annotation(tag_levels = 'a')
a = (p1 + p2 + p3) + 
  plot_layout(guides = "collect") + 
  plot_annotation(tag_levels = 'a') & 
  theme(legend.position = 'bottom',
        text=element_text(family="Liberation Sans Narrow"),
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

ggsave(here("analysis", "figures", "figure2.pdf"), device=cairo_pdf)

knitr::include_graphics(here::here("analysis", "figures", "figure2.pdf"))

Fig.3

pdat = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), by="sector_agg_id") %>% 
  select(five_sectors, eu_q_rank, contains("pae_co2eq"), contains("pae_energy")) %>%
  pivot_longer(cols=-c(five_sectors, eu_q_rank), names_to="indicator", values_to="value") %>%
  mutate(indicator_type = if_else(str_detect(indicator, "co2eq"), 
                                  "tCO2eq per adult eq", 
                                  "Energy GJ per adult eq")) %>%
  filter(indicator != "pae_co2eq_t", 
         indicator != "pae_energy_use_gj")
  
  
pal <- wes_palette("Zissou1", 4, type = "discrete")
pal = pal[c(2,1,3,4)]
  
p1 = ggplot(pdat %>% filter(indicator_type == "tCO2eq per adult eq"), 
       aes(x=(eu_q_rank), y=value, fill=indicator)) +
  geom_col(position = position_stack()) +
  facet_wrap(~five_sectors, ncol = 5) +
  scale_fill_manual(values = pal, name="Location", labels=c("Direct", "Domestic", "Europe", "non-Europe")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10), 
                     labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) +
  #theme_ipsum() +
  theme_minimal() +
   labs(x="", y="Carbon footprint tCO2eq per adult eq")# +
  # theme(legend.position = "bottom", 
  #       text=element_text(family="Liberation Sans Narrow"),
  #         axis.text.x = element_text(angle = 90, size = 12),
  #         axis.text.y = element_text(size = 12),
  #         strip.text = element_text(size = 12),
  #         axis.title.y = element_text(size = 10))


p2 = ggplot(pdat %>% filter(indicator_type == "Energy GJ per adult eq"), 
       aes(x=(eu_q_rank), y=value, fill=indicator)) +
  geom_col(position = position_stack()) +
  facet_wrap(~five_sectors, ncol = 5) +
  scale_fill_manual(values = pal, name="Location", labels=c("Direct", "Domestic", "Europe", "non-Europe")) +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10), 
                     labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) +
  #theme_ipsum() +
  theme_minimal() +
  #theme(text=element_text(family="Liberation Sans Narrow")) +
  labs(x="", y="Energy footprint GJ per adult eq") #+
  # theme(legend.position = "bottom", 
  #         axis.text.x = element_text(angle = 90, size = 12),
  #         axis.text.y = element_text(size = 12),
  #         strip.text = element_text(size = 12),
  #         axis.title.y = element_text(size = 10))

a = (p1 / p2) + 
  plot_layout(guides = "collect") & 
  theme(text=element_text(family="Liberation Sans Narrow"),
        legend.position = "bottom",
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(angle = 90, size = 11),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        strip.text = element_text(size = 13),
        legend.title = element_text(size=13),
        plot.margin =  margin(1, 1, 1, 1, "mm"),
        panel.spacing = unit(c(1, 1, 1, 1), "mm"))


# hm = pdat %>% filter(indicator_type == "Energy GJ per adult eq") %>%
#   group_by(eu_q_rank) %>%
#   summarise(value = sum(value))

ggsave(here("analysis", "figures", "figure3.pdf"))

knitr::include_graphics(here::here("analysis", "figures", "figure3.pdf"))

Fig.4

pal <- wes_palette("Zissou1", 4, type = "discrete")
pal = pal[c(2,1,3,4)]

pdat_int_energy = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>% 
  mutate(intensity_energy = (total_energy_use_tj)/(total_fd_me)) %>%
  select(five_sectors, eu_q_rank, 
         intensity_energy) %>%
  filter(eu_q_rank == 10) %>%
  select(-eu_q_rank)

pdat_final_demand = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>%
  left_join(pdat_int_energy, by="five_sectors") %>%
  mutate(total_energy_use_tj_new = (total_fd_me)*intensity_energy) %>%
  mutate(total_energy_use_tj_diff = total_energy_use_tj-total_energy_use_tj_new) %>%
  select(five_sectors, eu_q_rank, 
         total_energy_use_tj_diff, total_energy_use_tj_new) %>%
  group_by(eu_q_rank) %>%
  summarise(total_energy_use_tj_new = sum(total_energy_use_tj_new)*0.000001,
            total_energy_use_tj_diff = sum(total_energy_use_tj_diff)*0.000001) %>%
  pivot_longer(-c(eu_q_rank), names_to = "indicator", values_to = "value")


p1 = ggplot(pdat_final_demand, aes(x=eu_q_rank, y=value, fill=indicator, alpha=indicator)) +
  geom_col(position = position_stack()) +
  scale_fill_manual(values=c(pal[1], pal[2]), 
                    labels=c("2015", "Best technology"), name="Energy\nfootprint") +
  scale_alpha_manual(values=c(0.3,1),labels=c("2015", "Best technology"), name="Energy\nfootprint")+
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10), 
                     labels = c("D01","D02","D03","D04","D05","D06","D07","D08","D09","D10")) +
  labs(y="Energy footprint (EJ)", x="") +
  theme_minimal() +
  coord_flip() +
  theme(legend.position="bottom")

pdat_energy_country = get_sector_summary_by_country_quintile_direct() %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>%
  left_join(pdat_int_energy, by="five_sectors") %>%
  mutate(total_energy_use_tj_new = (total_fd_me)*intensity_energy) %>%
  #mutate(total_energy_use_tj_save = total_energy_use_tj_new/total_energy_use_tj*100) %>%
  select(iso2, 
         total_energy_use_tj_new, total_energy_use_tj) %>%
  group_by(iso2) %>%
  summarise(total_energy_use_tj_save = 100 - sum(total_energy_use_tj_new)/sum(total_energy_use_tj)*100)

library(rworldmap)
library(ggthemes)



df_country_mean = pdat_energy_country %>%
  mutate(iso2 = if_else(iso2 == "EL", "GR", iso2)) %>%
  mutate(iso2 = if_else(iso2 == "UK", "GB", iso2)) %>%
  left_join(ISO_3166_1 %>% select(iso2 = Alpha_2, iso3 = Alpha_3), by="iso2")


quantile_rank_map = joinCountryData2Map(df_country_mean, 
                                        joinCode = "ISO_A3", 
                                        nameJoinColumn = "iso3")
quantile_rank_map_poly = fortify(quantile_rank_map) #extract polygons 

quantile_rank_map_poly = merge(quantile_rank_map_poly, quantile_rank_map@data, by.x="id", by.y="ADMIN", all.x=T)
quantile_rank_map_poly = quantile_rank_map_poly %>% arrange(id, order)

map = ggplot() + 
  geom_polygon(data = quantile_rank_map_poly, aes(long, lat, group = group), color="black") +
  geom_polygon(data = quantile_rank_map_poly, aes(long, lat, group = group, 
                                        fill=total_energy_use_tj_save)) + 
  #scale_fill_gradient(name="Mean decile \nrank", low="skyblue", high = "skyblue4", na.value = "#EEEEEE") +
  scale_fill_viridis(direction = -1, discrete = F, na.value = "white", name="Energy\nsavings (%)") +
  theme_map() +
  labs(x="",y="")+
  coord_map("bonne", lat0 = 50,xlim = c(-9, 40), ylim = c(38, 68), clip="on") +
  theme(legend.position = c(0.85, 0.6))
  

a = p1 + map + plot_annotation(tag_levels = 'a') +
  theme(plot.margin =  margin(1, 1, 1, 1, "mm"),
        legend.position = 'bottom',
        text=element_text(family="Liberation Sans Narrow"),
        axis.title.y = element_text(size=13, hjust = 0.5),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size=12),
        legend.title = element_text(size=13))

ggsave(here("analysis", "figures", "figure4.pdf"))

knitr::include_graphics(here::here("analysis", "figures", "figure4.pdf"))

Fig.5

pdat_int_energy = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>% 
  mutate(intensity_energy = (total_energy_use_tj)/(total_fd_me)) %>%
  select(five_sectors, eu_q_rank, 
         intensity_energy) %>%
  filter(eu_q_rank == 10) %>%
  select(-eu_q_rank)

pdat_final_demand = get_sector_summary_by_eu_ntile_direct(eu_q_count) %>%
  ungroup() %>%
  filter(year==2015) %>%
  left_join(read_csv(here("analysis/data/derived/sectors_agg_method1_ixi.csv")), 
            by="sector_agg_id") %>%
  left_join(pdat_int_energy, by="five_sectors") %>%
  mutate(total_energy_use_tj_new = (total_fd_me)*intensity_energy) %>%
  mutate(total_energy_use_tj_diff = total_energy_use_tj-total_energy_use_tj_new) %>%
  select(eu_q_rank,total_energy_use_tj_new) %>%
  group_by(eu_q_rank) %>%
  summarise(total_energy_use_tj_new = sum(total_energy_use_tj_new)) %>%
  mutate(pae_energy_use_tj = total_energy_use_tj_new/33417583,
         pae_energy_use_gj = pae_energy_use_tj*1000)

df_energy_deciles = pdat_final_demand %>%
  select(eu_q_rank, pae_energy_use_gj)

ineq_curr = df_energy_deciles$pae_energy_use_gj[10]/df_energy_deciles$pae_energy_use_gj[1]

df_scenario_info = read_excel(here("analysis/data/raw/scenarios.xlsx"), sheet="overview") %>%
  select(scenario, fe_gj_aeu = final_energy_gj_per_aeu_2050,
         ccs_required = primary_energy_fossil_w_ccs2050_ej,
         description) %>%
  arrange(fe_gj_aeu) %>%
  mutate(fe_gj_aeu = round(fe_gj_aeu),
         ccs_required = round(ccs_required)) 

mea_seq = c(df_scenario_info$fe_gj_aeu, seq(from=50, to=300, by=50))
#mea = c(16,300)
mer = c(16,60)

# vectorized function that returns scaled quantiles given 
#quantile index column and column with quantile averages
#qidx: quantile index
#qavg_to_scale: column to scale
#first_target: target value of first quantile
#mean_target: target mean of scaled quantiles

scaled_quantiles <- function(.data, 
                          qidx, 
                          qavg_to_scale, 
                          first_target, 
                          mean_target) {
  
  # cumbersomely extract current quantile mean
  mean_current = .data %>%
    ungroup() %>%
    summarise(mean_cur = first(mean({{qavg_to_scale}}))) %>%
    pull(mean_cur)
    
  # cumbersomely extract current first wuantile value
  first_current = .data %>%
    ungroup() %>%
    arrange({{qidx}}) %>%
    summarise(first_cur = first({{qavg_to_scale}})) %>%
    pull(first_cur)
  
  df_tmp = .data %>%
    mutate(tmp = {{qavg_to_scale}}*mean_target/mean_current)
  
  first_tmp = df_tmp$tmp[1]
  
  df_tmp = df_tmp %>%
    mutate(scaled = mean_target-(mean_target-tmp) * (mean_target-first_target)/(mean_target-first_tmp)) %>%
    select({{qidx}}, scaled) %>%
    mutate(v_mean = mean_target, v_first = first_target)
}

## run once to save file
df_all = NULL
for (min_energy in seq(from=mer[1], to=mer[2], by=0.1)) {
  for (mean_energy in mea_seq) {
    if (min_energy <= mean_energy) {
      df_all = df_all %>%
      bind_rows(df_energy_deciles %>%
                  scaled_quantiles(eu_q_rank, pae_energy_use_gj, min_energy, mean_energy))
    }
  }
}

saveRDS(df_all, here("analysis/data/derived/scenarios_extrafine.rds"))

round_by = 10

df_all = readRDS(here("analysis/data/derived/scenarios_extrafine.rds")) %>%
  filter(eu_q_rank %in% c(1,10)) %>%
  group_by(v_mean, v_first) %>%
  summarise(ratio = last(scaled)/first(scaled)) %>%
  mutate(bin_ratio = if_else((ratio*100)%%round_by > round_by*0.5, 
                             ratio*100+(round_by-(ratio*100)%%round_by),
                             ratio*100-(ratio*100)%%round_by)) %>%
  group_by(bin_ratio, v_first) %>%
  summarise(v_mean = mean(v_mean)) %>%
  mutate(bin_ratio = bin_ratio*0.01)

df_grid = df_all %>%
  filter(!(v_mean %in% df_scenario_info$fe_gj_aeu))
df_scenario = df_all %>%
  filter((v_mean %in% df_scenario_info$fe_gj_aeu))

df_dle = tibble(bin_ratio = 1, v_first = 16, v_mean = 16)

scenario_label = df_scenario_info %>%
  mutate(labl = paste0(scenario, " (", fe_gj_aeu, ")")) %>%
  pull(labl)

library(ggsci)

a = df_grid %>%
  ggplot(aes(x=v_first, y=bin_ratio, group=v_mean)) +
  geom_smooth(aes(linetype="Maximum energy\nsupply (GJ/ae)"), se=FALSE, color="grey", size=0.5) +
  scale_linetype_manual(name="", values = c(2)) +
  geom_smooth(data=df_scenario, aes(color=factor(v_mean)), se=FALSE) +
  geom_point(data=df_dle, aes(color=factor(v_mean))) +
  geom_hline(yintercept = ineq_curr, color="grey") +
  annotate(geom="text",x=56,y=12,label="300", angle=-35, size=3, color="grey32") +
  annotate(geom="text",x=46,y=9.5,label="200", angle=-35, size=3, color="grey32") +
  annotate(geom="text",x=33.5,y=6,label="100", angle=-35, size=3, color="grey32") +
  annotate(geom="text",x=60,y=ineq_curr-0.28,
           label="energy use inequality 2015", size=3.75, hjust=1, color="grey40") +
  scale_color_npg(
                     name = "Scenario", 
                     labels = scenario_label) +
  lims(x=c(15.5,60), y=c(1,13)) +
  labs(x= "Minimum energy use (GJ/ae)", y="Energy use inequality (10:10 ratio)") +
  theme_minimal() +
  theme(text=element_text(family="Liberation Sans Narrow"),
      axis.title.y = element_text(size=13, hjust = 0.5),
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      legend.text = element_text(size=12),
      legend.title = element_text(size=13))


ggsave(here("analysis", "figures", "figure5.pdf"))

knitr::include_graphics(here::here("analysis", "figures", "figure5.pdf"))