figures.Rmd

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

```{r setup, echo = FALSE, include = FALSE, message = FALSE}
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()

```


```{r functions}
# load data wrangling functions
source(here("analysis", "R", "wrangler_functions.R"))

```

```{r load-data, include=FALSE}

## 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

```{r ntiles-total}

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

```

```{r ntiles-intensity-violin}


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"))
```

```{r , fig.width=12, fig.height=8}

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)

```

```{r figure1, out.width="98%", fig.cap="Household expenditure and environmental footprints and intensities across European expenditure deciles. Total expenditures (a), energy footprint (b), and carbon footprint (c) per decile. Energy intensity of consumption as energy footprint per expenditure (d), carbon intensity of consumption as carbon footprint per expenditure (e), and carbon intensity of energy as carbon footprint per energy footprint (f)."}
knitr::include_graphics(here::here("analysis", "figures", "figure1.pdf"))
```

# Fig.2

```{r , fig.width=8, fig.height=3}

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)
```

```{r figure2, out.width="100%", fig.cap="Sectoral expenditure shares and carbon intensities of European expenditure deciles. Share of expenditure per final demand sector of total spending per decile in percent (a) and carbon intensity per final demand sector in kgCO2eq/€."}
knitr::include_graphics(here::here("analysis", "figures", "figure2.pdf"))
```

# Fig.3

```{r, fig.width=10, fig.height=6}

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"))
```

```{r figure3, out.width="100%", fig.cap="Energy and carbon footprints by final consumption sector and European expenditure decile in 2015, further broken down by emission source location."}
knitr::include_graphics(here::here("analysis", "figures", "figure3.pdf"))
```

# Fig.4

```{r}
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")

```

```{r}

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)

```

```{r include=FALSE}

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)

```

```{r dle-country, fig.width=8}
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))
  
```

```{r fig.height=7, fig.width=11}

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"))
```

```{r figure4, out.width="100%", fig.cap="Energy savings through sectoral best current technology by expenditure deciles a) and country b)."}
knitr::include_graphics(here::here("analysis", "figures", "figure4.pdf"))
```

# Fig.5

```{r }

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)


```


```{r eval=FALSE, calc-scenarios}

# 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"))

```


```{r}

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)


```


```{r, fig.width=8, fig.height=4}

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 (90: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"))
```


```{r figure5, out.width="70%", fig.align="center", fig.cap="Mean energy available for Europe in decarbonisation scenarios, positioned in option space between a range of minimum energy requirements and range of associated maximum inequality. All expenditure deciles have 'best technology' already."}
knitr::include_graphics(here::here("analysis", "figures", "figure5.pdf"))
```