main.Rmd

---
title: "The Green Hydrogen Ambition and Implementation Gap"
author: "Adrian Odenweller"
date: "2024-11-05"
output: html_document
---

# Setup and configurations

```{r setup, include=FALSE}
# Load required packages
library(tidyverse)  # Data wrangling
library(quitte)  # Data wrangling
library(readxl)  # Reading Excel files
library(cowplot)  # Plot themes and arranging plots
library(ggnewscale)  # Separate legends
library(ggsci)  # Colour palette
library(ggpp)  # For function position_stacknudge
library(ggrepel)  # Adding annotations
library(ggforce)  # Zooming into sub-parts of plots
library(ggsankey)  # Sankey diagrams
library(ggbrace)  # Creating braces in plots
library(svglite)  # Vector graphics
library(ggpattern)  # Plot hatching
library(RColorBrewer)  # Colour palette
library(zoo)  # na.approx function
library(tools)  # file_ext

# Configure plot theme
font.size <- 6
theme_set(theme_cowplot(font_size = font.size))
theme_update(
  plot.title = element_text(size = font.size, hjust = 0.5, face = "plain"),
  axis.line = element_line(linewidth = 0.25),
  axis.ticks = element_line(linewidth = 0.25)
)

# Source functions
source("scripts/readIEAHydrogenProjectsDB.R")
source("scripts/readIEAHydrogenProjectsDBWithEndUse.R")
source("scripts/calcProjectTracking.R")
source("scripts/calcProjectTracking2023.R")
source("scripts/plotProjectTracking.R")
source("scripts/calcRates.R")
source("scripts/plotRates.R")
source("scripts/readScenarios.R")
source("scripts/readTechCalcLearning.R")
source("scripts/plotImplementationGap2023.R")
source("scripts/plotProjectTracking2023.R")
source("scripts/plotDatabaseDelta.R")
source("scripts/plotCostGap.R")
source("scripts/plotCostGapRange.R")
source("scripts/plotRequiredSubsidies.R")
source("scripts/plotSubsidyGap.R")
source("scripts/plotRequiredSubsidies1.5Cmedian.R")
source("scripts/readIAMscenarios.R")

# Don't include code when knitting
knitr::opts_chunk$set(echo = FALSE)

```

# Read IEA Hydrogen Projects Database

## Version 2021

### Original

```{r data v2021 orig, warning=FALSE, message=FALSE}

file.v2021.orig <-
  "data/IEA_H2_DB_2021-10/original/IEA Hydrogen Projects Database 2021 (revised).xlsx"

temp <- readIEAHydrogenProjectsDB(
  file.h2projects = file.v2021.orig,
  h2projects.range = "A5:AF994",
  end.year = 2021
)

data.v2021.orig.proj <- temp[[1]]
data.v2021.orig <- temp[[2]] %>%
  mutate(version = "v2021")

```

### Validated

```{r data v2021 validated, warning=FALSE, message=FALSE}

file.v2021.checked <-
  "data/IEA_H2_DB_2021-10/quality_checked/IEA_Hydrogen_Projects_Database_2021_quality_checked.xlsx"

temp <- readIEAHydrogenProjectsDB(
  file.h2projects = file.v2021.checked,
  h2projects.range = "A5:AG994",
  end.year = 2021,
  checked = "TRUE"
)

data.v2021.proj <- temp[[1]]
data.v2021 <- temp[[2]] %>%
  mutate(version = "v2021")

```

## Version 2022

### Original

```{r data v2022 orig, warning=FALSE, message=FALSE}

file.v2022.orig <-
  "data/IEA_H2_DB_2022-10/original/Hydrogen projects database public version.xlsx"

temp <- readIEAHydrogenProjectsDB(
  file.h2projects = file.v2022.orig,
  h2projects.range = "A5:AF1482",
  end.year = 2022
)

data.v2022.orig.proj <- temp[[1]]
data.v2022.orig <- temp[[2]] %>%
  mutate(version = "v2022")

```

### Validated

```{r data v2022 validated, warning=FALSE, message=FALSE}

file.v2022.checked <-
  "data/IEA_H2_DB_2022-10/quality_checked/IEA_Hydrogen_Projects_Database_2022_quality_checked.xlsx"

temp <- readIEAHydrogenProjectsDB(
  file.h2projects = file.v2022.checked,
  h2projects.range = "A5:AG1482",
  end.year = 2022,
  checked = TRUE
)

data.v2022.proj <- temp[[1]]
data.v2022 <- temp[[2]] %>%
  mutate(version = "v2022")

```

## Version 2023

### Original

```{r data v2023 orig, warning=FALSE, message=FALSE}

file.v2023.orig <-
  "data/IEA_H2_DB_2023-10/original/Hydrogen_production_projects_corrected.xlsx"

temp <- readIEAHydrogenProjectsDB(
  file.h2projects = file.v2023.orig,
  h2projects.range = "A5:AF2001",
  end.year = 2023
)

data.v2023.orig.proj <- temp[[1]]
data.v2023.orig <- temp[[2]] %>%
  mutate(version = "v2023")

```

### Validated

```{r data v2023 validated, warning=FALSE, message=FALSE}

file.v2023.checked <-
  "data/IEA_H2_DB_2023-10/quality_checked/IEA_Hydrogen_Projects_Database_2023_quality_checked.xlsx"

temp <- readIEAHydrogenProjectsDB(
  file.h2projects = file.v2023.checked,
  h2projects.range = "A5:AG2005",
  end.year = 2023,
  checked = TRUE
)

data.v2023.proj <- temp[[1]]
data.v2023 <- temp[[2]] %>%
  mutate(version = "v2023")

```

# Outcome of 2023 projects

```{r data outcome 2023, warning=FALSE, message=FALSE}

file.v2023outcome <- 
  "data/IEA_H2_DB_2023-10_2023outcome/IEA_Hydrogen_Projects_Database_2023_only2023_outcome.xlsx"

temp <- readIEAHydrogenProjectsDB(
  file.h2projects = file.v2023outcome,
  h2projects.range = "A5:AG128",
  end.year = 2024,
  checked = TRUE
)

data.v2023outcome.proj <- temp[[1]]
data.v2023outcome <- temp[[2]] %>% 
  mutate(version = "v2023outcome")

```

# Read electrolysis capacity requirements in 1.5°C scenarios

## Institutional scenarios

```{r read institutional scenarios}

data.scenarios <- readScenarios("data/various_scenarios_greenh2.xlsx")

```

## IAM scenarios (IPCC AR6 and NFGS)

```{r read iam scenarios}

# Get realised capacity in 2023 to filter scenarios
cumcap.2023 <- data.v2023 %>%
  filter(status == "Operational",
         year == 2022) %>% 
  pull(cumcap.sum) %>% 
  sum() +
  data.v2023outcome %>%
  filter(year == 2023) %>%
  summarise(value = sum(cap.sum)) %>% 
  pull()

# IPCC AR6
data.ar6 <- readIAMscenarios("data/ipcc_ar6_c1_greenh2.csv",
                             cap2023 = cumcap.2023)

# NGFS
data.ngfs <- readIAMscenarios("data/ngfs_greenh2.xlsx",
                              cap2023 = cumcap.2023)

```

# Overview

## FIGURE 1: Overview of ambition and implementation gap

```{r figure 1}

range.left <- seq(2020, 2024)
range.right <- seq(2024, 2030)

# Get installed capacity from IEA v2023 and outcome of 2023 projects
data.real <- data.v2023 %>% 
  filter(status == "Operational",
         year %in% seq(2020, 2022)) %>% 
  group_by(year) %>% 
  summarise(value = sum(cumcap.sum))

# Add 2023 value (cumulative sum doesn't work here as this only contains 2023 projects)
data.real <- data.real %>% 
  add_row(year = 2023,
          value =
            filter(data.real, year == 2022)$value + 
            data.v2023outcome %>% filter(year == 2023) %>% summarise(value = sum(cap.sum)) %>% pull())

# Announcements in 2021 (IEA v2021)
data.v2021.plot.left <- data.v2021 %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Announcements in 2022 (IEA v2022)
data.v2022.plot.left <- data.v2022 %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Announcements in 2023 (IEA v2023)
data.v2023.plot.left <- data.v2023 %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

legend.labels <- c("v2021" = "Announcements by 2021",
                   "v2022" = "Announcements by 2022",
                   "v2023" = "Announcements by 2023")

# Plot until 2024
p.left <- ggplot() +
  # Project announcements
  geom_line(data = data.v2021.plot.left,
            mapping = aes(x = year, y = sum, color = "v2021", alpha = "v2021")) +
  geom_point(data = data.v2021.plot.left,
             mapping = aes(x = year, y = sum, color = "v2021", alpha = "v2021"),
             size = 0.75) +
  geom_line(data = data.v2022.plot.left,
            mapping = aes(x = year, y = sum, color = "v2022", alpha = "v2022")) +
  geom_point(data = data.v2022.plot.left,
             mapping = aes(x = year, y = sum, color = "v2022", alpha = "v2022"),
             size = 0.75) +
  geom_line(data = data.v2023.plot.left,
            mapping = aes(x = year, y = sum, color = "v2023", alpha = "v2023")) +
  geom_point(data = data.v2023.plot.left,
             mapping = aes(x = year, y = sum, color = "v2023", alpha = "v2023"),
             size = 0.75) +
  scale_color_npg(name = "Project announcements",
                  guide = guide_legend(order = 1),
                  breaks = c("v2021",
                             "v2022",
                             "v2023"),
                  labels = legend.labels) +
  scale_alpha_manual(name = "Project announcements",
                     guide = guide_legend(order = 1),
                     breaks = c("v2021",
                                "v2022",
                                "v2023"),
                     values = c("v2021" = 0.4,
                                "v2022" = 0.7,
                                "v2023" = 1),
                     labels = legend.labels) +
  new_scale_color() +
  # Realised projects
  geom_line(data = data.real,
            mapping = aes(x = year, y = value, color = "Realised projects")) +
  geom_point(data = data.real,
             mapping = aes(x = year, y = value, color = "Realised projects"),
             size = 1) +
  scale_color_manual(name = NULL, values = c("Realised projects" = "black"),
                     guide = guide_legend(order = 2)) +
  # Horizontal dashed lines (2022)
  geom_segment(
    mapping = aes(
      x = 2022,
      xend = 2023.4,
      y = data.real %>% filter(year == 2022) %>% pull(value),
      yend = data.real %>% filter(year == 2022) %>% pull(value)),
    linetype = "dashed",
    linewidth = 0.25
  ) +
  geom_segment(
    mapping = aes(
      x = 2022,
      xend = 2023.4,
      y = data.v2021.plot.left %>% filter(year == 2022) %>% pull(sum),
      yend = data.v2021.plot.left %>% filter(year == 2022) %>% pull(sum)),
    linetype = "dashed",
    linewidth = 0.25
  ) +
  # Horizontal dashed lines (2023)
  geom_segment(
    mapping = aes(
      x = 2023,
      xend = 2023.9,
      y = data.real %>% filter(year == 2023) %>% pull(value),
      yend = data.real %>% filter(year == 2023) %>% pull(value)),
    linetype = "dashed",
    linewidth = 0.25
  ) +
  geom_segment(
    mapping = aes(
      x = 2023,
      xend = 2023.9,
      y = data.v2021.plot.left%>% filter(year == 2023) %>% pull(sum),
      yend = data.v2021.plot.left %>% filter(year == 2023) %>% pull(sum)),
    linetype = "dashed",
    linewidth = 0.25
  ) +
  # Implementation gap (2022)
  geom_segment(
    aes(x = 2023.25,
        xend = 2023.25,
        y = data.real %>% filter(year == 2022) %>% pull(value),
        yend = data.v2021.plot.left %>% filter(year == 2022) %>% pull(sum)),
    arrow = arrow(length = unit(1.5, "mm"), ends = "both"),
    linewidth = 0.5
  ) +
  geom_segment(
    aes(x = 2023.75,
        xend = 2023.75,
        y = data.real %>% filter(year == 2023) %>% pull(value),
        yend = data.v2021.plot.left %>% filter(year == 2023) %>% pull(sum)),
    arrow = arrow(length = unit(1.5, "mm"), ends = "both"),
    linewidth = 0.5
    ) +
  # Annotation
  geom_text(
    aes(x = 2023.25,
        y = mean(c(data.real %>% filter(year == 2022) %>% pull(value),
                   data.v2021.plot.left %>% filter(year == 2022) %>% pull(sum)))),
    label = "2022",
    size = 5 / .pt,
    hjust = 0,
    vjust = 0.5,
    nudge_x = 0.05
  ) +
  geom_text(
    aes(x = 2023.75,
        y = mean(c(data.real %>% filter(year == 2023) %>% pull(value),
                   data.v2021.plot.left %>% filter(year == 2023) %>% pull(sum)))),
    label = "2023",
    size = 5 / .pt,
    hjust = 0,
    vjust = 0.5,
    nudge_x = 0.05
  ) +
  geom_text(
    aes(x = 2024,
        y = 1),
    label = "\u2460 Past\nimplementation\ngap",
    size = 6 / .pt,
    hjust = 0,
    lineheight = 1,
    fontface = "bold"
  ) +
  # Formatting
  xlab("Year") +
  ylab("Capacity (GW)") +
  ggtitle("Past green hydrogen implementation gap") +
  theme(legend.position = "inside",
        legend.position.inside = c(0.1,0.9),
        legend.justification = c("left","top"),
        plot.margin = margin(t = 1, r = 10, b = 1, l = 1, unit = "mm")) +
  coord_cartesian(xlim = c(2020,2024), clip = "off")

# Expectation in 2021 (IEA v2021)
data.v2021.plot.right <- data.v2021 %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Expectation in 2022 (IEA v2022)
data.v2022.plot.right <- data.v2022 %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Expectation in 2023 (IEA v2023)
data.v2023.plot.right <- data.v2023 %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Scenarios in 2030
data.plot.scenarios <- data.scenarios %>%
  filter(organisation == "IEA", year == 2030)

# Upscaling backed by subsidies
data.backed <- read_csv("data/fig1_capacity_backed_by_subsidies.csv") %>% 
  filter(year <= 2030) %>% 
  pivot_wider(names_from = scenario, values_from = value)

# Plot from 2024 to 2030
p.right <- ggplot() +
  # Project announcements
  geom_line(data = data.v2021.plot.right,
            mapping = aes(x = year, y = sum, color = "v2021", alpha = "v2021")) +
  geom_point(data = data.v2021.plot.right,
             mapping = aes(x = year, y = sum, color = "v2021", alpha = "v2021"),
             size = 0.75) +
  geom_line(data = data.v2022.plot.right,
            mapping = aes(x = year, y = sum, color = "v2022", alpha = "v2022")) +
  geom_point(data = data.v2022.plot.right,
             mapping = aes(x = year, y = sum, color = "v2022", alpha = "v2022"),
             size = 0.75) +
  geom_line(data = data.v2023.plot.right,
            mapping = aes(x = year, y = sum, color = "v2023", alpha = "v2023")) +
  geom_point(data = data.v2023.plot.right,
             mapping = aes(x = year, y = sum, color = "v2023", alpha = "v2023"),
             size = 0.75) +
  scale_color_npg(name = "Project announcements",
                  guide = guide_legend(order = 2),
                  breaks = c("v2021",
                             "v2022",
                             "v2023"),
                  labels = legend.labels) +
  scale_alpha_manual(name = "Project announcements",
                     guide = guide_legend(order = 2),
                     breaks = c("v2021",
                                "v2022",
                                "v2023"),
                     values = c("v2021" = 0.4,
                                "v2022" = 0.7,
                                "v2023" = 1),
                     labels = legend.labels) +
  # Range of IEA scenarios
  geom_errorbar(aes(x = 2030,
                    ymin = data.plot.scenarios %>% pull(value) %>% min(),
                    ymax = data.plot.scenarios %>% pull(value) %>% max()),
                width = 0.1,
                linewidth = 0.25) +
  # Upscaling backed by subsidies
  new_scale_color() +
  geom_ribbon(data = data.backed,
              mapping = aes(x = year, ymin = Conservative, ymax = Progressive, fill = "Upscaling"),
              alpha = 0.2) +
  geom_line(data = data.backed,
            mapping = aes(x = year, y = Default, color = "Upscaling")) +
  scale_fill_manual(name = NULL,
                    values = c("Upscaling" = "grey30"),
                    labels = c("Upscaling" = "Supported by implemented demand policies\nand announced subsidies"),
                    guide = guide_legend(order = 3)) +
  scale_color_manual(name = NULL,
                     values = c("Upscaling" = "grey30"),
                     labels = c("Upscaling" = "Supported by implemented demand policies\nand announced subsidies"),
                     guide = guide_legend(order = 3)) +
  # Horizontal dashed line (v2023)
  geom_segment(
    mapping = aes(
      x = 2030,
      xend = 2030.4,
      y = data.v2023.plot.right %>% filter(year == 2030) %>% pull(sum),
      yend = data.v2023.plot.right %>% filter(year == 2030) %>% pull(sum)),
    linetype = "dashed",
    linewidth = 0.25
    ) +
  # Horizontal dashed line (backed by subsidies)
  geom_segment(
    mapping = aes(
      x = 2030,
      xend = 2030.3,
      y = data.backed %>% filter(year == 2030) %>% pull(Default),
      yend = data.backed %>% filter(year == 2030) %>% pull(Default)
    ),
    linetype = "dashed",
    linewidth = 0.25
  ) +
  # Horizontal dashed line (scenarios)
  geom_segment(
    mapping = aes(
      x = 2030,
      xend = 2030.4,
      y = data.plot.scenarios %>% filter(year.publication == 2023) %>% pull(value),
      yend = data.plot.scenarios %>% filter(year.publication == 2023) %>% pull(value)
    ),
    linetype = "dashed",
    linewidth = 0.25
  ) +
  # Annotations: Future implementation gap
  geom_segment(
    aes(x = 2030.15,
        xend = 2030.15,
        y = data.backed %>% filter(year == 2030) %>% pull(Default),
        yend = data.v2023.plot.right %>% filter(year == 2030) %>% pull(sum)),
    arrow = arrow(length = unit(1.5, "mm"), ends = "both"),
    linewidth = 0.5
  ) +
  geom_text(
    aes(x = 2030.15,
        y = mean(c(data.backed %>% filter(year == 2030) %>% pull(Default),
                   data.v2023.plot.right %>% filter(year == 2030) %>% pull(sum))),
        label = "\u2462 2030\nimplementation\ngap"),
    size = 6 / .pt,
    hjust = 0,
    vjust = 0.5,
    nudge_x = 0.08,
    lineheight = 1,
    fontface = "bold"
  ) +
  # Annotations: Future ambition gap
  geom_segment(
    aes(x = 2030.25,
        xend = 2030.25,
        y = data.v2023.plot.right %>% filter(year == 2030) %>% pull(sum),
        yend = data.plot.scenarios %>% filter(year.publication == 2023) %>% pull(value)),
    arrow = arrow(length = unit(1.5, "mm"), ends = "both"),
    linewidth = 0.5) +
  geom_text(
    aes(x = 2030.25,
        y = mean(c(data.v2023.plot.right %>% filter(year == 2030) %>% pull(sum),
                   data.plot.scenarios %>% filter(year.publication == 2023) %>% pull(value))),
        label = "\u2461 2030\nambition\ngap"),
    size = 6 / .pt,
    hjust = 0,
    vjust = 0.5,
    nudge_x = 0.08,
    lineheight = 1,
    fontface = "bold"
  ) +
  # Annotations: IEA scenarios in 2030
  geom_text(
    aes(x = 2030.05,
        y = mean(
          c(
            data.plot.scenarios %>% pull(value) %>% min(),
            data.plot.scenarios %>% pull(value) %>% max()
          )
        )),
    label = "1.5°C\nscenario\nrequirements",
    size = 5 / .pt,
    hjust = 0,
    lineheight = 1
  ) +
  # Formatting
  xlab("Year") +
  scale_x_continuous(breaks = seq(2024,2030)) +
  ylab("Capacity (GW)") +
  ggtitle("2030 green hydrogen ambition and implementation gap") +
  theme(legend.position = "inside",
        legend.position.inside = c(0.1,0.9),
        legend.justification = c("left","top"),
        plot.margin = margin(t = 1, r = 15, b = 1, l = 1, unit = "mm")) +
  coord_cartesian(xlim = c(2024,2030), clip = "off")

p <- plot_grid(
  p.left,
  p.right,
  ncol = 2,
  labels = "auto",
  label_size = font.size,
  rel_widths = c(1,1.45))

print(p)

ggsave("figures/FIGURE01_overview.png",
       width = 180, height = 70, units = "mm", bg = "white", dpi = 600)
ggsave("figures/FIGURE01_overview.svg",
       width = 180, height = 70, units = "mm", bg = "white")

```

## FIGURE 2: Comparison of database versions

```{r figure 2, warning=FALSE}

range.left <- seq(2020,2024)

# Data for left bar (v2021)
data.v2021.plot.left <- data.v2021 %>% 
  filter(year %in% range.left) %>% 
  mutate(status = case_when(status == "Decommissioned" ~ "Operational",
                            TRUE ~ status)) %>% 
  group_by(status, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

# Data for middle bar (v2022)
data.v2022.plot.left <- data.v2022 %>% 
  filter(year %in% range.left) %>% 
  mutate(status = case_when(status == "Decommissioned" ~ "Operational",
                            TRUE ~ status)) %>% 
  group_by(status, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

# Data for right bar (v2023)
data.v2023.plot.left <- data.v2023 %>% 
  filter(year %in% range.left) %>% 
  mutate(status = case_when(status == "Decommissioned" ~ "Operational",
                            TRUE ~ status)) %>% 
  group_by(status, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

data.plot.text <- tribble(
  ~year, ~value, ~text,
  2021.75, data.v2021.plot.left %>%
    filter(year == 2022) %>%
    pull(cumcap.sum) %>%
    sum(), "Announcements\nby 2021",
  2022, data.v2022.plot.left %>%
    filter(year == 2022) %>%
    pull(cumcap.sum) %>%
    sum(), "Announcements\nby 2022",
  2022.25, data.v2023.plot.left %>%
    filter(year == 2022) %>%
    pull(cumcap.sum) %>%
    sum(), "Outcome\nin 2023"
)

# Add box around year 2024
data.plot.2024boxheight <- data.v2023.plot.left %>% 
  filter(year == 2024) %>% 
  pull(cumcap.sum) %>% 
  sum()

p.status.left <- ggplot() + 
  # Add data snapshots
  geom_bar(data = data.v2021.plot.left,
           mapping = aes(x = year, y = cumcap.sum, fill = status),
           stat = "identity",
           position = position_stacknudge(x = -0.25),
           width = 0.2) +
  geom_bar(data = data.v2022.plot.left,
           mapping = aes(x = year, y = cumcap.sum, fill = status),
           stat = "identity",
           position = position_stacknudge(x = 0),
           width = 0.2) +
  geom_bar(data = data.v2023.plot.left,
           mapping = aes(x = year, y = cumcap.sum, fill = status),
           stat = "identity",
           position = position_stacknudge(x = 0.25),
           width = 0.2) +
  # Add labels to bars in 2022
  geom_label_repel(data = data.plot.text,
                   mapping = aes(x = year, y = value, label = text),
                   size = 5 / .pt,
                   nudge_x = c(-1, -0.3, 0.6),
                   nudge_y = c(4,7,10),
                   direction = c("y"),
                   segment.curvature = 0.1,
                   segment.shape = 0.5,
                   segment.square = FALSE,
                   segment.inflect = TRUE,
                   segment.ncp = 10,
                   segment.size = 0.25) +
  # Add box around 2024
  geom_rect(mapping = aes(xmin = 2023.6, xmax = 2024.4,
                          ymin = -0.2, ymax = data.plot.2024boxheight + 0.2),
            fill = "transparent",
            color = "grey40",
            linewidth = 0.25) +
  scale_fill_npg(name = "Status") +
  scale_x_continuous(name = "Year",
                     breaks = range.left) +
  scale_y_continuous(name = "Capacity (GW)") +
  ggtitle("Projects until 2024 by status") +
  background_grid(major = c("y"),
                  size.major = 0.25)

range.right <- seq(2024,2030)

data.v2021.plot.right <- data.v2021 %>% 
  filter(year %in% range.right) %>% 
  mutate(status = case_when(status == "Decommissioned" ~ "Operational",
                            TRUE ~ status)) %>% 
  group_by(status, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

data.v2022.plot.right <- data.v2022 %>% 
  filter(year %in% range.right) %>% 
  mutate(status = case_when(status == "Decommissioned" ~ "Operational",
                            TRUE ~ status)) %>% 
  group_by(status, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

data.v2023.plot.right <- data.v2023 %>% 
  filter(year %in% range.right) %>% 
  mutate(status = case_when(status == "Decommissioned" ~ "Operational",
                            TRUE ~ status)) %>% 
  group_by(status, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

data.plot.text <- tribble(
  ~year, ~value, ~text,
  2026.75, data.v2021.plot.right %>%
    filter(year == 2027) %>%
    pull(cumcap.sum) %>%
    sum(), "Announcements\nby 2021",
  2027, data.v2022.plot.right %>%
    filter(year == 2027) %>%
    pull(cumcap.sum) %>%
    sum(), "Announcements\nby 2022",
  2027.25, data.v2023.plot.right %>%
    filter(year == 2027) %>%
    pull(cumcap.sum) %>%
    sum(), "Announcements\nby 2023"
)

p.status.right <- ggplot() + 
  geom_bar(data = data.v2021.plot.right,
           mapping = aes(x = year, y = cumcap.sum, fill = status),
           stat = "identity",
           position = position_stacknudge(x = -0.25),
           width = 0.2) +
  geom_bar(data = data.v2022.plot.right,
           mapping = aes(x = year, y = cumcap.sum, fill = status),
           stat = "identity",
           position = position_stacknudge(x = 0),
           width = 0.2) +
  geom_bar(data = data.v2023.plot.right,
           mapping = aes(x = year, y = cumcap.sum, fill = status),
           stat = "identity",
           position = position_stacknudge(x = 0.25),
           width = 0.2) +
  # Add labels to bars in 2027
  geom_label_repel(data = data.plot.text,
                   mapping = aes(x = year, y = value, label = text),
                   size = 5 / .pt,
                   nudge_x = c(-0.7, -0.2, 0.6),
                   nudge_y = c(180, 190, 200),
                   direction = c("y"),
                   segment.curvature = 0.1,
                   segment.shape = 0.5,
                   segment.square = FALSE,
                   segment.inflect = TRUE,
                   segment.ncp = 10,
                   segment.size = 0.25) +
  # Add box around 2024
  geom_rect(mapping = aes(xmin = 2023.6, xmax = 2024.4,
                          ymin = -0.2, ymax = data.plot.2024boxheight + 0.2),
            fill = "transparent",
            color = "grey40",
            linewidth = 0.25) +
  scale_fill_npg(name = "Status") +
  scale_x_continuous(name = "Year",
                     breaks = range.right) +
  scale_y_continuous(name = "Capacity (GW)") +
  ggtitle("Projects from 2024 by status") +
  background_grid(major = c("y"),
                  size.major = 0.25)

data.v2021.plot.left <- data.v2021 %>% 
  filter(year %in% range.left) %>% 
  group_by(region, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

data.v2022.plot.left <- data.v2022 %>% 
  filter(year %in% range.left) %>% 
  group_by(region, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

data.v2023.plot.left <- data.v2023 %>% 
  filter(year %in% range.left) %>% 
  group_by(region, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

p.regi.left <- ggplot() + 
  geom_bar(data = data.v2021.plot.left,
           mapping = aes(x = year, y = cumcap.sum, fill = region),
           stat = "identity",
           position = position_stacknudge(x = -0.25),
           width = 0.2) +
  geom_bar(data = data.v2022.plot.left,
           mapping = aes(x = year, y = cumcap.sum, fill = region),
           stat = "identity",
           position = position_stacknudge(x = 0),
           width = 0.2) +
  geom_bar(data = data.v2023.plot.left,
           mapping = aes(x = year, y = cumcap.sum, fill = region),
           stat = "identity",
           position = position_stacknudge(x = 0.25),
           width = 0.2) +
  geom_rect(mapping = aes(xmin = 2023.6, xmax = 2024.4,
                          ymin = -0.2, ymax = data.plot.2024boxheight + 0.2),
            fill = "transparent",
            color = "grey40",
            linewidth = 0.25) +
  scale_fill_npg(name = "Region") +
  scale_x_continuous(name = "Year",
                     breaks = range.left) +
  scale_y_continuous(name = "Capacity (GW)") +
  ggtitle("Projects until 2024 by region") +
  background_grid(major = c("y"),
                  size.major = 0.25)

data.v2021.plot.right <- data.v2021 %>% 
  filter(year %in% range.right) %>% 
  group_by(region, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

data.v2022.plot.right <- data.v2022 %>% 
  filter(year %in% range.right) %>% 
  group_by(region, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

data.v2023.plot.right <- data.v2023 %>% 
  filter(year %in% range.right) %>% 
  group_by(region, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

p.regi.right <- ggplot() + 
  geom_bar(data = data.v2021.plot.right,
           mapping = aes(x = year, y = cumcap.sum, fill = region),
           stat = "identity",
           position = position_stacknudge(x = -0.25),
           width = 0.2) +
  geom_bar(data = data.v2022.plot.right,
           mapping = aes(x = year, y = cumcap.sum, fill = region),
           stat = "identity",
           position = position_stacknudge(x = 0),
           width = 0.2) +
  geom_bar(data = data.v2023.plot.right,
           mapping = aes(x = year, y = cumcap.sum, fill = region),
           stat = "identity",
           position = position_stacknudge(x = 0.25),
           width = 0.2) +
  geom_rect(mapping = aes(xmin = 2023.6, xmax = 2024.4,
                          ymin = -0.2, ymax = data.plot.2024boxheight + 0.2),
            fill = "transparent",
            color = "grey40",
            linewidth = 0.25) +
  scale_fill_npg(name = "Region") +
  scale_x_continuous(name = "Year",
                     breaks = range.right) +
  scale_y_continuous(name = "Capacity (GW)") +
  ggtitle("Projects from 2024 by region") +
  background_grid(major = c("y"),
                  size.major = 0.25)

p.panels <- plot_grid(
  p.status.left + theme(legend.position = "none"),
  p.status.right + theme(legend.position = "none"),
  p.regi.left + theme(legend.position = "none"),
  p.regi.right + theme(legend.position = "none"),
  labels = "auto",
  rel_widths = c(5,7),
  label_size = font.size)

p.leg.top <- get_legend(p.status.left)

p.leg.bottom <- get_legend(p.regi.left)

p.leg <- plot_grid(p.leg.top, p.leg.bottom, ncol = 1)

p.tot <- plot_grid(p.panels,
                   NULL,
                   p.leg,
                   ncol = 3,
                   rel_widths = c(1, 0.01, 0.13))

print(p.tot)

ggsave("figures/FIGURE02_database_comparison.png",
       width = 180, height = 120, units = "mm", bg = "white", dpi = 600)
ggsave("figures/FIGURE02_database_comparison.svg",
       width = 180, height = 120, units = "mm", bg = "white")

```

# Past green hydrogen implementation gap

## Calculate project tracking for 2023

```{r calc tracking 2023}

data.track.calc2023 <-
  calcProjectTracking2023(
    data.v2021.proj = data.v2021.proj,
    data.v2022.proj = data.v2022.proj,
    data.v2023.proj = data.v2023.proj,
    data.v2023outcome.proj = data.v2023outcome.proj
  )

```

## Calculate rates of success, delay and disappearance for 2023

```{r calc rates 2023}

# Comparing announcements in 2021 to outcome
data.rates.calc2023.v21out <-
  calcRates(
    data.track.calc = data.track.calc2023,
    year.tracking = 2023,
    compare = "v2021-out"
  )

# Print rates
data.rates.calc2023.v21out %>% 
  group_by(outcome) %>% 
  summarise(rate.tot = round(sum(cap.share.tot), 3))

data.rates.calc2023.v21out %>% 
  group_by(status, outcome) %>% 
  summarise(rate = round(sum(cap.share.disagg), 3))

# Comparing announcements in 2022 to outcome
data.rates.calc2023.v22out <-
  calcRates(
    data.track.calc = data.track.calc2023,
    year.tracking = 2023,
    compare = "v2022-out"
  )

# Print rates
data.rates.calc2023.v22out %>% 
  group_by(outcome) %>% 
  summarise(rate.tot = round(sum(cap.share.tot), 3))

data.rates.calc2023.v22out %>% 
  group_by(status, outcome) %>% 
  summarise(rate = round(sum(cap.share.disagg), 3))

# Comparing announcements in 2023 to outcome
data.rates.calc2023.v23out <-
  calcRates(
    data.track.calc = data.track.calc2023,
    year.tracking = 2023,
    compare = "v2023-out"
  )

# Print rates
data.rates.calc2023.v23out %>% 
  group_by(outcome) %>% 
  summarise(rate.tot = round(sum(cap.share.tot), 3))

data.rates.calc2023.v23out %>% 
  group_by(status, outcome) %>% 
  summarise(rate = round(sum(cap.share.disagg), 3))

```


## FIGURE 3: 2023 implementation gap

```{r figure 3, warning = FALSE}

regions.all <- data.track.calc2023 %>% pull(region) %>% unique()

p.gap2023 <- plotImplementationGap2023(
  region.plot = regions.all,
  title = "Tracking global green hydrogen projects announced for 2023",
  shift.incoming = 8500,
  shift.outgoing = 900,
  y.upper = 5000,
  legend = "v23")

print(p.gap2023)

ggsave("figures/FIGURE03_implementation_gap_2023.png",
       width = 180, height = 150, units = "mm", bg = "white", dpi = 600)
ggsave("figures/FIGURE03_implementation_gap_2023.svg",
       width = 180, height = 150, units = "mm", bg = "white")

```

# Future green hydrogen ambition gap

## FIGURE 4: 1.5°C scenarios and ambition gap in 2030

```{r figure 4}

# Institutional scenarios
data.scenarios.2030 <- data.scenarios %>% 
  filter(year == 2030) %>% 
  # Highlight IEA scenarios
  mutate(highlight = ifelse(organisation == "IEA", "yes", "no")) %>% 
  # Remove Deloitte as outlier (1700 GW in 2030)
  filter(organisation != "Deloitte")

# IPCC AR6 C1 scenarios
data.scenarios.ar6.2030 <- data.ar6 %>% 
  filter(period == 2030) %>% 
  select(model, scenario, region, period, value)

# NGFS scenarios
data.scenarios.ngfs.2030 <- data.ngfs %>% 
  filter(period == 2030) %>% 
  select(model, scenario, region, period, value)

# Combine IAM
data.plot.iam.2030 <- bind_rows(
  data.scenarios.ar6.2030,
  data.scenarios.ngfs.2030
)

# Scenario plot
set.seed(16)
p.scen <- ggplot() +
  # IAM scenarios
  geom_boxplot(data = data.plot.iam.2030,
               mapping = aes(x = "IAM", y = value),
               linewidth = 0.25,
               outliers = FALSE,
               width = 0.15) +
  geom_sina(data = data.plot.iam.2030,
            mapping = aes(x = "IAM", y = value),
            color = "grey50",
            size = 0.5,
            alpha = 0.5) +
  geom_violin(data = data.plot.iam.2030,
              mapping = aes(x = "IAM", y = value),
              linewidth = 0.25,
              fill = NA) +
  # Institutional scenarios
  geom_boxplot(data = data.scenarios.2030,
               mapping = aes(x = "Institutional", y = value.mw),
               linewidth = 0.25,
               width = 0.15) +
  geom_sina(data = data.scenarios.2030,
              mapping = aes(x = "Institutional", y = value.mw, color = highlight, group = factor("Institutional")),
            size = 1,
            alpha = 0.5,
            show.legend = FALSE) +
  scale_color_manual(values = c("yes" = "darkred",
                                "no" = "grey50")) +
  geom_violin(data = data.scenarios.2030,
              mapping = aes(x = "Institutional", y = value.mw),
              linewidth = 0.25,
              fill = NA) +
  scale_x_discrete(name = NULL,
                   labels = c("IAM" = "IAM\nscenarios",
                              "Institutional" = "Institutional and\ncorporate scenarios")) +
  scale_y_continuous(
    name = "Capacity in 2030 (GW)",
    limits = c(0, 1080)) +
  ggtitle("Green hydrogen in 1.5°C scenarios in 2030")

# Electrolysis projects announced for 2030 (v2021)
data.v2021.proj.2030 <- data.v2021.proj %>% 
  filter(year <= 2030) %>% 
  group_by(reference) %>% 
  filter(!any(status == "Decommissioned")) %>% 
  complete(status, year = 2010:2030, name, region, fill = list(capacity = 0)) %>% 
  group_by(reference, name, status) %>% 
  arrange(year) %>% 
  mutate(cumcap.sum = 1E-3*cumsum(capacity)) %>%  # in GW
  filter(year == 2030) %>% 
  arrange(cumcap.sum) %>% 
  ungroup() %>% 
  order.levels(status = c("Concept", "Feasibility study", "FID/Construction", "Operational"))

# Electrolysis projects announced for 2030 (v2022)
data.v2022.proj.2030 <- data.v2022.proj %>% 
  filter(year <= 2030) %>% 
  group_by(reference) %>% 
  filter(!any(status == "Decommissioned")) %>% 
  complete(status, year = 2010:2030, name, region, fill = list(capacity = 0)) %>% 
  group_by(reference, name, status) %>% 
  arrange(year) %>% 
  mutate(cumcap.sum = 1E-3*cumsum(capacity)) %>%  # in GW
  filter(year == 2030) %>% 
  arrange(cumcap.sum) %>% 
  ungroup() %>% 
  order.levels(status = c("Concept", "Feasibility study", "FID/Construction", "Operational"))

# Electrolysis projects announced for 2030 (v2023)
data.v2023.proj.2030 <- data.v2023.proj %>% 
  filter(year <= 2030) %>% 
  group_by(reference) %>% 
  filter(!any(status == "Decommissioned")) %>% 
  complete(status, year = 2010:2030, name, region, fill = list(capacity = 0)) %>% 
  group_by(reference, name, status) %>% 
  arrange(year) %>% 
  mutate(cumcap.sum = 1E-3*cumsum(capacity)) %>%  # in GW
  filter(year == 2030) %>% 
  arrange(cumcap.sum) %>% 
  ungroup() %>% 
  order.levels(status = c("Concept", "Feasibility study", "FID/Construction", "Operational"))

# Distance between arrow head and value
d <- 5

# IEA NZE 2030 values
iea.2030.v2021 <- data.scenarios.2030 %>%
  filter(scen.name == "IEA NZE 2021",
         year == 2030) %>%
  pull(value.mw)

iea.2030.v2022 <- data.scenarios.2030 %>%
  filter(scen.name == "IEA NZE 2022",
         year == 2030) %>%
  pull(value.mw)

iea.2030.v2023 <- data.scenarios.2030 %>%
  filter(scen.name == "IEA NZE 2023",
         year == 2030) %>%
  pull(value.mw)

# Arrow data
data.arrows <- tribble(
  ~ iea.version, ~ from, ~ to,
  "IEA v2021", data.v2021.proj.2030 %>% pull(cumcap.sum) %>% sum() + d,
  iea.2030.v2021 - d,
  "IEA v2022", data.v2022.proj.2030 %>% pull(cumcap.sum) %>% sum() + d,
  iea.2030.v2022 - d,
  "IEA v2023", data.v2023.proj.2030 %>% pull(cumcap.sum) %>% sum() + d,
  iea.2030.v2023 - d,
)

# Number of projects
data.number <- tribble(
  ~ iea.version, ~ cap, ~ number,
  "IEA v2021", data.v2021.proj.2030 %>% pull(cumcap.sum) %>% sum(),
  data.v2021.proj.2030 %>% pull(reference) %>% unique() %>% length(),
  "IEA v2022", data.v2022.proj.2030 %>% pull(cumcap.sum) %>% sum(),
  data.v2022.proj.2030 %>% pull(reference) %>% unique() %>% length(),
  "IEA v2023", data.v2023.proj.2030 %>% pull(cumcap.sum) %>% sum(),
  data.v2023.proj.2030 %>% pull(reference) %>% unique() %>% length()
)

# Scenario range
data.plot.scenarios.range <- data.scenarios.2030 %>% 
  summarise(min = min(value.mw),
            max = max(value.mw),
            lower = quantile(value.mw, 0.25),
            upper = quantile(value.mw, 0.75),
            median = median(value.mw))

# Plot closing green hydrogen ambition gap in 2030
p.gap <- ggplot() + 
  # Plot announcements
  geom_bar(data = data.v2021.proj.2030,
           mapping = aes(x = "IEA v2021", y = cumcap.sum, fill = status, color = status),
           stat = "identity",
           alpha = 0.5,
           linewidth = 0.15) +
  geom_bar(data = data.v2022.proj.2030,
           mapping = aes(x = "IEA v2022", y = cumcap.sum, fill = status, color = status),
           stat = "identity",
           alpha = 0.5,
           linewidth = 0.15,
           show.legend = FALSE) +
  geom_bar(data = data.v2023.proj.2030,
           mapping = aes(x = "IEA v2023", y = cumcap.sum, fill = status, color = status),
           stat = "identity",
           alpha = 0.5,
           linewidth = 0.15,
           show.legend = FALSE) +
  scale_fill_npg(name = "Project status") +
  scale_color_npg(name = "Project status") +
  # Plot number of projects
  geom_text(data = data.number,
            mapping = aes(x = iea.version, y = cap, label = paste0("n=", number)),
            size = 5 / .pt,
            position = position_nudge(x = 0.43, y = 5),
            vjust = 0,
            hjust = 1) +
  # Plot IEA NZE scenarios and dashed horizontal lines
  geom_segment(aes(x = 0.55, xend = 1.45,
                   y = iea.2030.v2021, yend = iea.2030.v2021),
               linewidth = 0.25) +
  geom_segment(aes(x = 1.45, xend = 3.45,
                   y = iea.2030.v2021, yend = iea.2030.v2021),
               linetype = "dashed",
               linewidth = 0.25,
               alpha = 0.3) +
  geom_text(aes(x = 1, y = iea.2030.v2021, label = "NZE 2021"),
            size = 5 / .pt,
            nudge_y = -10,
            vjust = -1) +
  geom_segment(aes(x = 1.55, xend = 2.45,
                   y = iea.2030.v2022, yend = iea.2030.v2022),
               linewidth = 0.25) +
  geom_segment(aes(x = 2.45, xend = 3.45,
                   y = iea.2030.v2022, yend = iea.2030.v2022),
               linetype = "dashed",
               linewidth = 0.25,
               alpha = 0.3) +
  geom_text(aes(x = 2, y = iea.2030.v2022, label = "NZE 2022"),
            size = 5 / .pt,
            nudge_y = -10,
            vjust = -1) +
  geom_segment(aes(x = 2.55, xend = 3.45,
                   y = iea.2030.v2023, yend = iea.2030.v2023),
               linewidth = 0.25) +
  geom_text(aes(x = 3, y = iea.2030.v2023, label = "NZE 2023"),
            size = 5 / .pt,
            nudge_y = -10,
            vjust = -1) +
  # Plot arrows
  geom_segment(
    data = data.arrows,
    mapping = aes(
      x = iea.version,
      y = from,
      xend = iea.version,
      yend = to
    ),
    arrow = arrow(length = unit(0.03, "npc"), ends = "both"),
    linewidth = 0.25) +
  # Plot braces
  stat_brace(aes(x = c(3.5, 3.6),
                 y = c(3, data.v2023.proj.2030 %>% pull(cumcap.sum) %>% sum() - 3)),
             rotate = 90,
             linewidth = 0.25,
             outside = FALSE) +
  stat_bracetext(aes(x = c(3.5, 3.6),
                     y = c(3, data.v2023.proj.2030 %>% pull(cumcap.sum) %>% sum() - 3)),
                 rotate = 90,
                 outside = FALSE,
                 label = "Increasing\ngreen hydrogen\nproject pipeline",
                 size = 6 / .pt) +
  stat_brace(aes(x = c(3.5, 3.6),
                 y = c(data.v2023.proj.2030 %>% pull(cumcap.sum) %>% sum() + 3, iea.2030.v2023 - 3)),
             rotate = 90,
             outside = FALSE,
             linewidth = 0.25) +
  stat_bracetext(aes(x = c(3.5, 3.6),
                     y = c(data.v2023.proj.2030 %>% pull(cumcap.sum) %>% sum() + 3, iea.2030.v2023 - 3)),
                 rotate = 90,
                 outside = FALSE,
                 label = "\u2461 Green hydrogen\nambition gap in 2030",
                 size = 6 /.pt,
                 fontface = "bold") +
  stat_brace(aes(x = c(3.5, 3.6),
                 y = c(iea.2030.v2023 + 3, iea.2030.v2021 - 3)),
             rotate = 90,
             outside = FALSE,
             linewidth = 0.25) +
  stat_bracetext(aes(x = c(3.5, 3.6),
                     y = c(iea.2030.v2023 + 3, iea.2030.v2021 - 3)),
                 rotate = 90,
                 outside = FALSE,
                 label = "Decreasing role of\ngreen hydrogen in\nIEA Net-Zero Emissions\nscenarios",
                 size = 6 /.pt) +
  coord_cartesian(clip = "off") +
  scale_x_discrete(name = "Year of announcement / publication",
                   labels = c("IEA v2021" = "2021",
                              "IEA v2022" = "2022",
                              "IEA v2023" = "2023")) +
  scale_y_continuous(name = "Capacity in 2030 (GW)",
                     limits = c(0, 1080)) +
  ggtitle("Green hydrogen ambition gap in 2030") +
  theme(plot.margin = unit(c(3,50,3,3), units = "pt"),
        axis.title.y = element_text(hjust = 0.25))

# Arrange plot panels
p.panels <- plot_grid(
  p.scen,
  NULL,
  p.gap + theme(legend.position = "none"),
  align = "h",
  labels = c("a", "", "b"),
  label_size = font.size,
  ncol = 3,
  rel_widths = c(0.5, 0.05, 1)
)

# Combine with legend to full plot
p.all <- plot_grid(
  p.panels,
  NULL,
  get_legend(p.gap),
  ncol = 3,
  rel_widths = c(1,0.07,0.15)
)

print(p.all)

ggsave("figures/FIGURE04_ambition_gap_2030.png",
       width = 160, height = 70, units = "mm", bg = "white", dpi = 600)
ggsave("figures/FIGURE04_ambition_gap_2030.svg",
       width = 160, height = 70, units = "mm", bg = "white")

```

# Future green hydrogen implementation gap

## Read IEA Database with end use

```{r read iea with enduse}

file.v2023.checked <-
  "data/IEA_H2_DB_2023-10/quality_checked/IEA_Hydrogen_Projects_Database_2023_quality_checked.xlsx"

temp <- readIEAHydrogenProjectsDBWithEndUse(
  file.h2projects = file.v2023.checked,
  h2projects.range = "A5:AG2005",
  file.h2projects.2023outcome = file.v2023outcome,
  h2projects.range.2023outcome = "A5:AG128",
  end.year = 2023,
  checked = TRUE
)

data.v2023.proj.enduse <- temp[[1]]
data.v2023.enduse <- temp[[2]]

```

## Read mapping of enduse to green-fossil competition and transport costs

```{r read mapping enduse}

data.enduse.transportcost <- read_excel("data/enduse_green_fossil_transportcost.xlsx") %>% 
  mutate(enduse = ifelse(enduse == "NA", NA_character_, enduse))

data.iea.enduse <- data.v2023.enduse %>% 
  full_join(data.enduse.transportcost)

data.iea.proj.enduse <- data.v2023.proj.enduse %>% 
  full_join(data.enduse.transportcost)

enduse.labels <- c(
  "ammonia" = "Ammonia",
  "biofuels" = "Biofuels",
  "ch4gridinj" = "CH4 grid injection",
  "ch4mobility" = "CH4 mobility",
  "chp" = "CHP",
  "domesticheat" = "Domestic heat",
  "gridinj" = "Grid injection",
  "ironsteel" = "Iron & Steel",
  "methanol" = "Methanol",
  "mobility" = "Mobility",
  "otherind" = "Other industry",
  "power" = "Power",
  "refining" = "Refining",
  "synfuels" = "Synthetic fuels")

```

## Read cost data and calculate technological learning

```{r read costs and calculate learning}

file.costs <- "data/cost_scenarios_enduse.xlsx"

data.costs <- readTechCalcLearning(
  file.costs = file.costs,
  data.iea.enduse = data.iea.enduse,
  data.scenarios = data.scenarios
)

```

## Calculate levelised cost of green hydrogen (LCOH) and LCOX of green hydrogen derivatives

```{r calculate lcoh}

data.lcox <- data.costs %>% 
  # Calculate annuities
  calc_addVariable(
    "`Annuity`" = "`Tech|Interest` / (1 - (1 + `Tech|Interest`)**(-`Tech|Payback period`))",
    "`Annuity|Stack`" = "`Tech|Interest` / (1 - (1 + `Tech|Interest`)**(-`Tech|Green hydrogen|Lifetime|Stack`))",
    units = "1") %>% 
  calc_addVariable(
    # GREEN HYDROGEN
    "`LCOX|Green hydrogen|CAPEX|Other`" = "( 1 / `Tech|Green hydrogen|Efficiency` ) * `Annuity` * `Tech|Green hydrogen|CAPEX|Other` / `Tech|Green hydrogen|FLH`",
    "`LCOX|Green hydrogen|CAPEX|Stack`" = "( 1 / `Tech|Green hydrogen|Efficiency` ) * `Annuity|Stack` * `Tech|Green hydrogen|CAPEX|Stack` / `Tech|Green hydrogen|FLH`",
    "`LCOX|Green hydrogen|FOM`" = "( 1 / `Tech|Green hydrogen|Efficiency` ) * `Tech|Green hydrogen|FOM` * `Tech|Green hydrogen|CAPEX` / `Tech|Green hydrogen|FLH`",
    "`LCOX|Green hydrogen|Electricity`" =  "( 1 / `Tech|Green hydrogen|Efficiency` ) * `Price|Electricity|Hydrogen`",
    "`LCOX|Green hydrogen`" =  "`LCOX|Green hydrogen|CAPEX|Other` + `LCOX|Green hydrogen|CAPEX|Stack` + `LCOX|Green hydrogen|FOM` + `LCOX|Green hydrogen|Electricity`",
    # E-METHANOL
    "`LCOX|E-Methanol|CAPEX`" = "`Annuity` * `Tech|E-Methanol|CAPEX` / `Tech|E-Methanol|FLH`",
    "`LCOX|E-Methanol|FOM`" = "`Tech|E-Methanol|FOM` * `Tech|E-Methanol|CAPEX` / `Tech|E-Methanol|FLH`",
    "`LCOX|E-Methanol|Hydrogen`" = "( 1 / `Tech|E-Methanol|Efficiency` ) * `LCOX|Green hydrogen`",
    "`LCOX|E-Methanol|Carbon`" = "`Price|CO2|Renewable` * `Tech|E-Methanol|CO2 intensity`",
    # E-METHANE
    "`LCOX|E-Methane|CAPEX`" = "`Annuity` * `Tech|E-Methane|CAPEX` / `Tech|E-Methane|FLH`",
    "`LCOX|E-Methane|FOM`" = "`Tech|E-Methane|FOM` * `Tech|E-Methane|CAPEX` / `Tech|E-Methane|FLH`",
    "`LCOX|E-Methane|Hydrogen`" = "( 1 / `Tech|E-Methane|Efficiency` ) * `LCOX|Green hydrogen`",
    "`LCOX|E-Methane|Carbon`" = "`Price|CO2|Renewable` * `Tech|E-Methane|CO2 intensity`",
    # E-KEROSENE
    "`LCOX|E-Kerosene|CAPEX`" = "`Annuity` * `Tech|E-Kerosene|CAPEX` / `Tech|E-Kerosene|FLH`",
    "`LCOX|E-Kerosene|FOM`" = "`Tech|E-Kerosene|FOM` * `Tech|E-Kerosene|CAPEX` / `Tech|E-Kerosene|FLH`",
    "`LCOX|E-Kerosene|Hydrogen`" = "( 1 / `Tech|E-Kerosene|Efficiency` ) * `LCOX|Green hydrogen`",
    "`LCOX|E-Kerosene|Carbon`" = "`Price|CO2|Renewable` * `Tech|E-Kerosene|CO2 intensity`",
    units = "USD/MWh",
    only.new = TRUE
  )

# Add transport and storage costs, depending on the enduse
data.lcox.vom <- data.enduse.transportcost %>% 
  mutate(variable = paste0("LCOX|", greenproduct, "|VOM")) %>% 
  rename(value = transportcost.green) %>% 
  crossing(scenario = c("Default", "Progressive", "Conservative"),
           period = seq(2024, 2050)) %>% 
  mutate(region = "World",
         unit = "USD/MWh") %>% 
  select(scenario, variable, region, unit, period, value, greenproduct, enduse)

# Combine with VOM and process data
data.lcox <- data.lcox %>% 
  rowwise() %>% 
  mutate(greenproduct = str_split_1(variable, "\\|")[2]) %>% 
  full_join(data.enduse.transportcost %>% select(enduse, greenproduct),
            relationship = "many-to-many") %>% 
  bind_rows(data.lcox.vom) %>%
  # Calculate total costs including transport and storage costs
  filter(variable != "LCOX|Green hydrogen") %>%   # Replace again
  calc_addVariable(
    "`LCOX|Green hydrogen`" = "`LCOX|Green hydrogen|CAPEX|Other` + `LCOX|Green hydrogen|CAPEX|Stack` + `LCOX|Green hydrogen|FOM` + `LCOX|Green hydrogen|Electricity` + `LCOX|Green hydrogen|VOM`",
    "`LCOX|E-Methanol`" = "`LCOX|E-Methanol|CAPEX` + `LCOX|E-Methanol|FOM` + `LCOX|E-Methanol|Hydrogen` + `LCOX|E-Methanol|Carbon` + `LCOX|E-Methanol|VOM`",
    "`LCOX|E-Methane`" = "`LCOX|E-Methane|CAPEX` + `LCOX|E-Methane|FOM` + `LCOX|E-Methane|Hydrogen` + `LCOX|E-Methane|Carbon` + `LCOX|E-Methane|VOM`",
    "`LCOX|E-Kerosene`" = "`LCOX|E-Kerosene|CAPEX` + `LCOX|E-Kerosene|FOM` + `LCOX|E-Kerosene|Hydrogen` + `LCOX|E-Kerosene|Carbon` + `LCOX|E-Kerosene|VOM`",
    units = "USD/MWh"
  ) %>% 
  mutate(component = paste(str_split_1(variable, "\\|")[-c(1,2)], collapse = "|")) %>% 
  mutate(component = ifelse(component == "", NA_character_, component))

```

## Calculate cost components of fossil competitors

```{r calculate fossil competitors}

# Calculate fossil cost components
data.fossilcosts <- data.costs %>%
  calc_addVariable(
    # NATURAL GAS
    # Price
    "`Cost|Natural gas|w/ CO2|Price`" = "`Price|Natural gas`",
    "`Cost|Natural gas|w/o CO2|Price`" = "`Price|Natural gas`",
    # Cost of CO2 price
    "`Cost|Natural gas|w/ CO2|CO2`" = "`Price|CO2` * `Tech|Natural gas|Emission intensity`",
    "`Cost|Natural gas|w/o CO2|CO2`" = "`Price|CO2` * 0",  # Hack to get variable
    # GREY HYDROGEN
    # Price
    "`Cost|Grey hydrogen|w/ CO2|Price`" = "`Price|Grey hydrogen`",
    "`Cost|Grey hydrogen|w/o CO2|Price`" = "`Price|Grey hydrogen`",
    # Cost of CO2 price
    "`Cost|Grey hydrogen|w/ CO2|CO2`" = "`Price|CO2` * `Tech|Grey hydrogen|Emission intensity`",
    "`Cost|Grey hydrogen|w/o CO2|CO2`" = "`Price|CO2` * 0",  # Hack to get variable
    # GREY METHANOL
    # Price
    "`Cost|Grey methanol|w/ CO2|Price`" = "`Price|Grey methanol`",
    "`Cost|Grey methanol|w/o CO2|Price`" = "`Price|Grey methanol`",
    # Cost of CO2 price
    "`Cost|Grey methanol|w/ CO2|CO2`" = "`Price|CO2` * `Tech|Grey methanol|Emission intensity`",
    "`Cost|Grey methanol|w/o CO2|CO2`" = "`Price|CO2` * 0",  # Hack to get variable
    # KEROSENE
    # Price
    "`Cost|Kerosene|w/ CO2|Price`" = "`Price|Kerosene`",
    "`Cost|Kerosene|w/o CO2|Price`" = "`Price|Kerosene`",
    # Cost of CO2 price
    "`Cost|Kerosene|w/ CO2|CO2`" = "`Price|CO2` * `Tech|Kerosene|Emission intensity`",
    "`Cost|Kerosene|w/o CO2|CO2`" = "`Price|CO2` * 0",  # Hack to get variable
    # Price
    "`Cost|Diesel|w/ CO2|Price`" = "`Price|Diesel`",
    "`Cost|Diesel|w/o CO2|Price`" = "`Price|Diesel`",
    # Cost of CO2 price
    "`Cost|Diesel|w/ CO2|CO2`" = "`Price|CO2` * `Tech|Diesel|Emission intensity`",
    "`Cost|Diesel|w/o CO2|CO2`" = "`Price|CO2` * 0",  # Hack to get variable
    units = "USD/MWh",
    only.new = TRUE)

# Add transport and storage costs, depending on the enduse
data.fossilcosts.vom <- data.enduse.transportcost %>% 
  select(enduse, fossilcompetitor, transportcost.fossil) %>% 
  mutate(variable = paste0("Cost|", fossilcompetitor, "|w/ CO2|VOM")) %>%
  rename(value = transportcost.fossil) %>% 
  crossing(scenario = c("Default", "Progressive", "Conservative"),
           period = seq(2024, 2050)) %>% 
  mutate(region = "World",
         unit = "USD/MWh") %>% 
  select(scenario, variable, region, unit, period, value, fossilcompetitor, enduse) %>% 
  # Copy VOM for w/o CO2 case
  calc_addVariable(
    "`Cost|Natural gas|w/o CO2|VOM`" = "`Cost|Natural gas|w/ CO2|VOM`",
    "`Cost|Grey hydrogen|w/o CO2|VOM`" = "`Cost|Grey hydrogen|w/ CO2|VOM`",
    "`Cost|Grey methanol|w/o CO2|VOM`" = "`Cost|Grey methanol|w/ CO2|VOM`",
    "`Cost|Kerosene|w/o CO2|VOM`" = "`Cost|Kerosene|w/ CO2|VOM`",
    "`Cost|Diesel|w/o CO2|VOM`" = "`Cost|Diesel|w/ CO2|VOM`",
    units = "USD/MWh"
  )

# Combine with VOM and process data
data.fossilcosts <- data.fossilcosts %>%
  rowwise() %>% 
  mutate(fossilcompetitor = str_split_1(variable, "\\|")[2]) %>% 
  full_join(data.enduse.transportcost %>% select(enduse, fossilcompetitor),
            relationship = "many-to-many") %>% 
  bind_rows(data.fossilcosts.vom) %>% 
  # Calculate total prices including transport and storage costs
  calc_addVariable(
    # Natural gas
    "`Cost|Natural gas|w/ CO2`" = "`Cost|Natural gas|w/ CO2|Price` + `Cost|Natural gas|w/ CO2|VOM` + `Cost|Natural gas|w/ CO2|CO2`",
    "`Cost|Natural gas|w/o CO2`" = "`Cost|Natural gas|w/o CO2|Price` + `Cost|Natural gas|w/o CO2|VOM`",
    # Grey hydrogen
    "`Cost|Grey hydrogen|w/ CO2`" = "`Cost|Grey hydrogen|w/ CO2|Price` + `Cost|Grey hydrogen|w/ CO2|VOM` + `Cost|Grey hydrogen|w/ CO2|CO2`",
    "`Cost|Grey hydrogen|w/o CO2`" = "`Cost|Grey hydrogen|w/o CO2|Price` + `Cost|Grey hydrogen|w/o CO2|VOM`",
    # Grey methanol
    "`Cost|Grey methanol|w/ CO2`" = "`Cost|Grey methanol|w/ CO2|Price` +  `Cost|Grey methanol|w/ CO2|VOM` + `Cost|Grey methanol|w/ CO2|CO2`",
    "`Cost|Grey methanol|w/o CO2`" = "`Cost|Grey methanol|w/o CO2|Price` + `Cost|Grey methanol|w/o CO2|VOM`",
    # Kerosene
    "`Cost|Kerosene|w/ CO2`" = "`Cost|Kerosene|w/ CO2|Price` + `Cost|Kerosene|w/ CO2|VOM` + `Cost|Kerosene|w/ CO2|CO2`",
    "`Cost|Kerosene|w/o CO2`" = "`Cost|Kerosene|w/o CO2|Price` + `Cost|Kerosene|w/o CO2|VOM`",
    # Diesel
    "`Cost|Diesel|w/ CO2`" = "`Cost|Diesel|w/ CO2|Price` + `Cost|Diesel|w/ CO2|VOM` + `Cost|Diesel|w/ CO2|CO2`",
    "`Cost|Diesel|w/o CO2`" = "`Cost|Diesel|w/o CO2|Price` + `Cost|Diesel|w/o CO2|VOM`",
    units = "USD/MWh") %>% 
  # Make carbon price an additional column
  mutate(component = str_split_1(variable, "\\|")[4],
         carbonprice = str_split_1(variable, "\\|")[3],
         carbonprice = ifelse(carbonprice == "w/ CO2", TRUE, FALSE))

```

## Read data for demand secured by implemented policies

```{r read demand data}

data.demand <- read_excel("data/iea_hydrogen_demand_regulation.xlsx") %>% 
  pivot_longer(`2020`:`2050`, names_to = "year", values_to = "value") %>% 
  crossing(scenario = c("Default", "Progressive", "Conservative")) %>%  
  mutate(period = as.numeric(year)) %>% 
  select(scenario, variable, region, unit, period, value) %>% 
  filter(!is.na(value))

# Calculate equiv. capacity that is backed by implemented demand-side policies
data.demand.cap <- data.costs %>% 
  filter(variable %in% c("Tech|Green hydrogen|Efficiency", "Tech|Green hydrogen|FLH")) %>% 
  bind_rows(data.demand) %>%
  ungroup() %>% 
  # From MtH2/yr to GW electrolysis capacity
  calc_addVariable(
  "`Capacity|Hydrogen|Demand policies`" = "1E3 * (`Final Energy|Hydrogen|Demand policies` * 33.33) / (`Tech|Green hydrogen|Efficiency` * `Tech|Green hydrogen|FLH`)",
  "`Capacity|Industry|Other|Hydrogen|Demand policies`" = "1E3 * (`Final Energy|Industry|Other|Hydrogen|Demand policies` * 33.33) / (`Tech|Green hydrogen|Efficiency` * `Tech|Green hydrogen|FLH`)",
  "`Capacity|Industry|Refining|Hydrogen|Demand policies`" = "1E3 * (`Final Energy|Industry|Refining|Hydrogen|Demand policies` * 33.33) / (`Tech|Green hydrogen|Efficiency` * `Tech|Green hydrogen|FLH`)",
  "`Capacity|Power|Hydrogen|Demand policies`" = "1E3 * (`Final Energy|Power|Hydrogen|Demand policies` * 33.33) / (`Tech|Green hydrogen|Efficiency` * `Tech|Green hydrogen|FLH`)",
  "`Capacity|Transportation|Hydrogen|Demand policies`" = "1E3 * (`Final Energy|Transportation|Hydrogen|Demand policies` * 33.33) / (`Tech|Green hydrogen|Efficiency` * `Tech|Green hydrogen|FLH`)",
  units = "GW",
  only.new = TRUE)

# Mapping between demand policy variables and enduse sectors
demandpolicy.enduse.mapping <- tribble(
  ~enduse, ~variable,
  # Other industry
  "otherind", "Capacity|Industry|Other|Hydrogen|Demand policies",
  # Refining
  "refining", "Capacity|Industry|Refining|Hydrogen|Demand policies",
  # Power
  "power", "Capacity|Power|Hydrogen|Demand policies",
  # Synfuels to transport
  # (demand policies primarily target e-fuels for aviation/maritime)
  "synfuels", "Capacity|Transportation|Hydrogen|Demand policies",
)

# Distribute cumulative capacity in 2030 backed by demand to 2024-2030
# by keeping the share of project announcements backed by demand constant
data.addedcap <- data.iea.enduse %>%
  filter(year %in% seq(2024, 2030)) %>%
  group_by(year, enduse) %>%
  summarise(cap = sum(cap.sum))

# Distribute secured demand by sector to 2024-2030 by keeping the share of
# project announcements backed by demand policies constant across time
data.cap.backed <- data.demand.cap %>% 
  full_join(demandpolicy.enduse.mapping) %>% 
  left_join(data.addedcap, relationship = "many-to-many") %>%
  filter(!is.na(enduse)) %>% 
  group_by(scenario, variable, enduse) %>%
  # Distribute backed capacity along years 
  mutate(cap.backed = value * cap/sum(cap)) %>%
  # Fill NA enduses with zero
  group_by(scenario, year, region, unit) %>% 
  complete(enduse = data.enduse.transportcost$enduse, fill = list(cap.backed = 0)) %>%
  group_by(scenario, enduse) %>% 
  complete(year = seq(2024, 2050), fill = list(cap.backed = 0)) %>% 
  # Clean up
  ungroup() %>% 
  mutate(cap.backed = 1E3*cap.backed) %>%   # GW to MW
  select(scenario, enduse, year, cap.backed)

# cap.backed can be larger than cap if the project pipeline is not lage enough
# This is considered later

```

## Read hydrogen end use in IEA NZE (2030-2050)

```{r read iea nze end use}

data.iea.nze <- read_excel("data/iea_nze_hydrogen_use.xlsx") %>% 
  pivot_longer(`2020`:`2050`, names_to = "year", values_to = "value") %>% 
  select(-comment) %>% 
  mutate(year = as.numeric(year))

# Mapping of IEA NZE scenario end uses to IEA Hydrogen Projects Database end uses
mapping.nze.h2projects <- tribble(
  ~ variable, ~ enduse,
  "Final Energy|Industry|Refining|Hydrogen", "refining",
  # Assigning chemicals to methanol
  "Final Energy|Industry|Chemicals|Hydrogen", "methanol",
  "Final Energy|Industry|Steel|Hydrogen", "ironsteel",
  "Final Energy|Industry|Other|Hydrogen", "otherind",
  "Final Energy|Transportation|Aviation and Marine|Hydrogen", "synfuels",
  "Secondary Energy Input|Hydrogen|Electricity", "power",
  "Final Energy|Transportation|Road|Hydrogen", "mobility",
  "Final Energy|Other|Hydrogen", NA_character_
)

```

## Calculate required subsidies

```{r calculate required subsidies}

# New capacity from project announcements until 2030
data.newcap.until2030 <- data.iea.enduse %>%
  filter(year %in% seq(2024, 2030),
         # Only calculate subsidies for projects without an FID
         status %in% c("Feasibility study", "Concept")) %>% 
  group_by(year, enduse) %>% 
  summarise(cap = sum(cap.sum))

# After 2030, use median
data.scenarios.median.after2030 <- data.scenarios %>%
  filter(year %in% c(2040, 2050)) %>% 
  group_by(year) %>% 
  summarise(cumcap = median(value.mw))

# Cumulative capacity of project announcements in 2030
data.cumcap.2030.total <- data.iea.enduse %>% 
  filter(year == 2030) %>% 
  summarise(sum(cumcap.sum)) %>% 
  pull()

# From 2030 use median of 1.5°C scenario
data.cumcap.1.5Cmedian <- tibble(year = 2030, cumcap = data.cumcap.2030.total) %>%
  bind_rows(data.scenarios.median.after2030) %>%
  select(year, cumcap) %>%
  complete(year = seq(2030, 2050)) %>%
  mutate(cumcap = zoo::na.approx(cumcap))

# Create sectoral newcap data
data.newcap.from2030 <- full_join(data.iea.nze, mapping.nze.h2projects) %>% 
  # Linear interpolation
  filter(year %in% seq(2030, 2050)) %>% 
  group_by(source, report, region, variable, unit, enduse) %>% 
  complete(year = seq(2030, 2050)) %>% 
  mutate(value = na.approx(value)) %>% 
  # Calculate share of enduse
  group_by(year) %>% 
  mutate(share = value / sum(value)) %>% 
  # Calculate hypothetical capacities if the end-use of project announcements
  # would perfectly match the IEA NZE scenario
  full_join(data.cumcap.1.5Cmedian) %>% 
  mutate(value.nze = share * cumcap) %>% 
  # Calculate year-on-year change
  group_by(enduse) %>% 
  mutate(cap = value.nze - lag(value.nze)) %>% 
  # Clean up
  select(enduse, year, cap) %>% 
  filter(!is.na(cap))

# Combine added capacities into one dataframe
data.newcap <- bind_rows(data.newcap.until2030, data.newcap.from2030) %>%
  group_by(enduse) %>% 
  # Fill cap with zeros if enduse doesn't exist in IEA NZE shares
  complete(year = seq(2024, 2050), fill = list(cap = 0)) %>% 
  # Transform from GW to MW
  mutate(cap = 1E3 * cap)

# Get payback period
data.payback <- data.costs %>% 
  filter(variable == "Tech|Payback period") %>% 
  rename(payback = value,
         year = period) %>% 
  select(scenario, year, payback)

# Get full load hours of electrolysis
data.flh <- data.costs %>% 
  filter(variable == "Tech|Green hydrogen|FLH") %>% 
  select(scenario, period, value) %>% 
  rename(year = period,
         flh.electrolysis = value)

# Get efficiencies
data.eta <- data.costs %>% 
  filter(variable %in% c("Tech|Green hydrogen|Efficiency",
                         "Tech|E-Methanol|Efficiency",
                         "Tech|E-Kerosene|Efficiency",
                         "Tech|E-Methane|Efficiency")) %>% 
  # For E-Fuels get multiply electrolysis efficiency with synthesis efficiency 
  calc_addVariable(
    "`Tech|E-Methanol|Efficiency`" = "`Tech|Green hydrogen|Efficiency` * `Tech|E-Methanol|Efficiency`",
    "`Tech|E-Kerosene|Efficiency`" = "`Tech|Green hydrogen|Efficiency` * `Tech|E-Kerosene|Efficiency`",
    "`Tech|E-Methane|Efficiency`" = "`Tech|Green hydrogen|Efficiency` * `Tech|E-Methane|Efficiency`",
    units = "p.u."
  ) %>% 
  # Get energy carrier
  rowwise() %>% 
  mutate(greenproduct = str_split_1(variable, "\\|")[2]) %>% 
  select(scenario, greenproduct, period, value) %>% 
  rename(year = period,
         eta = value)

# Get LCOX
data.lcox.tot <- data.lcox %>% 
  # Get total, not individual components
  filter(is.na(component)) %>% 
  select(scenario, enduse, greenproduct, period, value) %>% 
  rename(year = period,
         lcox = value)

# Get cost of fossil competitor
data.fossil.tot <- data.fossilcosts %>% 
  # Get total, not individual components
  filter(is.na(component)) %>% 
  select(scenario, enduse, fossilcompetitor, carbonprice, period, value) %>% 
  rename(year = period,
         fossilcosts = value)


# Start with fossil costs, these are the only where future changes are taken into account
# For other LCOX-related parameters, the LCOX must be recovered over the payback period,
# so the cost gap only reduces for new projects, not for projects built in the past
data.subsidies.calc <- data.fossil.tot %>% 
  # Create additional column "year.cap" to explicitly track subsidies of projects over time
  group_by(scenario, fossilcompetitor, carbonprice, enduse, year, fossilcosts) %>%
  mutate(year.cap = year) %>% 
  complete(year.cap = seq(2024, 2050)) %>% 
  filter(year.cap <= year) %>% 
  # Respect payback period for year.cap column
  full_join(data.payback, by = c("scenario", "year")) %>% 
  filter(year - year.cap < payback) %>% 
  select(-payback) %>% 
  # Join with mapping
  full_join(data.enduse.transportcost %>% select(enduse, greenproduct, fossilcompetitor, relefficiency)) %>% 
  # Join with new capacity
  full_join(data.newcap, by = c("enduse", "year.cap" = "year")) %>% 
  # Join with capacity backed by demand policies
  full_join(data.cap.backed, by = c("scenario", "enduse", "year.cap" = "year")) %>%
  # Join with full load hours
  full_join(data.flh, by = c("scenario", "year.cap" = "year")) %>% 
  # Join with efficiencies
  full_join(data.eta, by = c("scenario", "greenproduct", "year.cap" = "year")) %>% 
  # Join with LCOX
  full_join(data.lcox.tot, by = c("scenario", "enduse", "greenproduct", "year.cap" = "year"))

# Calculate required annual subsidies
data.subsidies.annual <- data.subsidies.calc %>% 
  rowwise() %>% 
  mutate(subsidies.annual = 1E-9 * max(0, cap - cap.backed) * flh.electrolysis * eta * max(0, (lcox/relefficiency) - fossilcosts)) %>% 
  select(scenario, carbonprice, enduse, greenproduct, fossilcompetitor, year, year.cap, cap, cap.backed, flh.electrolysis, eta, relefficiency, lcox, fossilcosts, subsidies.annual) %>% 
  ungroup()

```

## Read tracked subsidies

```{r read tracked subsidies}

# BNEF data, see file
file.subsidies.tracker <- "data/subsidies_bnef.xlsx"

data.subsidies.tracker <- read_excel(file.subsidies.tracker) %>% 
  select(-c("comments", "link"))

```

## FIGURE 5: 2030 implementation gap

```{r figure 5}

p.greenh2.natgas.woco2 <- plotCostGap(
  scenario.select = "Default",
  carbonprice.select = FALSE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Natural gas",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Natural gas / Grey hydrogen",
  title = "Green hydrogen vs. Natural gas"
)

p.greenh2.natgas.wco2 <- plotCostGap(
  scenario.select = "Default",
  carbonprice.select = TRUE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Natural gas",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Natural gas / Grey hydrogen",
  title = "Green hydrogen vs. Natural gas"
)

p.greenh2.greyh2.woco2 <- plotCostGap(
  scenario.select = "Default",
  carbonprice.select = FALSE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Grey hydrogen",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Natural gas / Grey hydrogen",
  title = "Green hydrogen vs. Grey hydrogen"
)

p.greenh2.greyh2.wco2 <- plotCostGap(
  scenario.select = "Default",
  carbonprice.select = TRUE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Grey hydrogen",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Natural gas / Grey hydrogen",
  title = "Green hydrogen vs. Grey hydrogen"
)

# Plot required subsidies
temp <- plotRequiredSubsidies(
  scenario.select = "Default",
  carbonprice.select = FALSE,
  title = "Required subsidies for all projects by 2030")

p.subsidies.woco2 <- temp[[1]]
secaxis.factor.woco2 <- temp[[2]]
subsidies.annual.max.woco2 <- temp[[3]]
subsidies.cumulative.2030.woco2 <- temp[[4]]

temp <- plotRequiredSubsidies(
  scenario.select = "Default",
  carbonprice.select = TRUE,
  title = "Required subsidies for all projects by 2030",
  secaxis.factor = secaxis.factor.woco2,
  ymax.annual = subsidies.annual.max.woco2)

p.subsidies.wco2 <- temp[[1]]
secaxis.factor.wco2 <- temp[[2]]
subsidies.annual.max.wco2 <- temp[[3]]
subsidies.cumulative.2030.wco2 <- temp[[4]]

# Compare required cumulative subsidies to announced subsidies
p.subsidygap <- plotSubsidyGap(
  scenario.select = "Default",
  subsidies.cumulative.2030.woco2 = subsidies.cumulative.2030.woco2,
  subsidies.cumulative.2030.wco2 = subsidies.cumulative.2030.wco2,
  ymax = subsidies.cumulative.2030.woco2,
  title = "Required cumulative subsidies\nfor all projects by 2030"
  )

# Cumulative subsidies for conservative scenario
temp <- plotRequiredSubsidies(
  scenario.select = "Conservative",
  carbonprice.select = FALSE,
  title = "Required subsidies for all projects by 2030\nwith carbon price",
  secaxis.factor = secaxis.factor.woco2,
  ymax.annual = subsidies.annual.max.woco2)

subsidies.cumulative.2030.conservative.woco2 <- temp[[4]]

temp <- plotRequiredSubsidies(
  scenario.select = "Conservative",
  carbonprice.select = TRUE,
  title = "Required subsidies for all projects by 2030\nwith carbon price",
  secaxis.factor = secaxis.factor.woco2,
  ymax.annual = subsidies.annual.max.woco2)

subsidies.cumulative.2030.conservative.wco2 <- temp[[4]]

# Cumulative subsidies for progressive scenario
temp <- plotRequiredSubsidies(
  scenario.select = "Progressive",
  carbonprice.select = FALSE,
  title = "Required subsidies for all projects by 2030\nwith carbon price",
  secaxis.factor = secaxis.factor.woco2,
  ymax.annual = subsidies.annual.max.woco2)

subsidies.cumulative.2030.progressive.woco2 <- temp[[4]]

temp <- plotRequiredSubsidies(
  scenario.select = "Progressive",
  carbonprice.select = TRUE,
  title = "Required subsidies for all projects by 2030\nwith carbon price",
  secaxis.factor = secaxis.factor.woco2,
  ymax.annual = subsidies.annual.max.woco2)

subsidies.cumulative.2030.progressive.wco2 <- temp[[4]]

data.subsidy.range <- tribble(
  ~x, ~central, ~min, ~max,
  "Without carbon price", subsidies.cumulative.2030.woco2, subsidies.cumulative.2030.progressive.woco2, subsidies.cumulative.2030.conservative.woco2,
  "With carbon price", subsidies.cumulative.2030.wco2, subsidies.cumulative.2030.progressive.wco2, subsidies.cumulative.2030.conservative.wco2,
) %>% 
  mutate(central.round = sprintf("%.1f", round(central/1E3, 1)),
         min.round = sprintf("%.1f", round(min/1E3, 1)),
         max.round = sprintf("%.1f", round(max/1E3, 1))) %>% 
  mutate(label = paste0("$", central.round, " trillion\n($", min.round, "-", max.round, " trillion)"))

# Add labels
p.subsidygap <- p.subsidygap +
  geom_text(data = data.subsidy.range,
            mapping = aes(x = x, y = central, label = label),
            size = 5 / .pt,
            lineheight = 1,
            hjust = 0.5,
            vjust = 0,
            nudge_y = 30,
            color = c("#0e6684", "#7C4569")) +
  scale_y_continuous(name = "Cumulative subsidies (Billion $)",
                     limits = c(0, subsidies.cumulative.2030.woco2 + 100))

# Recreate panels e-f
temp <- plotRequiredSubsidies(
  scenario.select = "Default",
  carbonprice.select = FALSE,
  title = "Required subsidies for all projects by 2030",
  secaxis.factor = secaxis.factor.woco2,
  ymax.annual = (subsidies.cumulative.2030.woco2 + 100) / secaxis.factor.woco2)

p.subsidies.woco2 <- temp[[1]]

temp <- plotRequiredSubsidies(
  scenario.select = "Default",
  carbonprice.select = TRUE,
  title = "Required subsidies for all projects by 2030",
  secaxis.factor = secaxis.factor.woco2,
  ymax.annual = (subsidies.cumulative.2030.woco2 + 100) / secaxis.factor.woco2)

p.subsidies.wco2 <- temp[[1]]

# Align plots on the left
p.left <- align_plots(
  p.greenh2.natgas.woco2 + theme(legend.position = "none"),
  p.greenh2.greyh2.woco2 + theme(legend.position = "none"),
  p.subsidies.woco2 + theme(legend.position = "none"),
  align = "v")

# Align plots in the middle
p.middle <- align_plots(
  p.greenh2.natgas.wco2 + theme(legend.position = "none"),
  p.greenh2.greyh2.wco2 + theme(legend.position = "none"),
  p.subsidies.wco2 + theme(legend.position = "none"),
  align = "v")

# Plot enduses and capacity with curly brace
p.brace <- ggplot() +
  stat_brace(aes(x = c(0, 1),
                 y = c(0, 1)),
             rotate = 90,
             linewidth = 0.25,
             outside = FALSE) +
  theme_void()

# Helper function for end-use labels
drawEndUseLabel <- function(green, fossil) {
  enduse.select <- data.enduse.transportcost %>%
    filter(greenproduct == green, fossilcompetitor == fossil) %>%
    pull(enduse)

  cap <- data.iea.enduse %>% 
    filter(enduse %in% enduse.select, year == 2030) %>% 
    summarise(cumcap.sum = sum(cumcap.sum))

  label <- paste0("End-uses (", round(cap, 0), " GW):\n",
                  enduse.labels[enduse.select] %>% paste(collapse = "\n"))
  
  p <- ggdraw() + draw_label(label, size = 5, x = 0.05, hjust = 0)
}

p.row1.enduse <- drawEndUseLabel(green = "Green hydrogen", fossil = "Natural gas")

# First row: Green hydrogen vs. natural gas
p.row1 <- plot_grid(
  p.left[[1]],
  p.middle[[1]],
  NULL,
  p.brace,
  p.row1.enduse,
  ncol = 5,
  rel_widths = c(1, 1, -0.1, 0.1, 0.3),
  labels = c("a", "b", "", ""),
  label_size = font.size
)

# Plot enduses and capacity with curly brace
p.row2.enduse <- drawEndUseLabel(green = "Green hydrogen", fossil = "Grey hydrogen")

# Second row: Green hydrogen vs. grey hydrogen
p.row2 <- plot_grid(
  p.left[[2]],
  p.middle[[2]],
  NULL,
  p.brace,
  p.row2.enduse,
  ncol = 5,
  rel_widths = c(1, 1, -0.1, 0.1, 0.3),
  labels = c("c", "d", "", ""),
  label_size = font.size
)

# Cost gap panels
p.costpanels <- plot_grid(
  p.row1,
  p.row2,
  ncol = 1
)

# Cost gap panels with legend
p.costpanels.legend <- plot_grid(
  p.costpanels,
  get_legend(p.greenh2.natgas.woco2),
  rel_widths = c(2.3, 0.45)
)

# Third row: Required subsidies
p.row3 <- plot_grid(
  p.left[[3]],
  p.middle[[3]],
  p.subsidygap,
  ncol = 3,
  rel_widths = c(1, 1, 0.75),
  align = "h",
  labels = c("e", "f", "g"),
  label_size = font.size
)

# Fourth row: Legend for subsidies
p.row4 <- plot_grid(
  get_legend(p.subsidies.woco2 + theme(legend.box = "horizontal",
                                       legend.direction = "vertical",
                                       legend.justification = "centre")),
  ncol = 2,
  rel_widths = c(2, 0.75)
  )

# Add labels at the very top
p.woco2 <- ggdraw() +
  draw_label(aes(label = "Without carbon price"),
             size = 6,
             fontface = "bold")

p.wco2 <- ggdraw() +
  draw_label(aes(label = "With ambitious carbon price"),
             size = 6,
             fontface = "bold")
  
p.colnames <- plot_grid(
  p.woco2,
  p.wco2,
  NULL,
  ncol = 3,
  rel_widths = c(1,1, 0.75)
)

# Total plot
p <- plot_grid(
  p.colnames,
  p.costpanels.legend,
  p.row3,
  p.row4,
  rel_heights = c(0.1, 2, 1.2, 0.22),
  ncol = 1
)

ggsave("figures/FIGURE05_future_implementation_gap.png",
       width = 180, height = 160, units = "mm", bg = "white", dpi = 600)
ggsave("figures/FIGURE05_future_implementation_gap.svg",
       width = 180, height = 160, units = "mm", bg = "white")

```

## TABLE 1: Required cumulative subsidies

```{r table 1, message = FALSE}

# Subsidies required for projects announced until 2030
data.subsidies.annual %>% 
  filter(year.cap <= 2030) %>% 
  group_by(scenario, carbonprice) %>% 
  summarise(subsidies.required = sum(subsidies.annual)) %>%
  mutate(subsidies.announced = data.subsidies.tracker %>% 
           filter(region == "World",
                  publication == 2023) %>% 
           pull(value)) %>%
  mutate(subsidies.gap = subsidies.required - subsidies.announced) %>% 
  # Round numbers
  mutate(subsidies.required.trillion = round(1E-3*subsidies.required, 1),
         subsidies.announced.trillion = round(1E-3*subsidies.announced, 1),
         subsidies.gap.trillion = round(1E-3*subsidies.gap, 1))

# Subsidies required for 1.5°C scenario until 2050
data.subsidies.annual %>% 
  filter(year.cap <= 2050) %>% 
  group_by(scenario, carbonprice) %>% 
  summarise(subsidies.required = sum(subsidies.annual)) %>% 
  # Round numbers
  mutate(subsidies.required.trillion = round(1E-3*subsidies.required, 1))

# Specific subsidies required for 1.5°C scenario in 2050
data.subsidies.annual %>% 
  filter(year == 2050) %>% 
  mutate(prod = cap * flh.electrolysis * eta) %>% 
  group_by(scenario, carbonprice, year) %>%
  summarise(subsidies.relative = 1E9*sum(subsidies.annual) / sum(prod)) %>% 
  # Round numbers
  mutate(subsidies.relative.rounded = round(subsidies.relative, 0))

```

# EXTENDED DATA FIGURES

## EXTENDED DATA FIGURE 1: Electrolysis in 1.5°C scenarios

```{r extended data figure 1}

# Formatting
size.point <- 0.75
size.point.hist <- 1
lwd <- 0.5

# Institutional scenarios labels
data.label <- data.scenarios %>% 
  filter(year == 2050)

# Calculate required exponential growth rates for 2030 1.5°C scenarios
cap.2023 <- data.real %>% 
  filter(year == 2023) %>% 
  pull(value)

data.growth.inst <- data.scenarios %>% 
  filter(year == 2030) %>% 
  mutate(growth.rate = (value.mw/cap.2023)**(1/7) - 1) %>% 
  select(organisation, report, scen.name, region, value.mw, year, growth.rate) %>% 
  group_by(organisation, report, scen.name, region, growth.rate) %>% 
  complete(year = seq(2023, 2030, 0.25)) %>% 
  mutate(value.mw = cap.2023 * (1 + growth.rate)**(year - 2023))

data.growth.inst.minmax <- data.growth.inst %>% 
  group_by(year) %>% 
  summarise(min = min(value.mw),
            max = max(value.mw))

# Historical growth rate
growth.hist <- data.real %>% 
  summarise(growth.hist = (last(value)/first(value))**(1/(last(year - first(year)))) - 1) %>% 
  pull(growth.hist)

# Left panel until 2030
data.plot.scenarios <- data.scenarios %>% 
  filter(year %in% c(2025, 2030))

# Institutional
p.inst.left <- ggplot() +
  # Historical data
  geom_line(data = data.real,
            mapping = aes(x = year, y = value, color = "Realised projects"),
            linewidth = lwd) +
  geom_point(data = data.real,
             mapping = aes(x = year, y = value, color = "Realised projects"),
             size = size.point.hist) +
  scale_color_manual(name = NULL,
                     values = c("Realised projects" = "black")) +
  new_scale_color() +
  # Required growth for 2030
  geom_ribbon(data = data.growth.inst.minmax,
              mapping = aes(x = year, ymin = min, ymax = max),
              fill = "grey95") +
  geom_line(data = data.growth.inst,
            mapping = aes(x = year, y = value.mw, color = organisation, group = scen.name),
            linetype = "dotted",
            linewidth = 0.25) +
  # Institutional scenarios
  geom_line(data = data.plot.scenarios,
            mapping = aes(x = year, y = value.mw, color = organisation, group = scen.name),
            linewidth = lwd) +
  geom_point(data = data.plot.scenarios,
            mapping = aes(x = year, y = value.mw, color = organisation, group = scen.name),
            size = size.point) +
  # Annotation (2020-2023)
  geom_label_repel(
    mapping = aes(x = 2021.5, y = 0,
                  label = paste0("Historical growth rate\nfrom 2020\u20132023:\n",
                                 growth.hist %>% round(2) * 100,
                                 " %/yr")),
    size = 6 / .pt,
    nudge_y = 500,
    segment.size = 0.25,
    segment.curvature = -0.1,
    segment.shape = -1,
    segment.square = FALSE,
    min.segment.length = 0
  ) +
  # Annotation (2023-2030)
  geom_label_repel(
    mapping = aes(x = 2029, y = 400, 
                  label = paste0("Required growth rate \nfrom 2023\u20132030:\n",
                                 data.growth.inst %>% pull(growth.rate) %>% min %>% round(2) * 100,
                                 "\u2013",
                                 data.growth.inst %>% pull(growth.rate) %>% max %>% round(2) * 100,
                                 " %/yr")
                  ),
    size = 6 / .pt,
    nudge_x = -3,
    nudge_y = 1000,
    segment.size = 0.25,
    segment.curvature = 0.2,
    segment.shape = -1,
    segment.square = FALSE,
    segment.angle = 90,
    min.segment.length = 0) +
  scale_color_npg(name = "Scenario source") +
  scale_x_continuous(name = "Year", breaks = c(2020, 2025, 2030)) +
  scale_y_continuous(name = "Capacity (GW)", limits = c(0, 1730)) +
  ggtitle("Institutional and corporate 1.5°C scenarios until 2030") +
  theme(axis.line = element_blank(),
        panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5))

# Left panel until 2030
data.plot.scenarios <- data.scenarios %>% 
  filter(year %in% seq(2025, 2050))

p.inst.right <- ggplot() +
  geom_rect(mapping = aes(xmin = 2020, xmax = 2030, ymin = 0, ymax = 1730),
          fill = "transparent",
          color = "grey40",
          linewidth = 0.25) +
  # Historical data
  geom_line(data = data.real,
            mapping = aes(x = year, y = value, color = "Realised projects"),
            linewidth = lwd) +
  geom_point(data = data.real,
             mapping = aes(x = year, y = value, color = "Realised projects"),
             size = size.point.hist) +
  scale_color_manual(name = NULL,
                     values = c("Realised projects" = "black")) +
  new_scale_color() +
  # Required growth for 2030
  geom_line(data = data.growth.inst,
            mapping = aes(x = year, y = value.mw, color = organisation, group = scen.name),
            linetype = "dotted",
            linewidth = 0.25) +
  # Scenarios
  geom_line(data = data.plot.scenarios,
            mapping = aes(x = year, y = value.mw, color = organisation, group = scen.name),
            linewidth = lwd) +
  geom_point(data = data.plot.scenarios,
            mapping = aes(x = year, y = value.mw, color = organisation, group = scen.name),
            size = size.point) +
  scale_color_npg(name = "Scenario source") +
  scale_x_continuous(name = "Year",
                     breaks = seq(2020, 2050, 5),
                     limits = c(2020, 2055)) +
  scale_y_continuous(name = "Capacity (GW)",
                     limits = c(0, 7600)) +
  geom_text_repel(
    data = data.label,
    mapping = aes(x = 2050.3, y = value.mw, label = scen.name, color = organisation),
    size = 5 /.pt,
    nudge_x = 4,
    hjust = 0,
    direction = "y",
    segment.size = 0.25,
    min.segment.length = 0,
    xlim = c(2050, 2080)
  ) +
  ggtitle("Institutional and corporate 1.5°C scenarios until 2050") +
  coord_cartesian(clip = "off") +
  theme(plot.margin = unit(c(1,20,1,1), units = "mm"))

# IAMs: 2020-2030
data.plot.ar6 <- data.ar6 %>% 
  filter(period %in% seq(2020, 2030)) %>% 
  mutate(period = as.numeric(period))

data.plot.ar6.imp <- data.plot.ar6 %>% 
  filter(!is.na(imp_marker))

data.plot.ngfs <- data.ngfs %>% 
  filter(period %in% seq(2020, 2030)) %>% 
  mutate(modelscen = paste(model, scenario)) %>% 
  mutate(period = as.numeric(period))

ar6.colors <- pal_npg()(3)
names(ar6.colors) = c("LD", "Ren", "SP")

ar6.colors = c(ar6.colors, "C1" = "grey")

ar6.labels <- c("C1" = "C1 (1.5°C)",
                "LD" = "IMP-LD",
                "Ren" = "IMP-Ren",
                "SP" = "IMP-SP")

ngfs.colors <- pal_npg()(5)[4:5]
names(ngfs.colors) = c("Low demand", "Net Zero 2050")

# Calculate required growth
data.growth.iam <- data.plot.ar6 %>% 
  filter(period == 2030) %>% 
  mutate(growth.rate = (value/cap.2023)**(1/7) - 1) %>% 
  group_by(model, scenario, region, variable, growth.rate) %>% 
  complete(period = seq(2023, 2030, 0.25)) %>% 
  mutate(value = cap.2023 * (1 + growth.rate)**(period - 2023))

data.growth.iam.minmax <- data.growth.iam %>% 
  group_by(period) %>% 
  summarise(min = min(value),
            max = max(value))

# Historical growth rate
growth.hist <- data.real %>% 
  summarise(growth.hist = (last(value)/first(value))**(1/(last(year - first(year)))) - 1) %>% 
  pull(growth.hist)

# 2020-2050
p.iam.left <- ggplot() +
  # Required growth for 2030
  geom_ribbon(data = data.growth.iam.minmax,
              mapping = aes(x = period, ymin = min, ymax = max),
              fill = "grey95") +
  # AR6 spaghetti plot
  geom_line(data = data.plot.ar6,
            mapping = aes(x = period, y = value, color = "C1", group = modelscen),
            linewidth = 0.1) +
  # IMP scenarios
  geom_line(data = data.plot.ar6.imp,
            mapping = aes(x = period, y = value, color = imp_marker, group = modelscen),
            linewidth = lwd) +
  geom_point(data = data.plot.ar6.imp,
             mapping = aes(x = period, y = value, color = imp_marker, group = modelscen),
             size = size.point) +
  scale_color_manual(name = "IPCC AR6",
                     values = ar6.colors,
                     labels = ar6.labels,
                     guide = guide_legend(order = 1)) +
  # NGFS scenarios
  new_scale_color() +
  geom_line(data = data.plot.ngfs,
            mapping = aes(x = period, y = value, color = scenario, group = modelscen),
            linewidth = lwd) +
  geom_point(data = data.plot.ngfs,
             mapping = aes(x = period, y = value, color = scenario, group = modelscen),
             shape = 17,
             size = size.point) +
  scale_color_manual(name = "NGFS",
                     values = ngfs.colors,
                     guide = guide_legend(order = 2)) +
  new_scale_color() +
  # Historical data
  geom_line(data = data.real,
            mapping = aes(x = year, y = value, color = "Realised projects"),
            linewidth = lwd) +
  geom_point(data = data.real,
             mapping = aes(x = year, y = value, color = "Realised projects"),
             size = size.point.hist) +
  # Annotation (2020-2023)
  geom_label_repel(
    mapping = aes(x = 2021.5, y = 0,
                  label = paste0("Historical growth rate\nfrom 2020\u20132023:\n",
                                 growth.hist %>% round(2) * 100,
                                 " %/yr")),
    size = 6 / .pt,
    nudge_y = 500,
    segment.size = 0.25,
    segment.curvature = -0.1,
    segment.shape = -1,
    segment.square = FALSE,
    min.segment.length = 0
  ) +
  # Annotation (2023-2030)
  geom_label_repel(
    mapping = aes(x = 2029.5, y = 300, 
                  label = paste0("Required growth rate \nfrom 2023\u20132030:\n",
                                 data.growth.iam %>% pull(growth.rate) %>% min %>% round(2) * 100,
                                 "\u2013",
                                 data.growth.iam %>% pull(growth.rate) %>% max %>% round(2) * 100,
                                 " %/yr")
                  ),
    size = 6 / .pt,
    nudge_x = -3,
    nudge_y = 1000,
    segment.size = 0.25,
    segment.curvature = 0.2,
    segment.shape = -1,
    segment.square = FALSE,
    segment.angle = 90,
    min.segment.length = 0) +
  scale_color_manual(name = NULL,
                     values = c("Realised projects" = "black")) +
  scale_x_continuous(name = "Year", breaks = c(2020, 2025, 2030)) +
  scale_y_continuous(name = "Capacity (GW)", limits = c(0, 1730)) +
  ggtitle("IAM 1.5°C scenarios until 2030") +
  theme(axis.line = element_blank(),
        panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5))

# IAMs: 2020-2050
data.plot.ar6.right <- data.ar6 %>% 
  filter(period %in% seq(2020, 2050)) %>% 
  mutate(period = as.numeric(period))

data.plot.ar6.imp.right <- data.plot.ar6.right %>% 
  filter(!is.na(imp_marker))

data.plot.ngfs.right <- data.ngfs %>% 
  filter(period %in% seq(2020, 2050)) %>% 
  mutate(modelscen = paste(model, scenario)) %>% 
  mutate(period = as.numeric(period))

p.iam.right <- ggplot() +
  geom_rect(mapping = aes(xmin = 2020, xmax = 2030, ymin = 0, ymax = 1730),
          fill = "transparent",
          color = "grey40",
          linewidth = 0.25) +
  # AR6 spaghetti plot
  geom_line(data = data.plot.ar6.right,
            mapping = aes(x = period, y = value, color = "C1", linewidth = imp_marker, group = modelscen),
            linewidth = 0.1) +
  # IMP scenarios
  geom_line(data = data.plot.ar6.imp.right,
            mapping = aes(x = period, y = value, color = imp_marker, group = modelscen),
            linewidth = lwd) +
  geom_point(data = data.plot.ar6.imp.right,
             mapping = aes(x = period, y = value, color = imp_marker, group = modelscen),
             size = size.point) +
  scale_color_manual(name = "IPCC AR6",
                     values = ar6.colors,
                     labels = ar6.labels,
                     guide = guide_legend(order = 1)) +
  # NGFS scenarios
  new_scale_color() +
  geom_line(data = data.plot.ngfs.right,
            mapping = aes(x = period, y = value, color = scenario, group = modelscen),
            linewidth = lwd) +
  geom_point(data = data.plot.ngfs.right,
             mapping = aes(x = period, y = value, color = scenario, group = modelscen),
             shape = 17,
             size = size.point) +
  scale_color_manual(name = "NGFS",
                     values = ngfs.colors,
                     guide = guide_legend(order = 2)) +
  new_scale_color() +
  # Historical data
  geom_line(data = data.real,
            mapping = aes(x = year, y = value, color = "Realised projects"),
            linewidth = lwd) +
  geom_point(data = data.real,
             mapping = aes(x = year, y = value, color = "Realised projects"),
             size = size.point.hist) +
  scale_color_manual(name = NULL,
                     values = c("Realised projects" = "black")) +
  scale_x_continuous(name = "Year",
                     breaks = seq(2020, 2050, 5),
                     limits = c(2020, 2055)) +
  scale_y_continuous(name = "Capacity (GW)",
                     limits = c(0, 7600)) +
  ggtitle("IAM 1.5°C scenarios until 2050") +
  coord_cartesian(clip = "off") +
  theme(plot.margin = unit(c(1,20,1,1), units = "mm"))

p.panels <- plot_grid(
  p.inst.left + theme(legend.position = "none"),
  p.inst.right + theme(legend.position = "none"),
  p.iam.left + theme(legend.position = "none"),
  p.iam.right + theme(legend.position = "none"),
  ncol = 2,
  align = "h",
  labels = "auto",
  label_size = font.size,
  rel_widths = c(6, 10)
)

p.legend <- plot_grid(
  get_legend(p.inst.left),
  get_legend(p.iam.left),
  ncol = 1
)

# Total plot
p <- plot_grid(
  p.panels,
  p.legend,
  ncol = 2,
  rel_widths = c(1, 0.15)
)

ggsave("figures/EXT_DATA_FIGURE01_electrolysis_1.5C_scenarios.png",
       width = 180, height = 120, units = "mm", bg = "white", dpi = 600)
ggsave("figures/EXT_DATA_FIGURE01_electrolysis_1.5C_scenarios.svg",
       width = 180, height = 120, units = "mm", bg = "white")

```

## EXTENDED DATA FIGURE 2: Designated end-uses of green hydrogen projects

```{r extended data figure 2}

# Colour palette with 15 colours
enduse.colors <- colorRampPalette(brewer.pal(8, "Set2"))(15)

# Arrange by size of fossil competitor
fossil.order <- data.iea.enduse %>%
  filter(year == 2030) %>%
  group_by(fossilcompetitor) %>%
  summarise(cumcap.sum = sum(cumcap.sum)) %>%
  arrange(-cumcap.sum) %>%
  pull(fossilcompetitor)

# Arrange by size of green product
enduse.order <- data.iea.enduse %>% 
  filter(year == 2030) %>% 
  group_by(enduse) %>% 
  summarise(cumcap.sum = sum(cumcap.sum)) %>% 
  arrange(-cumcap.sum) %>% 
  pull(enduse)

# Arrange by both
enduse.order.both <- data.iea.enduse %>% 
  filter(year == 2030) %>% 
  arrange(match(fossilcompetitor, fossil.order),
          match(enduse, enduse.order)) %>% 
  pull(enduse) %>% 
  unique() %>%
  rev()

# Plot until 2024
data.plot <- data.iea.enduse %>% 
  filter(year %in% seq(2020, 2024)) %>% 
  mutate(enduse = factor(enduse, levels = enduse.order.both)) %>% 
  group_by(enduse, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

p.left <- ggplot() + 
  geom_col(data = data.plot,
           mapping = aes(x = year, y = cumcap.sum, fill = enduse)) +
  scale_fill_manual(name = "End use",
                    values = enduse.colors,
                    labels = enduse.labels) +
  scale_x_continuous(name = "Year",
                     breaks = seq(2020, 2024)) +
  ylab("Capacity (GW)") +
  ggtitle("Projects until 2024 by end use")

# Plot from 2024
data.plot <- data.iea.enduse %>% 
  filter(year %in% seq(2024, 2030)) %>% 
  mutate(enduse = factor(enduse, levels = enduse.order.both)) %>% 
  group_by(enduse, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

# Anchor points for braces
natgas.2030.min <- 0
natgas.2030.max <- data.iea.enduse %>% 
  filter(year == 2030,
         fossilcompetitor == "Natural gas") %>% 
  pull(cumcap.sum) %>% 
  sum()

greyh2.2030.min <- natgas.2030.max
greyh2.2030.max <- data.iea.enduse %>% 
  filter(year == 2030,
         fossilcompetitor %in% c("Natural gas", "Grey hydrogen")) %>% 
  pull(cumcap.sum) %>% 
  sum()

diesel.2030.min <- greyh2.2030.max
diesel.2030.max <- data.iea.enduse %>% 
  filter(year == 2030,
         fossilcompetitor %in% c("Natural gas", "Grey hydrogen", "Diesel")) %>% 
  pull(cumcap.sum) %>% 
  sum()

kerosene.2030.min <- diesel.2030.max
kerosene.2030.max <- data.iea.enduse %>% 
  filter(year == 2030,
         fossilcompetitor %in% c("Natural gas", "Grey hydrogen", "Diesel", "Kerosene")) %>% 
  pull(cumcap.sum) %>% 
  sum()

greymeoh.2030.min <- kerosene.2030.max
greymeoh.2030.max <- data.iea.enduse %>% 
  filter(year == 2030) %>% 
  pull(cumcap.sum) %>% 
  sum()

data.label <- tribble(
   ~ value.mw, ~ label,
   mean(c(diesel.2030.min, diesel.2030.max)), "Fossil competitor: Diesel",
   mean(c(kerosene.2030.min, kerosene.2030.max)), "Fossil competitor: Kerosene",
   mean(c(greymeoh.2030.min, greymeoh.2030.max)), "Fossil competitor: Grey methanol"
)

p.right <- ggplot() +
  geom_col(data = data.plot,
           mapping = aes(x = year, y = cumcap.sum, fill = enduse)) +
  scale_fill_manual(name = "End use",
                    values = enduse.colors,
                    labels = enduse.labels) +
  # Add braces for fossil competitors
  stat_brace(aes(x = c(2030.5, 2031),
                 y = c(0, natgas.2030.max)),
             rotate = 90,
             linewidth = 0.25,
             outside = FALSE) +
  stat_bracetext(aes(x = c(2030.5, 2031),
                     y = c(0, natgas.2030.max)),
                 rotate = 90,
                 outside = FALSE,
                 label = "Fossil competitor:\nNatural gas",
                 size = 5 / .pt) +
  stat_brace(aes(x = c(2030.5, 2031),
                 y = c(greyh2.2030.min, greyh2.2030.max)),
             rotate = 90,
             linewidth = 0.25,
             outside = FALSE) +
  stat_bracetext(aes(x = c(2030.5, 2031),
                     y = c(greyh2.2030.min, greyh2.2030.max)),
                 rotate = 90,
                 outside = FALSE,
                 label = "Fossil competitor:\nGrey hydrogen",
                 size = 5 / .pt) +
  # Add text labels for other end uses
  geom_text_repel(
    data = data.label,
    mapping = aes(x = 2030.5, y = value.mw, label = label),
    hjust = 0,
    size = 5 /.pt,
    nudge_x = 0.55,
    direction = "y",
    box.padding = 0.2,
    segment.size = 0.25,
    min.segment.length = 0
  ) +
  scale_x_continuous(name = "Year",
                     breaks = seq(2024, 2030),
                     limits = c(NA, 2034.6)) +
  ylab("Capacity (GW)") +
  ggtitle("Projects from 2024 by end use")    

p <- plot_grid(
  p.left + theme(legend.position = "none"),
  p.right + theme(legend.position = "none"),
  get_legend(p.right),
  ncol = 3,
  labels = c("a", "b"),
  label_size = font.size,
  rel_widths = c(4, 8.5, 2.2)
  )

print(p)

ggsave("figures/EXT_DATA_FIGURE02_enduse.png",
       width = 180, height = 70, units = "mm", bg = "white", dpi = 600)
ggsave("figures/EXT_DATA_FIGURE02_enduse.svg",
       width = 180, height = 70, units = "mm", bg = "white")

```

## EXTENDED DATA FIGURE 3: Illustration of modelling pay-as-bid market premium

```{r extended data figure 3}

# Use green hydrogen vs. natural gas competition
data.plot.green <- data.lcox %>% 
  filter(scenario == "Default",
         variable == "LCOX|Green hydrogen",
         enduse == "otherind")

data.plot.fossil <- data.fossilcosts %>% 
  filter(scenario == "Default",
         variable == "Cost|Natural gas|w/ CO2",
         enduse == "otherind")

# For area plot
data.plot.gap <- bind_rows(
  data.plot.green %>% 
    select(scenario, variable, region, unit, period, value),
  data.plot.fossil %>% 
    select(scenario, variable, region, unit, period, value)
) %>% 
  pivot_wider(names_from = variable,
              values_from = value) %>% 
  rename(lcoh = `LCOX|Green hydrogen`,
         fossilcost = `Cost|Natural gas|w/ CO2`)

# Cost gap for projects built in 2024
data.plot.gap.2024 <- data.plot.gap %>% 
  filter(period %in% seq(2024, 2039)) %>% 
  mutate(lcoh = first(lcoh))

# Cost gap for projects built in 2035
data.plot.gap.2035 <- data.plot.gap %>% 
  filter(period %in% seq(2035, 2050)) %>% 
  mutate(lcoh = first(lcoh))

# Helper function
# Find intersection between two lines (L1, L2), where L1: (x1,y1)-(x2,y2) and
# L2: (x3,y3)-(x4,y4) and x1 = x3 and x2 = x4 (in years)
findIntersection <- function(period, L1, L2) {
  x1 <- period[1]
  x2 <- period[2]
  y1 <- L1[1]
  y2 <- L1[2]
  y3 <- L2[1]
  y4 <- L2[2]
  x <- (x1*y2 - y1*x2 - x1*y4 + y3*x2)/(y3 - y4 - y1 + y2)
  y <- ((x1*y2 - y1*x2)*(y3-y4) - (y1-y2)*(x1*y4-y3*x2))/((x1-x2)*(y3-y4-y1+y2))
  return(c(x,y))
}

# Find intersection of green with fossil
temp <- bind_rows(
  data.plot.gap %>% filter(lcoh>fossilcost) %>% slice(n()),
  data.plot.gap %>% filter(lcoh<fossilcost) %>% slice(1))

intersect.greenfossil <- findIntersection(temp$period,
                                          temp$lcoh,
                                          temp$fossilcost)

# Find intersection of green LCOH in 2035 with fossil
temp <- bind_rows(
  data.plot.gap.2035 %>% filter(lcoh>fossilcost) %>% slice(n()),
  data.plot.gap.2035 %>% filter(lcoh<fossilcost) %>% slice(1)
)

intersect.green2035fossil <- findIntersection(temp$period,
                                              temp$lcoh,
                                              temp$fossilcost)

p <- ggplot() +
  # Vertical lines
  geom_line(aes(x = c(2024, 2024),
                y = c(-5, data.plot.green %>% filter(period == 2024) %>% pull(value) + 40)),
            color = "grey80",
            linetype = "dashed",
            linewidth = 0.25) +
  geom_line(aes(x = c(2039, 2039),
                y = c(-5, data.plot.green %>% filter(period == 2024) %>% pull(value) + 40)),
            color = "grey80",
            linetype = "dashed",
            linewidth = 0.25) +
  geom_line(aes(x = c(2035, 2035),
                y = c(-5, data.plot.green %>% filter(period == 2024) %>% pull(value) + 20)),
            color = "grey80",
            linetype = "dashed",
            linewidth = 0.25) +
  geom_line(aes(x = c(2050, 2050),
                y = c(-5, data.plot.green %>% filter(period == 2024) %>% pull(value) + 20)),
            color = "grey80",
            linetype = "dashed",
            linewidth = 0.25) +
  geom_line(aes(x = c(intersect.greenfossil[1], intersect.greenfossil[1]),
                y = c(-5, intersect.greenfossil[2])),
            color = "grey80",
            linetype = "dashed",
            linewidth = 0.25) +
  geom_line(aes(x = c(intersect.green2035fossil[1], intersect.green2035fossil[1]),
                y = c(-5, intersect.green2035fossil[2])),
            color = "grey80",
            linetype = "dashed",
            linewidth = 0.25) +
  # LCOH and natural gas costs
  geom_line(data = data.plot.green,
            mapping = aes(x = period, y = value, color = "Green hydrogen")) +
  geom_line(data = data.plot.fossil,
            mapping = aes(x = period, y = value, color = "Natural gas")) +
  scale_color_manual(
    name = NULL,
    values = c("Green hydrogen" = "#76AF63",
               "Natural gas" = "#6A5541"),
    labels = c("Green hydrogen" = "Green hydrogen costs (LCOH)",
               "Natural gas" = "Natural gas prices (with carbon price)")
  ) + 
  # Add constant LCOH from 2024 to 2038
  geom_line(mapping = aes(x = c(2024, 2039),
                          y = c(data.plot.green %>% filter(period == 2024) %>% pull(value),
                                data.plot.green %>% filter(period == 2024) %>% pull(value))),
            color = "#76AF63",
            linetype = "dashed") +
  geom_text(mapping = aes(x = 2039.5, y = data.plot.green %>% filter(period == 2024) %>% pull(value)),
            label = "LCOH of projects\nbuilt in 2024",
            lineheight = 1,
            size = 6 / .pt,
            hjust = 0) +
  # Add horizontal errorbar to show payback period
  geom_errorbarh(mapping = aes(xmin = 2024, xmax = 2039, y = data.plot.green %>% filter(period == 2024) %>% pull(value) + 40),
                 height = 10,
                 linewidth = 0.25) +
  geom_text(mapping = aes(x = 2031.5, y = data.plot.green %>% filter(period == 2024) %>% pull(value) + 45),
            label = "Payback period: 15 years",
            size = 6 / .pt,
            hjust = 0.5) +
  geom_errorbarh(mapping = aes(xmin = 2035, xmax = 2050, y = data.plot.green %>% filter(period == 2024) %>% pull(value) + 20),
                 height = 10,
                 linewidth = 0.25) +
  geom_text(mapping = aes(x = 2042.5, y = data.plot.green %>% filter(period == 2024) %>% pull(value) + 25),
            label = "Payback period: 15 years",
            size = 6 / .pt,
            hjust = 0.5) +
  # Add constant LCOH from 2035 to 2049
  geom_line(mapping = aes(x = c(2035, 2050),
                          y = c(data.plot.green %>% filter(period == 2035) %>% pull(value),
                                data.plot.green %>% filter(period == 2035) %>% pull(value))),
            color = "#76AF63",
            linetype = "dashed") +
  geom_text(mapping = aes(x = 2050.5, y = data.plot.green %>% filter(period == 2035) %>% pull(value)),
            label = "LCOH of projects\nbuilt in 2035",
            lineheight = 1,
            size = 6 / .pt,
            hjust = 0) +
  # Add ribbon for cost gap
  geom_ribbon(data = data.plot.gap.2024,
              mapping = aes(x = period, ymin = lcoh, ymax = fossilcost),
              fill = "#C10505",
              alpha = 0.1) +
  geom_ribbon(data = data.plot.gap.2035,
              mapping = aes(x = period, ymin = lcoh, ymax = pmin(fossilcost, lcoh)),
              fill = "#892F71",
              alpha = 0.1) +
  # Add arrows for cost gap (2024 projects)
  geom_segment(
    aes(x = 2024+2.5,
        xend = 2024+2.5,
        y = data.plot.green %>% filter(period == 2024) %>% pull(value) - 3,
        yend = data.plot.fossil %>% filter(period == 2027) %>% pull(value) + 3),
    arrow = arrow(length = unit(1, "mm"), ends = "both"),
    color = "#C10505",
    linewidth = 0.25
    ) +
  geom_segment(
    aes(x = 2024+7.5,
        xend = 2024+7.5,
        y = data.plot.green %>% filter(period == 2024) %>% pull(value) - 3,
        yend = data.plot.fossil %>% filter(period == 2032) %>% pull(value) + 3),
    arrow = arrow(length = unit(1, "mm"), ends = "both"),
    color = "#C10505",
    linewidth = 0.25
    ) + 
  geom_segment(
    aes(x = 2024+12.5,
        xend = 2024+12.5,
        y = data.plot.green %>% filter(period == 2024) %>% pull(value) - 3,
        yend = data.plot.fossil %>% filter(period == 2037) %>% pull(value) + 3),
    arrow = arrow(length = unit(1, "mm"), ends = "both"),
    color = "#C10505",
    linewidth = 0.25
    ) +
  # Add arrows for cost gap (2035 projects)
  geom_segment(
    aes(x = 2035+2.5,
        xend = 2035+2.5,
        y = data.plot.green %>% filter(period == 2035) %>% pull(value) - 3,
        yend = data.plot.fossil %>% filter(period == 2038) %>% pull(value) + 3),
    arrow = arrow(length = unit(1, "mm"), ends = "both"),
    color = "#892F71",
    linewidth = 0.25
    ) +
  geom_segment(
    aes(x = 2035+7.5,
        xend = 2035+7.5,
        y = data.plot.green %>% filter(period == 2035) %>% pull(value) - 3,
        yend = data.plot.fossil %>% filter(period == 2043) %>% pull(value) + 3),
    arrow = arrow(length = unit(1, "mm"), ends = "both"),
    color = "#892F71",
    linewidth = 0.25
    ) +
  # Annotations
  geom_label(
    aes(x = 2045, y = 50),
    label = "Projects built after 2043\ndo not require any subsidies",
    size = 6 / .pt,
    lineheight = 1
  ) +
  geom_label(
    aes(x = 2030, y = 30),
    label = "Projects built in 2024\nrequire subsidies throughout\ntheir payback period",
    size = 6 / .pt,
    lineheight = 1,
    color = "#C10505"
  ) +
  geom_label(
    aes(x = 2045, y = 150),
    label = "Projects built in 2035\nrequire subsidies until 2047",
    size = 6 / .pt,
    lineheight = 1,
    color =  "#892F71",
  ) +
  # Formatting
  scale_x_continuous(name = "Year",
                     breaks = c(2024, seq(2030, 2050, 5)),
                     limits = c(2024, 2054)) +
  scale_y_continuous(name = expression("$/MWh"[LHV]),
                     limits = c(-5, 240),
                     expand = c(0, 0)) +
  ggtitle("Modelling the pay-as-bid market premium")

ggsave("figures/EXT_DATA_FIGURE03_illustration_model.png",
       width = 160, height = 80, units = "mm", bg = "white", dpi = 600)
ggsave("figures/EXT_DATA_FIGURE03_illustration_model.svg",
       width = 160, height = 80, units = "mm", bg = "white")

```

## EXTENDED DATA FIGURE 4: Green LCOH in 2030 compared to recent studies

```{r extended data figure 4}

data.lcoh <- read_excel("data/lcoh_studies.xlsx") %>% 
  pivot_longer(`2023`:`2050`, names_to = "year", values_to = "value") %>% 
  # Transform to same real 2023 USD
  mutate(
    value = ifelse(`Currency base year` == 2020, 1.18*value, value),
    value = ifelse(`Currency base year` == 2022, 1.04*value, value),
    value = ifelse(Unit == "EUR/kg", 1.09*value, value),
    unit = "2023USD/kg")

data.plot.lcoh <- data.lcoh %>% 
  filter(year == 2030) %>% 
  select(Organisation, `Short name`, Variable, value) %>%
  # Transform into $/MWh
  mutate(value = 1E3 * value / 33.33) %>% 
  pivot_wider(names_from = Variable, values_from = value)
  
data.plot.lcoh.own <- data.lcox %>% 
  filter(variable == "LCOX|Green hydrogen",
         period == 2030) %>% 
  # For default, use value with transport costs
  group_by(scenario) %>% 
  reframe(value = case_when(scenario == "Default" ~ value[is.na(enduse)],
                           scenario == "Progressive" ~ min(value),
                           scenario == "Conservative" ~ max(value))) %>% 
  group_by(scenario) %>% 
  summarise(value = unique(value)) %>% 
  pivot_wider(names_from = scenario, values_from = value)

p.lcoh <- ggplot() +
  # Grey background for range
  geom_rect(data = data.plot.lcoh.own,
            mapping = aes(xmin = 0.4, xmax = 14.6, ymin = Progressive, ymax = Conservative),
            fill = "grey",
            alpha = 0.3) +
  # Our own analysis
  geom_point(data = data.plot.lcoh.own,
             mapping = aes(x = "A", y = Default)) +
  geom_linerange(data = data.plot.lcoh.own,
                 mapping = aes(x = "A", ymin = Progressive, ymax = Conservative)) +
  geom_point(data = data.plot.lcoh,
             mapping = aes(x = `Short name`, y = `LCOH|Green|central`, color = Organisation)) +
  geom_linerange(data = data.plot.lcoh,
                 mapping = aes(x = `Short name`, ymin = `LCOH|Green|min`, ymax = `LCOH|Green|max`, color = Organisation)) +
  scale_color_npg() +
  geom_vline(xintercept = 1.5, linetype = "dashed", linewidth = 0.25) +
  scale_x_discrete(name = NULL,
                   labels = c("A" = "This study",
                              "BCG Europe" = "BCG\nEurope",
                              "Hydrogen Council Global" = "Hydrogen Council\nGlobal",
                              "BNEF Germany" = "BNEF\nGermany",
                              "BNEF US" = "BNEF\nUS",
                              "Capgemini Survey Global" = "Capgemini\nSurvey\nGlobal",
                              "DNV Dedicated Global" = "DNV\nDedicated\nGlobal",
                              "DNV Grid Global" = "DNV\nGrid\nGlobal",
                              "ICCT Central EU" = "ICCT\nCentral\nEU",
                              "ICCT Central US" = "ICCT\nCentral\nUS",
                              "IEA NZE Onshore Wind Global" = "IEA NZE\nOnshore Wind\nGlobal",
                              "IEA NZE Solar Global" = "IEA NZE\nSolar\nGlobal",
                              "NPC NZ2050 Delivered Costs US Gulf Coast" = "NPC NZ2050\nDelivered Costs\nUS Gulf Coast",
                              "CE Delft & TNO Dutch Wind Offshore" = "CE Delft & TNO\nOffshore Wind\nNetherlands")) +
  scale_y_continuous(name = expression("LCOH in 2030 ($/MWh"[LHV]*")"),
                     sec.axis = sec_axis(~.*33.33/1E3,
                                         name = expression("LCOH in 2030 ($/kgH"[2]*")")),
                     limits = c(0, NA)) +
  ggtitle("Levelised cost of green hydrogen in 2030") +
  theme(legend.position = "none")

ggsave("figures/EXT_DATA_FIGURE04_lcoh_comparison.png",
       width = 180, height = 80, units = "mm", bg = "white", dpi = 600)
ggsave("figures/EXT_DATA_FIGURE04_lcoh_comparison.svg",
       width = 180, height = 80, units = "mm", bg = "white")

```

## EXTENDED DATA FIGURE 5: Sectoral competition of green products with fossil competitors in all scenarios

```{r extended data figure 5}

# Green hydrogen vs. natural gas
p.greenh2.natgas.range.woco2 <- plotCostGapRange(
  carbonprice.select = FALSE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Natural gas",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE)

p.greenh2.natgas.range.wco2 <- plotCostGapRange(
  carbonprice.select = TRUE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Natural gas",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE)

# Green hydrogen vs. grey hydrogen
p.greenh2.greyh2.range.woco2 <- plotCostGapRange(
  carbonprice.select = FALSE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Grey hydrogen",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE)

p.greenh2.greyh2.range.wco2 <- plotCostGapRange(
  carbonprice.select = TRUE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Grey hydrogen",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE)

# Green hydrogen vs. Diesel
p.greenh2.diesel.range.woco2 <- plotCostGapRange(
  carbonprice.select = FALSE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Diesel",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Fossil competitor",,
  title.noco2 = TRUE,
  adjust.lcox = 1.1)  # Account for improved end-use efficiency in mobility

p.greenh2.diesel.range.wco2 <- plotCostGapRange(
  carbonprice.select = TRUE,
  greenproduct.select = "Green hydrogen",
  fossilcompetitor.select = "Diesel",
  greenproduct.label = "Green hydrogen",
  fossilcompetitor.label = "Fossil competitor",,
  title.noco2 = TRUE,
  adjust.lcox = 1.1)  # Account for improved end-use efficiency in mobility

# E-Kerosene vs. Kerosene
p.ekerosene.kerosene.range.woco2 <- plotCostGapRange(
  carbonprice.select = FALSE,
  greenproduct.select = "E-Kerosene",
  fossilcompetitor.select = "Kerosene",
  greenproduct.label = "Electrofuel",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE)

p.ekerosene.kerosene.range.wco2 <- plotCostGapRange(
  carbonprice.select = TRUE,
  greenproduct.select = "E-Kerosene",
  fossilcompetitor.select = "Kerosene",
  greenproduct.label = "Electrofuel",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE)

# E-Methanol vs. Methanol
p.emethanol.greymethanol.range.woco2 <- plotCostGapRange(
  carbonprice.select = FALSE,
  greenproduct.select = "E-Methanol",
  fossilcompetitor.select = "Grey methanol",
  greenproduct.label = "Electrofuel",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE)

p.emethanol.greymethanol.range.wco2 <- plotCostGapRange(
  carbonprice.select = TRUE,
  greenproduct.select = "E-Methanol",
  fossilcompetitor.select = "Grey methanol",
  greenproduct.label = "Electrofuel",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE)

# E-Methane vs. natural gas
p.emethane.natgas.range.woco2 <- plotCostGapRange(
  carbonprice.select = FALSE,
  greenproduct.select = "E-Methane",
  fossilcompetitor.select = "Natural gas",
  greenproduct.label = "Electrofuel",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE,
  remove.xaxis.labels = FALSE)

p.emethane.natgas.range.wco2 <- plotCostGapRange(
  carbonprice.select = TRUE,
  greenproduct.select = "E-Methane",
  fossilcompetitor.select = "Natural gas",
  greenproduct.label = "Electrofuel",
  fossilcompetitor.label = "Fossil competitor",
  title.noco2 = TRUE,
  remove.xaxis.labels = FALSE)

# Arrange plots

# First row
p.brace <- ggplot() +
  stat_brace(aes(x = c(0, 1),
                 y = c(0, 1)),
             rotate = 90,
             linewidth = 0.25,
             outside = FALSE) +
  theme_void()

p.row1.enduse <- drawEndUseLabel(green = "Green hydrogen", fossil = "Natural gas")

p.row1 <- plot_grid(
  p.greenh2.natgas.range.woco2 + theme(legend.position = "none"),
  p.greenh2.natgas.range.wco2 + theme(legend.position = "none"),
  NULL,
  p.brace,
  p.row1.enduse,
  ncol = 5,
  rel_widths = c(1, 1, 0, 0.1, 0.3),
  labels = c("a", "b"),
  label_size = font.size
)

# Second row
p.row2.enduse <- drawEndUseLabel(green = "Green hydrogen", fossil = "Grey hydrogen")

p.row2 <- plot_grid(
  p.greenh2.greyh2.range.woco2 + theme(legend.position = "none"),
  p.greenh2.greyh2.range.wco2 + theme(legend.position = "none"),
  NULL,
  p.brace,
  p.row2.enduse,
  ncol = 5,
  rel_widths = c(1, 1, 0, 0.1, 0.3),
  labels = c("c", "d"),
  label_size = font.size
)

# Third row
p.row3.enduse <- drawEndUseLabel(green = "Green hydrogen", fossil = "Diesel")

p.row3 <- plot_grid(
  p.greenh2.diesel.range.woco2 + theme(legend.position = "none"),
  p.greenh2.diesel.range.wco2 + theme(legend.position = "none"),
  NULL,
  p.brace,
  p.row3.enduse,
  ncol = 5,
  rel_widths = c(1, 1, 0, 0.1, 0.3),
  labels = c("e", "f"),
  label_size = font.size
)

# Fourth row
p.row4.enduse <- drawEndUseLabel(green = "E-Kerosene", fossil = "Kerosene")

p.row4 <- plot_grid(
  p.ekerosene.kerosene.range.woco2 + theme(legend.position = "none"),
  p.ekerosene.kerosene.range.wco2 + theme(legend.position = "none"),
  NULL,
  p.brace,
  p.row4.enduse,
  ncol = 5,
  rel_widths = c(1, 1, 0, 0.1, 0.3),
  labels = c("g", "h"),
  label_size = font.size
)

# Fifth row
p.row5.enduse <- drawEndUseLabel(green = "E-Methanol", fossil = "Grey methanol")

p.row5 <- plot_grid(
  p.emethanol.greymethanol.range.woco2 + theme(legend.position = "none"),
  p.emethanol.greymethanol.range.wco2 + theme(legend.position = "none"),
  NULL,
  p.brace,
  p.row5.enduse,
  ncol = 5,
  rel_widths = c(1, 1, 0, 0.1, 0.3),
  labels = c("i", "j"),
  label_size = font.size
)

# Sixth row
p.row6.enduse <- drawEndUseLabel(green = "E-Methane", fossil = "Natural gas")

p.row6 <- plot_grid(
  p.emethane.natgas.range.woco2 + theme(legend.position = "none"),
  p.emethane.natgas.range.wco2 + theme(legend.position = "none"),
  NULL,
  p.brace,
  p.row6.enduse,
  ncol = 5,
  rel_widths = c(1, 1, 0, 0.1, 0.3),
  labels = c("k", "l"),
  label_size = font.size
)

# Add labels at the very top
p.woco2 <- ggdraw() +
  draw_label(aes(label = "Without carbon price"),
             size = 6,
             fontface = "bold")

p.wco2 <- ggdraw() +
  draw_label(aes(label = "With ambitious carbon price"),
             size = 6,
             fontface = "bold")
  
p.colnames <- plot_grid(
  p.woco2,
  p.wco2,
  NULL,
  ncol = 3,
  rel_widths = c(1, 1, 0.4)
)

# Arrange panels
p.panels <- plot_grid(
  p.colnames,
  p.row1,
  p.row2,
  p.row3,
  p.row4,
  p.row5,
  p.row6,
  ncol = 1,
  rel_heights = c(0.15, 1, 1, 1, 1, 1, 1.07)
)

# x label
p.xlabel <- ggdraw() +
  draw_label("Year", size = font.size, x = 0.43)

p.panels.axis <- plot_grid(
  p.panels,
  p.xlabel,
  ncol = 1,
  rel_heights = c(1, 0.02)
)

# y label
p.ylabel <- ggdraw() +
  draw_label(expression("Specific cost ($/MWh"[LHV]*")"), size = font.size,
             angle = 90)

# Legend column
leg.h2 <- get_legend(p.greenh2.natgas.range.woco2)
leg.efuels <- get_legend(p.ekerosene.kerosene.range.woco2)

p.leg <- plot_grid(
  leg.h2,
  leg.efuels,
  ncol = 1
)

p <- plot_grid(
  p.ylabel,
  p.panels.axis,
  NULL,
  p.leg,
  ncol = 4,
  rel_widths = c(0.03, 1, 0.03, 0.19)
)

ggsave("figures/EXT_DATA_FIGURE05_cost_gap_by_enduse.png",
       width = 180, height = 180, units = "mm", bg = "white", dpi = 600)
ggsave("figures/EXT_DATA_FIGURE05_cost_gap_by_enduse.svg",
       width = 180, height = 180, units = "mm", bg = "white")

```

## EXTENDED DATA FIGURE 6: Required subsidies in all scenarios

```{r extended data figure 6}

# Conservative scenario
p.conservative.woco2 <- plotRequiredSubsidies(
  scenario.select = "Conservative",
  carbonprice.select = FALSE,
  title = "Required subsidies for all projects by 2030\n(conservative)")

p.conservative.wco2 <- plotRequiredSubsidies(
  scenario.select = "Conservative",
  carbonprice.select = TRUE,
  title = "Required subsidies for all projects by 2030\n(conservative)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

p.conservative.subsidygap <- plotSubsidyGap(
  scenario.select = "Conservative",
  subsidies.cumulative.2030.woco2 = p.conservative.woco2[[4]],
  subsidies.cumulative.2030.wco2 = p.conservative.wco2[[4]],
  ymax = p.conservative.woco2[[4]],
  title = "Required cumulative subsidies\nfor all projects by 2030 (conservative)")

# Default scenario
p.default.woco2 <- plotRequiredSubsidies(
  scenario.select = "Default",
  carbonprice.select = FALSE,
  title = "Required subsidies for all projects by 2030\n(central)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

p.default.wco2 <- plotRequiredSubsidies(
  scenario.select = "Default",
  carbonprice.select = TRUE,
  title = "Required subsidies for all projects by 2030\n(central)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

p.default.subsidygap <- plotSubsidyGap(
  scenario.select = "Default",
  subsidies.cumulative.2030.woco2 = p.default.woco2[[4]],
  subsidies.cumulative.2030.wco2 = p.default.wco2[[4]],
  ymax = p.conservative.woco2[[4]],
  title = "Required cumulative subsidies\nfor all projects by 2030 (central)")

# Progressive scenario
p.progressive.woco2 <- plotRequiredSubsidies(
  scenario.select = "Progressive",
  carbonprice.select = FALSE,
  title = "Required subsidies for all projects by 2030\n(progressive)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

p.progressive.wco2 <- plotRequiredSubsidies(
  scenario.select = "Progressive",
  carbonprice.select = TRUE,
  title = "Required subsidies for all projects by 2030\n(progressive)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

p.progressive.subsidygap <- plotSubsidyGap(
  scenario.select = "Progressive",
  subsidies.cumulative.2030.woco2 = p.progressive.woco2[[4]],
  subsidies.cumulative.2030.wco2 = p.progressive.wco2[[4]],
  ymax = p.conservative.woco2[[4]],
  title = "Required cumulative subsidies\nfor all projects by 2030 (progressive)")

# Arrange
p.panels <- plot_grid(
  p.default.woco2[[1]] + theme(legend.position = "none"),
  p.default.wco2[[1]] + theme(legend.position = "none"),
  p.default.subsidygap,
  p.progressive.woco2[[1]] + theme(legend.position = "none"),
  p.progressive.wco2[[1]] + theme(legend.position = "none"),
  p.progressive.subsidygap,
  p.conservative.woco2[[1]] + theme(legend.position = "none"),
  p.conservative.wco2[[1]] + theme(legend.position = "none"),
  p.conservative.subsidygap,
  align = "hv",
  ncol = 3,
  rel_widths = c(1,1,0.9),
  labels = "auto",
  label_size = font.size)

p.legend <- get_legend(p.default.woco2[[1]] + theme(legend.box = "horizontal",
                                               legend.direction = "vertical",
                                               legend.justification = "centre"))

p.row.legend <- plot_grid(
  p.legend,
  NULL,
  ncol = 2,
  rel_widths = c(2, 0.9)
)

# Add labels at the very top
p.woco2 <- ggdraw() +
  draw_label(aes(label = "Without carbon price"),
             size = 6,
             fontface = "bold")

p.wco2 <- ggdraw() +
  draw_label(aes(label = "With ambitious carbon price"),
             size = 6,
             fontface = "bold")
  
p.colnames <- plot_grid(
  p.woco2,
  p.wco2,
  NULL,
  ncol = 3,
  rel_widths = c(1, 1, 0.9)
)

# Total plot
p <- plot_grid(
  p.colnames,
  p.panels,
  p.row.legend,
  ncol = 1,
  rel_heights = c(0.03, 1, 0.07)
)

ggsave("figures/EXT_DATA_FIGURE06_subsidies_scenarios.png",
       width = 180, height = 180, units = "mm", bg = "white", dpi = 600)
ggsave("figures/EXT_DATA_FIGURE06_subsidies_scenarios.svg",
       width = 180, height = 180, units = "mm", bg = "white")

```

## EXTENDED DATA FIGURE 7: Required subsidies compared with solar PV and wind support

```{r extended data figure 7}

data.iea.support <- read_excel("data/iea_support_wind_solar.xlsx") %>% 
  pivot_longer(`2010`:`2030`, names_to = "year", values_to = "value") %>% 
  mutate(year = as.numeric(year)) %>% 
  select(-comment)

# Total support
total.vars <- c(
  "Support|Wind onshore|LCOE|Total" = "Wind",
  "Support|Wind offshore|LCOE|Total" = "Wind",
  "Support|PV res|LCOE|Total" = "Solar PV",
  "Support|PV utility|LCOE|Total" = "Solar PV",
  "Support|PV comm|LCOE|Total" = "Solar PV"
)

total.colors <- c(
  "Wind" = "#3E7AD3",
  "Solar PV" = "#F1A800"
)

# Total support for wind and PV according to IEA
data.plot.iea.support.total <- data.iea.support %>% 
  filter(variable %in% names(total.vars)) %>% 
  revalue.levels(variable = total.vars) %>% 
  group_by(source, variable, region, unit, year) %>% 
  summarise(value = sum(value))

data.iea.shares.2015 <- data.plot.iea.support.total %>%
  filter(year == 2015) %>%
  ungroup() %>%
  mutate(share = value/sum(value)) %>% 
  select(source, region, unit, variable, share)

data.plot.iea.support.total.before2015 <- data.iea.support %>% 
  filter(variable == "Support|Wind+PV|LCOE|Total",
         !is.na(value)) %>% 
  select(-variable) %>% 
  # Distribute to PV and Wind according to share in 2015
  full_join(data.iea.shares.2015, relationship = "many-to-many") %>% 
  mutate(value = value * share)

# Total required support for hydrogen
data.plot.h2.subsidies <- data.subsidies.annual %>% 
  mutate(group = ifelse(greenproduct == "Green hydrogen", "Green hydrogen", "Electrofuels")) %>% 
  group_by(scenario, carbonprice, group, year) %>% 
  summarise(subsidies.annual = sum(subsidies.annual))

data.plot.h2.subsidies.centre <- data.plot.h2.subsidies %>% 
  filter(scenario == "Default")

data.plot.h2.subsidies.range <- data.plot.h2.subsidies %>% 
  filter(scenario %in% c("Progressive", "Conservative")) %>% 
  pivot_wider(names_from = scenario, values_from = subsidies.annual)

ymin <- data.plot.iea.support.total %>% filter(!is.na(value)) %>% pull(value) %>% min()
ymax <- data.plot.h2.subsidies.range %>% filter(carbonprice == FALSE) %>% pull(Conservative) %>% max()

colors.group <- c(
  "Green hydrogen" = "#76AF63",
  "Electrofuels" = "#8E5FBE"
)

p.total.woco2 <- ggplot() +
  # Line at zero
  geom_hline(aes(yintercept = 0),
             linewidth = 0.25,
             color = "grey") +
  # Support for wind and solar
  geom_line(data = data.plot.iea.support.total,
            mapping = aes(x = year, y = value, color = variable)) +
  scale_color_manual(name = "Historical and projected support",
                     values = total.colors) +
  # Approximated support for wind and solar before 2015
  geom_line(data = data.plot.iea.support.total.before2015,
            mapping = aes(x = year, y = value, color = variable),
            linetype = "dotted") +
  new_scale_color() +
  # Support for green hydrogen
  geom_ribbon(data = data.plot.h2.subsidies.range %>% filter(carbonprice == FALSE),
              mapping = aes(x = year, ymin = Progressive, ymax = Conservative, fill = group),
              alpha = 0.2) +
  geom_line(data = data.plot.h2.subsidies.centre %>% filter(carbonprice == FALSE),
            mapping = aes(x = year, y = subsidies.annual, color = group)) +
  scale_fill_manual(name = "Required support",
                   values = colors.group,
                   guide = guide_legend(order = 1)) +
  scale_color_manual(name = "Required support",
                     values = colors.group,
                     guide = guide_legend(order = 1)) +
  scale_x_continuous(name = "Year",
                     limits = c(2010, 2050)) +
  scale_y_continuous(name = "Billion $/yr",
                     limits = c(ymin, ymax)) +
  ggtitle("Total annual support without carbon price")

p.total.wco2 <- ggplot() +
  # Line at zero
  geom_hline(aes(yintercept = 0),
             linewidth = 0.25,
             color = "grey") +
  # Support for wind and solar
  geom_line(data = data.plot.iea.support.total,
            mapping = aes(x = year, y = value, color = variable)) +
  scale_color_manual(name = "Historical and projected support",
                     values = total.colors) +
  # Approximated support for wind and solar before 2015
  geom_line(data = data.plot.iea.support.total.before2015,
            mapping = aes(x = year, y = value, color = variable),
            linetype = "dotted") +
  new_scale_color() +
  # Support for green hydrogen
  geom_ribbon(data = data.plot.h2.subsidies.range %>% filter(carbonprice == TRUE),
              mapping = aes(x = year, ymin = Progressive, ymax = Conservative, fill = group),
              alpha = 0.2) +
  geom_line(data = data.plot.h2.subsidies.centre %>% filter(carbonprice == TRUE),
            mapping = aes(x = year, y = subsidies.annual, color = group)) +
  scale_fill_manual(name = "Required support",
                   values = colors.group,
                   guide = guide_legend(order = 1)) +
  scale_color_manual(name = "Required support",
                     values = colors.group,
                     guide = guide_legend(order = 1)) +
  scale_x_continuous(name = "Year",
                     limits = c(2010, 2050)) +
  scale_y_continuous(name = "Billion $/yr",
                     limits = c(ymin, ymax)) +
  ggtitle("Total annual support with ambitious carbon price")

# Relative support
relative.vars <- c(
  "Support|Wind onshore|LCOE|Relative" = "Wind onshore",
  "Support|Wind offshore|LCOE|Relative" = "Wind offshore",
  "Support|PV res|LCOE|Relative" = "PV residential",
  "Support|PV utility|LCOE|Relative" = "PV utility",
  "Support|PV comm|LCOE|Relative" = "PV commercial"
)

relative.colors <- c(
  "Wind onshore" = "#49D3FF",
  "Wind offshore" = "#3E7AD3",
  "PV residential" = "#E34946",
  "PV utility" = "#FDD324",
  "PV commercial" = "#F1A800"
)

# Relative support for wind and PV according to IEA
data.plot.iea.relative <- data.iea.support %>% 
  filter(variable %in% names(relative.vars)) %>% 
  revalue.levels(variable = relative.vars)

# Relative required support for hydrogen
data.plot.h2.relsubsidies <- data.subsidies.annual %>% 
  rowwise() %>% 
  mutate(prod = max(cap - cap.backed, 0) * flh.electrolysis * eta) %>% 
  mutate(group = ifelse(greenproduct == "Green hydrogen", "Green hydrogen", "Electrofuels")) %>% 
  group_by(scenario, carbonprice, group, year) %>% 
  summarise(subsidies.relative = 1E9*sum(subsidies.annual) / sum(prod))

data.plot.h2.relsubsidies.centre <- data.plot.h2.relsubsidies %>% 
  filter(scenario == "Default")

data.plot.h2.relsubsidies.range <- data.plot.h2.relsubsidies %>% 
  filter(scenario %in% c("Progressive", "Conservative")) %>% 
  pivot_wider(names_from = scenario, values_from = subsidies.relative)

p.relative.woco2 <- ggplot() +
  # Line at zero
  geom_hline(aes(yintercept = 0),
             linewidth = 0.25,
             color = "grey") +
  # Support for wind and solar
  geom_line(data = data.plot.iea.relative,
            mapping = aes(x = year, y = value, color = variable)) +
  scale_color_manual(name = "Historical and projected support",
                     values = relative.colors,
                     guide = guide_legend(order = 2)) +
  new_scale_color() +
  # Support for green hydrogen
  geom_ribbon(data = data.plot.h2.relsubsidies.range %>% filter(carbonprice == FALSE),
              mapping = aes(x = year, ymin = Progressive, ymax = Conservative, fill = group),
              alpha = 0.2) +
  geom_line(data = data.plot.h2.relsubsidies.centre %>% filter(carbonprice == FALSE),
            mapping = aes(x = year, y = subsidies.relative, color = group)) +
  scale_fill_manual(name = "Required support",
                   values = colors.group,
                   guide = guide_legend(order = 1)) +
  scale_color_manual(name = "Required support",
                     values = colors.group,
                   guide = guide_legend(order = 1)) +
  scale_x_continuous(name = "Year",
                     limits = c(2010, 2050)) +
  ylab("$/MWh") +
  ggtitle("Relative support without carbon price")

p.relative.wco2 <- ggplot() +
  # Line at zero
  geom_hline(aes(yintercept = 0),
             linewidth = 0.25,
             color = "grey") +
  # Support for wind and solar
  geom_line(data = data.plot.iea.relative,
            mapping = aes(x = year, y = value, color = variable)) +
  scale_color_manual(name = "Historical and projected support",
                     values = relative.colors,
                     guide = guide_legend(order = 2)) +
  new_scale_color() +
  # Support for green hydrogen
  geom_ribbon(data = data.plot.h2.relsubsidies.range %>% filter(carbonprice == TRUE),
              mapping = aes(x = year, ymin = Progressive, ymax = Conservative, fill = group),
              alpha = 0.2) +
  geom_line(data = data.plot.h2.relsubsidies.centre %>% filter(carbonprice == TRUE),
            mapping = aes(x = year, y = subsidies.relative, color = group)) +
  scale_fill_manual(name = "Required support",
                   values = colors.group,
                   guide = guide_legend(order = 1)) +
  scale_color_manual(name = "Required support",
                     values = colors.group,
                   guide = guide_legend(order = 1)) +
  scale_x_continuous(name = "Year",
                     limits = c(2010, 2050)) +
  ylab("$/MWh") +
  ggtitle("Relative support with ambitious carbon price")

# Arrange
p.panels <- plot_grid(
  p.total.woco2 + theme(legend.position = "none"),
  p.total.wco2 + theme(legend.position = "none"),
  p.relative.woco2 + theme(legend.position = "none"),
  p.relative.wco2 + theme(legend.position = "none"),
  ncol = 2,
  labels = "auto",
  label_size = font.size
)

leg1 <- get_legend(p.total.woco2)

leg2 <- get_legend(p.relative.woco2)

legend <- plot_grid(
  leg1,
  leg2,
  ncol = 1
)

p <- plot_grid(
  p.panels,
  legend,
  ncol = 2,
  rel_widths = c(1, 0.22)
)

ggsave("figures/EXT_DATA_FIGURE07_required_subsidies_comparison.png",
       width = 180, height = 120, units = "mm", bg = "white", dpi = 600)
ggsave("figures/EXT_DATA_FIGURE07_required_subsidies_comparison.svg",
       width = 180, height = 120, units = "mm", bg = "white", dpi = 600)

```

# SUPPLEMENTARY FIGURES

## SUPPLEMENTARY FIGURE 1: IEA data quality validation (overview)

```{r supplementary figure 1}

range.left <- seq(2020, 2024)
range.right <- seq(2024, 2030)

# Get installed capacity from IEA v2023 and outcome of 2023 projects
data.real <- data.v2023 %>% 
  filter(status == "Operational",
         year %in% seq(2020, 2022)) %>% 
  group_by(year) %>% 
  summarise(value = sum(cumcap.sum)) %>% 
  # Add 2023 value (cumulative sum doesn't work here as this only contains 2023 projects)
  add_row(year = 2023,
          value =
            filter(data.real, year == 2022)$value + 
            data.v2023outcome %>% filter(year == 2023) %>% summarise(value = sum(cap.sum)) %>% pull())

# Announcements in 2021 (IEA v2021)
data.v2021.plot.left <- data.v2021 %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

data.v2021.orig.plot.left <- data.v2021.orig %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Announcements in 2022 (IEA v2022)
data.v2022.plot.left <- data.v2022 %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

data.v2022.orig.plot.left <- data.v2022.orig %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Announcements in 2023 (IEA v2023)
data.v2023.plot.left <- data.v2023 %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

data.v2023.orig.plot.left <- data.v2023.orig %>% 
  filter(year %in% range.left) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

legend.labels <- c("v2021" = "Announcements by 2021",
                   "v2022" = "Announcements by 2022",
                   "v2023" = "Announcements by 2023")

# Plot until 2024
p.left <- ggplot() +
  # Project announcements
  geom_line(data = data.v2021.plot.left,
            mapping = aes(x = year, y = sum, color = "v2021", alpha = "v2021", linetype = "Validated")) +
  geom_line(data = data.v2021.orig.plot.left,
            mapping = aes(x = year, y = sum, color = "v2021", alpha = "v2021", linetype = "Original")) +  
  geom_line(data = data.v2022.plot.left,
            mapping = aes(x = year, y = sum, color = "v2022", alpha = "v2022", linetype = "Validated")) +
  geom_line(data = data.v2022.orig.plot.left,
            mapping = aes(x = year, y = sum, color = "v2022", alpha = "v2022", linetype = "Original")) +  
  geom_line(data = data.v2023.plot.left,
            mapping = aes(x = year, y = sum, color = "v2023", alpha = "v2023", linetype = "Validated")) +
  geom_line(data = data.v2023.orig.plot.left,
            mapping = aes(x = year, y = sum, color = "v2023", alpha = "v2023", linetype = "Original")) +  
  scale_color_npg(name = "Project announcements",
                  guide = guide_legend(order = 1),
                  breaks = c("v2021",
                             "v2022",
                             "v2023"),
                  labels = legend.labels) +
  scale_alpha_manual(name = "Project announcements",
                     guide = guide_legend(order = 1),
                     breaks = c("v2021",
                                "v2022",
                                "v2023"),
                     values = c("v2021" = 0.4,
                                "v2022" = 0.7,
                                "v2023" = 1),
                     labels = legend.labels) +
  scale_linetype_manual(name = "Database",
                        guide = guide_legend(order = 2),
                        values = c("Validated" = "solid",
                                   "Original" = "dotted")) +
  new_scale_color() +
  # Realised projects
  geom_line(data = data.real,
            mapping = aes(x = year, y = value, color = "Realised projects")) +
  geom_point(data = data.real,
             mapping = aes(x = year, y = value, color = "Realised projects"),
             size = 1) +
  scale_color_manual(name = NULL, values = c("Realised projects" = "black"),
                     guide = guide_legend(order = 3)) +
  # Formatting
  xlab("Year") +
  ylab("Capacity (GW)") +
  ggtitle("Green hydrogen projects until 2024") +
  theme(legend.position = "inside",
        legend.position.inside = c(0.1,0.9),
        legend.justification = c("left","top"))

# Expectation in 2021 (IEA v2021)
data.v2021.plot.right <- data.v2021 %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

data.v2021.orig.plot.right <- data.v2021.orig %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Expectation in 2022 (IEA v2022)
data.v2022.plot.right <- data.v2022 %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

data.v2022.orig.plot.right <- data.v2022.orig %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Expectation in 2023 (IEA v2023)
data.v2023.plot.right <- data.v2023 %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

data.v2023.orig.plot.right <- data.v2023.orig %>% 
  filter(year %in% range.right) %>% 
  group_by(year) %>% 
  summarise(sum = sum(cumcap.sum))

# Plot from 2024 to 2030
p.right <- ggplot() +
  # Project announcements
  geom_line(data = data.v2021.plot.right,
            mapping = aes(x = year, y = sum, color = "v2021", alpha = "v2021", linetype = "Validated")) +
  geom_line(data = data.v2021.orig.plot.right,
            mapping = aes(x = year, y = sum, color = "v2021", alpha = "v2021", linetype = "Original")) +  
  geom_line(data = data.v2022.plot.right,
            mapping = aes(x = year, y = sum, color = "v2022", alpha = "v2022", linetype = "Validated")) +
  geom_line(data = data.v2022.orig.plot.right,
            mapping = aes(x = year, y = sum, color = "v2022", alpha = "v2022", linetype = "Original")) +  
  geom_line(data = data.v2023.plot.right,
            mapping = aes(x = year, y = sum, color = "v2023", alpha = "v2023", linetype = "Validated")) +
  geom_line(data = data.v2023.orig.plot.right,
            mapping = aes(x = year, y = sum, color = "v2023", alpha = "v2023", linetype = "Original")) +  
  scale_color_npg(name = "Project announcements",
                  guide = guide_legend(order = 1),
                  breaks = c("v2021",
                             "v2022",
                             "v2023"),
                  labels = legend.labels) +
  scale_alpha_manual(name = "Project announcements",
                     guide = guide_legend(order = 1),
                     breaks = c("v2021",
                                "v2022",
                                "v2023"),
                     values = c("v2021" = 0.4,
                                "v2022" = 0.7,
                                "v2023" = 1),
                     labels = legend.labels) +
  scale_linetype_manual(name = "Database",
                        guide = guide_legend(order = 2),
                        values = c("Validated" = "solid",
                                   "Original" = "dotted")) +
  # Formatting
  xlab("Year") +
  scale_x_continuous(breaks = seq(2024,2030)) +
  ylab("Capacity (GW)") +
  ggtitle("Green hydrogen projects from 2024") +
  theme(legend.position = "inside",
        legend.position.inside = c(0.1,0.9),
        legend.justification = c("left","top"))

p <- plot_grid(
  p.left,
  p.right,
  ncol = 2,
  labels = "auto",
  label_size = font.size,
  rel_widths = c(1,1.45))

print(p)

ggsave("figures/SUPP_FIGURE01_database_quality_validation_overview.png",
       width = 180, height = 70, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE01_database_quality_validation_overview.svg",
       width = 180, height = 70, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 2: IEA data quality validation (2021 project detail)

```{r supplementary figure 2}

p.status <- plotDatabaseDelta(data.new.proj = data.v2021.proj,
                              data.old.proj = data.v2021.orig.proj,
                              version = "2021",
                              by = "status")

p.region <- plotDatabaseDelta(data.new.proj = data.v2021.proj,
                              data.old.proj = data.v2021.orig.proj,
                              version = "2021",
                              by = "region")

leg.status <- get_legend(p.status[[1]])

leg.region <- get_legend(p.region[[1]])

p.panels <- plot_grid(
  p.status[[1]] + theme(legend.position = "none"),
  p.status[[2]] + theme(legend.position = "none"),
  leg.status,
  p.region[[1]] + theme(legend.position = "none"),
  p.region[[2]] + theme(legend.position = "none"),
  leg.region,
  ncol = 3,
  labels = c("a", "b", "", "c", "d", ""),
  label_size = font.size,
  rel_widths = c(5,7,1.5))

toplabel <- ggdraw() +
  draw_label(
    x = 0.5,
    y = 0.8,
    "Data quality validation of IEA Hydrogen Projects Database 2021",
    hjust = 0.5,
    vjust = 1,
    size = font.size)

p <- plot_grid(toplabel,
               p.panels,
               ncol = 1,
               rel_heights = c(0.08, 1))

ggsave("figures/SUPP_FIGURE02_data_quality_validation_v2021.png",
       width = 180, height = 110, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE02_data_quality_validation_v2021.svg",
       width = 180, height = 110, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 3: IEA data quality validation (2022 project detail)

```{r supplementary figure 3}

p.status <- plotDatabaseDelta(data.new.proj = data.v2022.proj,
                              data.old.proj = data.v2022.orig.proj,
                              version = "2022",
                              by = "status")

p.region <- plotDatabaseDelta(data.new.proj = data.v2022.proj,
                              data.old.proj = data.v2022.orig.proj,
                              version = "2022",
                              by = "region")

leg.status <- get_legend(p.status[[1]])

leg.region <- get_legend(p.region[[1]])

p.panels <- plot_grid(
  p.status[[1]] + theme(legend.position = "none"),
  p.status[[2]] + theme(legend.position = "none"),
  leg.status,
  p.region[[1]] + theme(legend.position = "none"),
  p.region[[2]] + theme(legend.position = "none"),
  leg.region,
  ncol = 3,
  labels = c("a", "b", "", "c", "d", ""),
  label_size = font.size,
  rel_widths = c(5,7,1.5))

toplabel <- ggdraw() +
  draw_label(
    x = 0.5,
    y = 0.8,
    "Data quality validation of IEA Hydrogen Projects Database 2022",
    hjust = 0.5,
    vjust = 1,
    size = font.size)

p <- plot_grid(toplabel,
               p.panels,
               ncol = 1,
               rel_heights = c(0.08, 1))

ggsave("figures/SUPP_FIGURE03_data_quality_validation_v2022.png",
       width = 180, height = 110, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE03_data_quality_validation_v2022.svg",
       width = 180, height = 110, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 4: IEA data quality validation (2023 project detail)

```{r supplementary figure 4}

p.status <- plotDatabaseDelta(data.new.proj = data.v2023.proj,
                              data.old.proj = data.v2023.orig.proj,
                              version = "2023",
                              by = "status")

p.region <- plotDatabaseDelta(data.new.proj = data.v2023.proj,
                              data.old.proj = data.v2023.orig.proj,
                              version = "2023",
                              by = "region")

leg.status <- get_legend(p.status[[1]])

leg.region <- get_legend(p.region[[1]])

p.panels <- plot_grid(
  p.status[[1]] + theme(legend.position = "none"),
  p.status[[2]] + theme(legend.position = "none"),
  leg.status,
  p.region[[1]] + theme(legend.position = "none"),
  p.region[[2]] + theme(legend.position = "none"),
  leg.region,
  ncol = 3,
  labels = c("a", "b", "", "c", "d", ""),
  label_size = font.size,
  rel_widths = c(5,7,1.5))

toplabel <- ggdraw() +
  draw_label(
    x = 0.5,
    y = 0.8,
    "Data quality validation of IEA Hydrogen Projects Database 2023",
    hjust = 0.5,
    vjust = 1,
    size = font.size)

p <- plot_grid(toplabel,
               p.panels,
               ncol = 1,
               rel_heights = c(0.08, 1))

ggsave("figures/SUPP_FIGURE04_data_quality_validation_v2023.png",
       width = 180, height = 110, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE04_data_quality_validation_v2023.svg",
       width = 180, height = 110, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 5: 2022 implementation gap

```{r supplementary figure 5}

# Calculate project tracking for 2022
data.track.calc2022 <-
  calcProjectTracking(
    data.v2021proj = data.v2021.proj,
    data.v2022proj = data.v2022.proj,
    data.v2023proj = data.v2023.proj,
    year.tracking = 2022
  )

# Calculate rates of success, delay and disappearance for 2022
# Comparing announcements in 2021 to outcome
data.rates.calc2022.v21v23 <-
  calcRates(
    data.track.calc = data.track.calc2022,
    year.tracking = 2022,
    compare = "v2021-v2023"
  )

# Print rates
data.rates.calc2022.v21v23 %>% 
  group_by(outcome) %>% 
  summarise(rate.tot = round(sum(cap.share.tot), 3))

# Comparing announcements in 2022 to outcome
data.rates.calc2022.v22v23 <-
  calcRates(
    data.track.calc = data.track.calc2022,
    year.tracking = 2022,
    compare = "v2022-v2023"
  )

# Print rates
data.rates.calc2022.v22v23 %>% 
  group_by(outcome) %>% 
  summarise(rate.tot = round(sum(cap.share.tot), 3))

# Plot sankey diagram
p.track.2022 <-
  plotProjectTracking(
    data.track.calc = data.track.calc2022,
    shift.incoming = 4000,
    shift.outgoing = 500,
    y.upper = 3000,
    year.tracking = 2022,
    annotate.labels = TRUE,
    annotate.gap = TRUE
  )

# Plot rates of success, delay and disappearance (v2021 vs outcome)
temp <- plotRates(
  data.rates = data.rates.calc2022.v21v23,
  year.tracking = 2022
  )

p.rates.2022.v21v23.tot <- temp[[1]]
p.rates.2022.v21v23.tot.range <- temp[[2]]
p.rates.2022.v21v23.disagg <- temp[[3]]
p.rates.2022.v21v23.disagg.range <- temp[[4]]

# Plot rates of success, delay and disappearance (v2022 vs outcome)
temp <- plotRates(
  data.rates = data.rates.calc2022.v22v23,
  year.tracking = 2022)

p.rates.2022.v22v23.tot <- temp[[1]]
p.rates.2022.v22v23.tot.range <- temp[[2]]
p.rates.2022.v22v23.disagg <- temp[[3]]
p.rates.2022.v22v23.disagg.range <- temp[[4]]

# Left-align panels
plot.align <- align_plots(p.track.2022 + theme(legend.position = "none"),
                          p.rates.2022.v21v23.tot + theme(legend.position = "none"),
                          align = "v",
                          axis = "l")

# Build bottom row left
# Horizontally align panels
plot.align.left <- align_plots(plot.align[[2]],
                               p.rates.2022.v21v23.disagg + theme(legend.position = "none"),
                               align = "h")

bottom_row.left.panels.tot <- plot_grid(
  plot.align.left[[1]],
  p.rates.2022.v21v23.tot.range,
  ncol = 1,
  rel_heights = c(1, 0.15),
  align = "v"
)

bottom_row.left.panels.disagg <- plot_grid(
  plot.align.left[[2]],
  p.rates.2022.v21v23.disagg.range,
  ncol = 1,
  rel_heights = c(1, 0.15),
  align = "v"
)

bottom_row.left.panels <- plot_grid(
  bottom_row.left.panels.tot,
  bottom_row.left.panels.disagg,
  ncol = 2,
  align = "h"
)

bottom_row.left.title <- ggdraw() +
  draw_label(
    "Announcements by 2021 vs. outcome in 2023",
    x = 0.55,
    hjust = 0.5,
    vjust = 0,
    size = font.size
  )

bottom_row.left <- plot_grid(
  bottom_row.left.title,
  bottom_row.left.panels,
  ncol = 1,
  rel_heights = c(0.15, 2)
)

# Build bottom row right
# Horizontally align panels
plot.align.right <- align_plots(p.rates.2022.v22v23.tot + theme(legend.position = "none"),
                                p.rates.2022.v22v23.disagg + theme(legend.position = "none"),
                                align = "h")

bottom_row.right.panels.tot <- plot_grid(
  plot.align.right[[1]],
  p.rates.2022.v22v23.tot.range,
  ncol = 1,
  rel_heights = c(1, 0.15),
  align = "v"
)

bottom_row.right.panels.disagg <- plot_grid(
  plot.align.right[[2]],
  p.rates.2022.v22v23.disagg.range,
  ncol = 1,
  rel_heights = c(1, 0.15),
  align = "v"
)

bottom_row.right.panels <- plot_grid(
  bottom_row.right.panels.tot,
  bottom_row.right.panels.disagg,
  ncol = 2,
  align = "h")

bottom_row.right.title <- ggdraw() +
  draw_label(
    "Announcements by 2022 vs. outcome in 2023",
    x = 0.55,
    hjust = 0.5,
    vjust = 0,
    size = font.size)

bottom_row.right <- plot_grid(
  bottom_row.right.title,
  bottom_row.right.panels,
  ncol = 1,
  rel_heights = c(0.15,2)
)

# Complete bottom row
bottom_row <- plot_grid(
  bottom_row.left,
  bottom_row.right,
  NULL,
  get_legend(p.rates.2022.v21v23.tot),
  ncol = 4,
  rel_widths = c(1, 1, 0.02, 0.2),
  labels = c("b", "c", "", ""),
  label_size = font.size
)

# Build top row
top_row <- plot_grid(
  plot.align[[1]],
  NULL,
  ncol = 2,
  rel_widths = c(1, 0),
  labels = "a",
  label_size = font.size)

# Build entire plot
p <- plot_grid(top_row,
               NULL,
               bottom_row,
               ncol = 1,
               rel_heights = c(1.3,0.2,1))

print(p)

ggsave("figures/SUPP_FIGURE05_2022_implementation_gap.png",
       width = 180, height = 150, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE05_2022_implementation_gap.svg",
       width = 180, height = 150, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 6: Project tracking for 2024

```{r supplementary figure 6}

# Calculate project tracking for 2022
data.track.calc2024 <-
  calcProjectTracking(
    data.v2021proj = data.v2021.proj,
    data.v2022proj = data.v2022.proj,
    data.v2023proj = data.v2023.proj,
    year.tracking = 2024
  )

# Plot sankey diagram
p.track.2024 <-
  plotProjectTracking(
    data.track.calc = data.track.calc2024,
    shift.incoming = 23000,
    shift.outgoing = 1500,
    y.upper = 13000,
    year.tracking = 2024,
    annotate.labels = TRUE,
    annotate.gap = FALSE
  )

ggsave("figures/SUPP_FIGURE06_project_tracking_2024.png",
       width = 180, height = 100, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE06_project_tracking_2024.svg",
       width = 180, height = 100, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 7-10: Regional 2023 implementation gap

```{r supplementary figure 7-10}

# Projects in Asia
p.gap2023.asia <- plotImplementationGap2023(
  region.plot = "Asia",
  title = "Tracking green hydrogen projects announced for 2023 in Asia",
  shift.incoming = 2000,
  shift.outgoing = 300,
  y.upper = 1000,
  legend = "v22",
  annotate.gap = FALSE)

print(p.gap2023.asia)

ggsave("figures/SUPP_FIGURE07_implementation_gap_2023_Asia.png",
       width = 180, height = 150, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE07_implementation_gap_2023_Asia.svg",
       width = 180, height = 150, units = "mm", bg = "white")


# Projects in Australia
p.gap2023.australia <- plotImplementationGap2023(
  region.plot = "Australia",
  title = "Tracking green hydrogen projects announced for 2023 in Australia",
  shift.incoming = 1800,
  shift.outgoing = 250,
  y.upper = 1000,
  legend = "v23",
  annotate.gap = FALSE)

print(p.gap2023.australia)

ggsave("figures/SUPP_FIGURE08_implementation_gap_2023_Australia.png",
       width = 180, height = 150, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE08_implementation_gap_2023_Australia.svg",
       width = 180, height = 150, units = "mm", bg = "white")


# Projects in Europe
p.gap2023.europe <- plotImplementationGap2023(
  region.plot = "Europe",
  title = "Tracking green hydrogen projects announced for 2023 in Europe",
  shift.incoming = 4000,
  shift.outgoing = 500,
  y.upper = 2000,
  legend = "v22",
  annotate.gap = FALSE)

print(p.gap2023.europe)

ggsave("figures/SUPP_FIGURE09_implementation_gap_2023_Europe.png",
       width = 180, height = 150, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE09_implementation_gap_2023_Europe.svg",
       width = 180, height = 150, units = "mm", bg = "white")


# Projects in N. America
p.gap2023.namerica <- plotImplementationGap2023(
  region.plot = "N. America",
  title = "Tracking green hydrogen projects announced for 2023 in North America",
  shift.incoming = 2300,
  shift.outgoing = 300,
  y.upper = 1000,
  legend = "v23",
  annotate.gap = FALSE)

print(p.gap2023.namerica)

ggsave("figures/SUPP_FIGURE10_implementation_gap_2023_N_America.png",
       width = 180, height = 150, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE10_implementation_gap_2023_N_America.svg",
       width = 180, height = 150, units = "mm", bg = "white")

```

## SUPPLEMENRATY FIGURE 11: Project announcements and median 1.5°C scenario until 2050

```{r supplementary figure 11}

# Project announcements
data.plot.projects <- data.iea.enduse %>% 
  group_by(status, year) %>% 
  filter(year %in% seq(2020, 2050)) %>% 
  summarise(cap.sum = sum(cap.sum),
            cumcap.sum = sum(cumcap.sum))

# Cumulative project announcements in 2030
projects.cum.2030 <- data.plot.projects %>%
  filter(year == 2030) %>%
  ungroup() %>%
  summarise(sum(cumcap.sum)) %>%
  pull()

# 1.5C scenario median
data.plot.scenarios <- data.subsidies.annual %>%
  filter(year.cap %in% seq(2031, 2050)) %>% 
  group_by(year.cap, enduse) %>% 
  summarise(cap = unique(cap)) %>% 
  group_by(year.cap) %>% 
  summarise(cap = 1E-3*sum(cap)) %>% 
  mutate(cumcap = cumsum(cap) + projects.cum.2030)

# Capacity secured by demand policies
data.plot.capdemand <- data.subsidies.annual %>% 
  filter(scenario == "Default",
         year.cap %in% seq(2024, 2030)) %>% 
  group_by(year.cap, enduse) %>% 
  summarise(cap = unique(cap),
            cap.backed = unique(cap.backed)) %>% 
  # Capacities only backed if enough announcements available, i.e. demand-side
  # quantities are not shifted around between end-uses
  rowwise() %>% 
  mutate(cap.backed = 1E-3*min(cap.backed, cap)) %>% 
  group_by(year.cap) %>% 
  summarise(cap.backed = sum(cap.backed)) %>% 
  complete(year.cap = seq(2024, 2050), fill = list(cap.backed = 0)) %>% 
  mutate(cumcap.backed = cumsum(cap.backed))

# Capacity additions
p.added <- ggplot() +
  geom_col(data = data.plot.scenarios,
           mapping = aes(x = year.cap, y = cap, fill = "Median"),
           alpha = 0.5) +
  scale_fill_manual(name = "1.5°C scenarios",
                    values = c("Median" = "grey"),
                    guide = guide_legend(order = 1)) +
  new_scale_fill() +
  geom_col(data = data.plot.projects,
           mapping = aes(x = year, y = cap.sum, fill = status)) +
  scale_fill_npg(name = "Project announcement status",
                 guide = guide_legend(order = 2)) +
  geom_col_pattern(data = data.plot.capdemand,
                   mapping = aes(x = year.cap, y = cap.backed, pattern = "Implemented policies"),
                   fill = NA,
                   pattern_fill = "black",
                   pattern_size = 0.1,
                   pattern_density = 0.2,
                   pattern_spacing = 0.01) +
  scale_pattern_manual(name = "Capacity supported by\ndemand-side regulation",
                       values = c("Implemented policies" = "stripe"),
                       guide = guide_legend(order = 3)) +
  scale_x_continuous(
    name = "Year",
    breaks = seq(2020, 2050, 5)
  ) + 
  ylab("Added capacity (GW/yr)") +
  ggtitle("Annual capacity additions")

p.cumulative <- ggplot() +
  geom_col(data = data.plot.scenarios,
           mapping = aes(x = year.cap, y = cumcap, fill = "Median"),
           alpha = 0.5) +
  scale_fill_manual(name = "1.5°C scenarios",
                    values = c("Median" = "grey")) +
  new_scale_fill() +
  geom_col(data = data.plot.projects,
           mapping = aes(x = year, y = cumcap.sum, fill = status)) +
  scale_fill_npg(name = "Project announcement status") +#
  geom_col_pattern(data = data.plot.capdemand,
                   mapping = aes(x = year.cap, y = cumcap.backed, pattern = "Implemented policies"),
                   fill = NA,
                   pattern_fill = "black",
                   pattern_size = 0.1,
                   pattern_density = 0.2,
                   pattern_spacing = 0.01) +
  scale_pattern_manual(name = "Capacity supported by\ndemand-side regulation",
                       values = c("Implemented policies" = "stripe"),
                       guide = guide_legend(order = 3)) +
  scale_x_continuous(
    name = "Year",
    breaks = seq(2020,2050,5)
  ) + 
  ylab("Capacity (GW)") +
  ggtitle("Cumulative installed capacity")

p <- plot_grid(
  p.added + theme(legend.position = "none"),
  p.cumulative + theme(legend.position = "none"),
  get_legend(p.added),
  labels = c("a", "b"),
  label_size = font.size,
  ncol = 3,
  rel_widths = c(1,1,0.45)
)

print(p)

ggsave("figures/SUPP_FIGURE11_project_announcements_1.5C_median.png",
       width = 160, height = 60, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE11_project_announcements_1.5C_median.svg",
       width = 160, height = 60, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 12: Technological learning for electrolysis

```{r supplementary figure 12}

data.plot <- data.costs %>% 
  filter(variable %in% c("Tech|Green hydrogen|CAPEX|Stack", "Tech|Green hydrogen|CAPEX|Other")) %>% 
  order.levels(variable = c("Tech|Green hydrogen|CAPEX|Stack", "Tech|Green hydrogen|CAPEX|Other"),
               scenario = c("Default", "Progressive", "Conservative"))

colors.capex <- c(
    "Tech|Green hydrogen|CAPEX|Stack" = "#8F9F7C",
    "Tech|Green hydrogen|CAPEX|Other" = "#4B5C67")

labels.capex <- c(
    "Tech|Green hydrogen|CAPEX|Stack" = "Stack",
    "Tech|Green hydrogen|CAPEX|Other" = "Balance of plant")

p <- ggplot() +
  geom_col(data = data.plot,
           mapping = aes(x = period, y = value/1E3, fill = variable)) +
  scale_fill_manual(
    name = "Components",
    values = colors.capex,
    labels = labels.capex) +
  scale_x_continuous(name = "Year", breaks = c(2023, seq(2030, 2050, 5))) +
  scale_y_continuous(name = expression("Specific investment costs ($/kW"[el]*")"),
                     breaks = seq(0, 2500, 500)) +
  ggtitle("Technological learning of electrolyser investment costs") +
  facet_wrap(~scenario, labeller = labeller(scenario = c("Default" = "Central"))) +
  theme(strip.background =element_rect(fill = NA),
        axis.line = element_line(linewidth = 0),
        panel.border = element_rect(color = "black")) + 
  background_grid(major = c("y"), size.major = 0.25)

print(p)

ggsave("figures/SUPP_FIGURE12_technological_learning.png",
       width = 160, height = 60, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE12_technological_learning.svg",
       width = 160, height = 60, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 13: Green hydrogen in 1.5°C scenario under IEA NZE end-use shares

```{r supplementary figure 13}

# Mapping of IEA NZE scenario end uses to IEA Hydrogen Projects Database end uses
mapping.nze.h2projects.plot <- tribble(
  ~ variable, ~ enduse,
  "Final Energy|Industry|Refining|Hydrogen", "refining",
  "Final Energy|Industry|Chemicals|Hydrogen", "chemicals", # Approximation
  "Final Energy|Industry|Steel|Hydrogen", "ironsteel",
  "Final Energy|Industry|Other|Hydrogen", "otherind",
  "Final Energy|Transportation|Aviation and Marine|Hydrogen", "synfuels",
  "Secondary Energy Input|Hydrogen|Electricity", "power",
  "Final Energy|Transportation|Road|Hydrogen", "mobility",
  "Final Energy|Other|Hydrogen", NA_character_
)

# IEA H2 database enduses
data.plot.h2db <- data.iea.enduse %>% 
  filter(year == 2030) %>% 
  mutate(enduse = factor(enduse, levels = enduse.order.both)) %>% 
  group_by(enduse) %>% 
  summarise(cumcap.sum = sum(cumcap.sum))

# Create tibble with cumulative capacity in 2030 according to project announcements
# and afterwards according to median of 1.5°C scenarios
data.cumcap <- tibble(year = 2030, cumcap = data.plot.h2db %>% summarise(sum(cumcap.sum)) %>% pull()) %>% 
  bind_rows(data.scenarios.median.after2030) %>% 
  select(year, cumcap) %>% 
  complete(year = seq(2030, 2050)) %>% 
  mutate(cumcap = na.approx(cumcap))

# IEA NZE enduses
data.plot.nze <- full_join(data.iea.nze, mapping.nze.h2projects.plot) %>% 
  mutate(enduse = replace_na(enduse, "other")) %>% 
  # Linear interpolation
  filter(year %in% seq(2030, 2050)) %>% 
  group_by(source, report, region, variable, unit, enduse) %>% 
  complete(year = seq(2030, 2050)) %>% 
  mutate(value = na.approx(value)) %>% 
  # Calculate share of enduse
  group_by(year) %>% 
  mutate(share = value / sum(value)) %>% 
  # Calculate hypothetical capacities if the end-use of project announcements
  # would perfectly match the IEA NZE scenario
  full_join(data.cumcap) %>% 
  mutate(value.nze = share * cumcap) %>% 
  # Order levels
  order.levels(enduse = c("chemicals", enduse.order.both, "other"))

enduse.colors.fix <- enduse.colors
names(enduse.colors.fix) <- names(enduse.labels[enduse.order.both])

# Labels for end-uses in hydrogen projects database
data.label.h2db <- data.plot.h2db %>% 
  mutate(enduse = as.character(enduse)) %>% 
  mutate(enduse = replace_na(enduse, "NA")) %>% 
  order.levels(enduse = c(enduse.order.both, "NA")) %>% 
  mutate(label = enduse) %>% 
  revalue.levels(label = enduse.labels)

# Labels for end-uses in IEA NZE
data.label.nze <- data.plot.nze %>% 
  filter(year == 2030) %>% 
  mutate(label = as.character(enduse)) %>% 
  revalue.levels(label = c(enduse.labels, "chemicals" = "Chemicals", "other" = "Other"))

# 2030 point
p.left <- ggplot() + 
  # Project end use data
  geom_col(data = data.plot.h2db,
           mapping = aes(x = 1, y = cumcap.sum, fill = enduse)) +
  scale_fill_manual(values = enduse.colors.fix,
                    labels = enduse.labels,
                    aesthetics = c("fill", "colour")) +
  geom_text_repel(data = data.label.h2db,
                  mapping = aes(x = 0.5, y = cumcap.sum, color = enduse, label = label),
                  position = position_stacknudge(vjust = 0.5, x = -2),
                  size = 5 / .pt,
                  hjust = 0,
                  direction = "y",
                  min.segment.length = 0,
                  segment.size = 0.25) +
  new_scale_fill() +
  new_scale_color() +
  # IEA NZE share based end use data
  geom_col(data = data.plot.nze %>% filter(year == 2030),
           mapping = aes(x = 2, y = value.nze, fill = enduse)) +
  scale_fill_manual(values = c(enduse.colors.fix, "chemicals" = "#66C2A5", "other" = "#7E7E7E"),
                    labels = c(enduse.labels, "chemicals" = "Chemicals"),
                    aesthetics = c("colour", "fill")) +
  geom_text_repel(data = data.label.nze,
                  mapping = aes(x = 2.5, y = value.nze, color = enduse, label = label),
                  position = position_stacknudge(vjust = 0.5, x = 1.2),
                  size = 5 / .pt,
                  hjust = 0,
                  direction = "y",
                  min.segment.length = 0,
                  segment.size = 0.25) +
  scale_x_continuous(name = NULL,
                     breaks = c(1, 2),
                     labels = c("Project\nannouncements", "IEA NZE\nshares")) +
  ylab("Capacity (GW)") +
  ggtitle("Green hydrogen end use in 2030") +
  theme(legend.position = "none")

# 2030-2050 IEA NZE shares
p.right <- ggplot() +
  geom_col(data = data.plot.nze,
           mapping = aes(x = year, y = value.nze, fill = enduse)) +
  scale_fill_manual(name = "End use",
                    values = c(enduse.colors.fix, "chemicals" = "#66C2A5", "other" = "#7E7E7E"),
                    labels = c(enduse.labels, "chemicals" = "Chemicals", "other" = "Other")) +
  scale_x_continuous(name = "Year") +
  ylab("Capacity (GW)") +
  ggtitle("Green hydrogen in median 1.5°C scenario with IEA NZE end-use shares")

p <- plot_grid(
  p.left,
  p.right,
  ncol = 2,
  labels = "auto",
  label_size = font.size,
  rel_widths = c(3, 5),
  align = "h")

ggsave("figures/SUPP_FIGURE13_iea_nze_enduse.png",
       width = 180, height = 80, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE13_iea_nze_enduse.svg",
       width = 180, height = 80, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 14: Number of end-uses of green hydrogen projects

```{r supplementary figure 14}

data.plot <- data.v2023.proj.enduse %>% 
  filter(year <= 2030,
         status != "Decommissioned") %>% 
  group_by(reference) %>% 
  mutate(enduse.unique = ifelse(all(is.na(enduse)), 0, n())) %>% 
  group_by(enduse.unique, status) %>% 
  summarise(capacity = sum(capacity))
  
p <- ggplot() +
  geom_col(data = data.plot,
           mapping = aes(x = enduse.unique, y = capacity/1E3, fill = status)) +
  scale_fill_npg(name = "Project status") +
  scale_x_continuous(name = "Number of designated end uses",
                     breaks = seq(0, 8),
                     labels = c("0 (NA)", seq(1,8))) +
  ylab("Capacity in 2030 (GW)") +
  ggtitle("Number of end uses of green hydrogen project announements by 2030")

print(p)

ggsave("figures/SUPP_FIGURE14_enduse_numbers.png",
       width = 100, height = 60, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE14_enduse_numbers.svg",
       width = 100, height = 60, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 15: End-use specific demand-side regulation

```{r supplementary figure 15}

# Capacity secured by demand policies
data.plot.demandbacked <- data.subsidies.annual %>% 
  filter(year.cap %in% seq(2024, 2030)) %>% 
  group_by(scenario, year.cap, enduse) %>% 
  summarise(cap.backed = 1E-3*unique(cap.backed)) %>% 
  group_by(scenario, enduse) %>% 
  filter(!all(cap.backed == 0)) %>% 
  group_by(scenario, enduse) %>% 
  arrange(year.cap) %>%
  mutate(cumcap.backed = cumsum(cap.backed)) %>% 
  order.levels(scenario = c("Default", "Progressive", "Conservative"))

data.plot.projects <- data.iea.enduse %>% 
  filter(year %in% seq(2024, 2030),
         enduse %in% unique(data.plot.demandbacked$enduse)) %>% 
  group_by(status, enduse, year) %>% 
  summarise(cumcap.sum = sum(cumcap.sum)) %>% 
  crossing(scenario = unique(data.plot.demandbacked$scenario)) %>% 
  order.levels(scenario = c("Default", "Progressive", "Conservative"))
    
p <- ggplot() +
  geom_col(data = data.plot.projects,
           mapping = aes(x = year, y = cumcap.sum, fill = status)) +
  scale_fill_npg(name = "Project status") +
  geom_col_pattern(data = data.plot.demandbacked,
                   mapping = aes(x = year.cap, y = cumcap.backed, pattern = "Implemented policies"),
                   fill = NA,
                   pattern_fill = "black",
                   pattern_size = 0.05,
                   pattern_density = 0.2,
                   pattern_spacing = 0.05) +
  scale_pattern_manual(name = "Capacity supported by\ndemand-side regulation",
                       values = c("Implemented policies" = "stripe")) +
  scale_x_continuous(name = "Year",
                     breaks = seq(2024, 2030, 2)) +
  ylab("Capacity (GW)") +
  ggtitle("Announced capacity and capacity backed by demand-side regulation") +
  facet_grid(scenario ~ enduse,
             labeller = labeller(enduse = enduse.labels,
                                 scenario = c("Default" = "Central"))) +
  theme(strip.background =element_rect(fill = NA),
        axis.line = element_line(linewidth = 0),
        panel.border = element_rect(color = "black")) 

print(p)

ggsave("figures/SUPP_FIGURE15_demand_regulation_by_enduse.png",
       width = 180, height = 120, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE15_demand_regulation_by_enduse.svg",
       width = 180, height = 120, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 16: Required subsidies for 1.5°C scenario median

```{r supplementary figure 16}

# Conservative scenario
p.conservative.woco2 <- plotRequiredSubsidies1.5Cmedian(
  scenario.select = "Conservative",
  carbonprice.select = FALSE,
  title = "Required subsidies for 1.5°C scenario median\n(conservative)")

p.conservative.wco2 <- plotRequiredSubsidies1.5Cmedian(
  scenario.select = "Conservative",
  carbonprice.select = TRUE,
  title = "Required subsidies for 1.5°C scenario median\n(conservative)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

# Default scenario
p.default.woco2 <- plotRequiredSubsidies1.5Cmedian(
  scenario.select = "Default",
  carbonprice.select = FALSE,
  title = "Required subsidies for 1.5°C scenario median\n(central)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

p.default.wco2 <- plotRequiredSubsidies1.5Cmedian(
  scenario.select = "Default",
  carbonprice.select = TRUE,
  title = "Required subsidies for 1.5°C scenario median\n(central)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

# Progressive scenario
p.progressive.woco2 <- plotRequiredSubsidies1.5Cmedian(
  scenario.select = "Progressive",
  carbonprice.select = FALSE,
  title = "Required subsidies for 1.5°C scenario median\n(progressive)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

p.progressive.wco2 <- plotRequiredSubsidies1.5Cmedian(
  scenario.select = "Progressive",
  carbonprice.select = TRUE,
  title = "Required subsidies for 1.5°C scenario median\n(progressive)",
  ymax.annual = p.conservative.woco2[[3]],
  secaxis.factor = p.conservative.woco2[[2]])

# Arrange panels
p.panels <- plot_grid(
  p.default.woco2[[1]] + theme(legend.position = "none"),
  p.default.wco2[[1]] + theme(legend.position = "none"),
  p.progressive.woco2[[1]] + theme(legend.position = "none"),
  p.progressive.wco2[[1]] + theme(legend.position = "none"),
  p.conservative.woco2[[1]] + theme(legend.position = "none"),
  p.conservative.wco2[[1]] + theme(legend.position = "none"),
  ncol = 2,
  labels = "auto",
  label_size = font.size
)

legend <- get_legend(p.default.woco2[[1]])

# Add labels at the very top
p.woco2 <- ggdraw() +
  draw_label(aes(label = "Without carbon price"),
             size = 6,
             fontface = "bold")

p.wco2 <- ggdraw() +
  draw_label(aes(label = "With ambitious carbon price"),
             size = 6,
             fontface = "bold")
  
p.colnames <- plot_grid(
  p.woco2,
  p.wco2,
  NULL,
  ncol = 2
)

p.panels <- plot_grid(
  p.colnames,
  p.panels,
  ncol = 1,
  rel_heights = c(0.03, 1)
)

# Total plot
p <- plot_grid(
  p.panels,
  NULL,
  legend,
  ncol = 3,
  rel_widths = c(1, 0.03, 0.22)
  )

ggsave("figures/SUPP_FIGURE16_subsidies_1.5C.png",
       width = 180, height = 180, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE16_subsidies_1.5C.svg",
       width = 180, height = 180, units = "mm", bg = "white")

```

## SUPPLEMENTARY FIGURE 17: Capacity backed by announced subsidies and demand-side regulation

```{r supplementary figure 17}

# Subsidy tracker from BNEF
data.subsidies.world <- data.subsidies.tracker %>% 
  filter(publication == 2023,
         region == "World") %>% 
  pull(value)

# Get capacity either operational of FID/Construction
data.cumcap.certain <- data.iea.enduse %>%
  filter(year %in% seq(2020, 2030),
         status %in% c("Decommissioned", "Operational", "FID/Construction")) %>%
  group_by(year) %>%
  summarise(cumcap.certain = sum(cumcap.sum))

# Get demand-backed capacity
data.cap.supported <- data.subsidies.annual %>%
  # First get end-use specific unique capacity along all year.cap
  group_by(scenario, carbonprice, enduse, year.cap) %>% 
  summarise(cap = unique(cap),
            cap.backed.demand = unique(cap.backed)) %>% 
  # Now summarise across all end-uses
  group_by(scenario, carbonprice, year.cap) %>% 
  summarise(cap = 1E-3*sum(cap),
            cap.backed.demand = 1E-3*sum(cap.backed.demand)) 

# Get total required subsidies for projects constructed in year.cap
data.subsidies.tot <- data.subsidies.annual %>% 
  group_by(scenario, carbonprice, year.cap) %>% 
  summarise(subsidies.total.required = sum(subsidies.annual))
  
# Continue
data.cap.supported <- data.cap.supported %>% 
  full_join(data.subsidies.tot) %>% 
  filter(year.cap %in% seq(2024, 2030)) %>% 
  rowwise() %>% 
  # Calculate capacity not backed by demand regulation
  mutate(cap.notbacked = max(0, cap - cap.backed.demand),
         cap.backed.demand = min(cap, cap.backed.demand)) %>% 
  group_by(scenario, carbonprice) %>% 
  arrange(year.cap) %>% 
  # Cumulative sums
   mutate(cumcap = cumsum(cap)) %>%
  # Join with capacity assumed as certain (operational/FID/construction)
  left_join(data.cumcap.certain, by = c("year.cap" = "year")) %>% 
  mutate(cumcap = cumcap + cumcap.certain) %>% 
  # Include total subsidies
  mutate(subsidies.2030.announced = data.subsidies.world) %>% 
  # Calculate remaining subsidies
  group_by(scenario, carbonprice) %>% 
  mutate(subsidies.left = subsidies.2030.announced - cumsum(subsidies.total.required)) %>% 
  # Calculate if projects can be supported from announced subsidies
  mutate(cap.backed.subsidies = case_when(
    # If subsidies are left, all new project announcement capacity can be built
    subsidies.left > 0 ~ cap.notbacked,
    # If already for project announcements in previous year no subsidies were
    # left, nothing can be built
    lag(subsidies.left) < 0 ~ 0,
    # If some subsidies are left (but not enough), only a part of new project
    # announcement capacity can be built
    subsidies.left < 0 ~ lag(subsidies.left)/subsidies.total.required * cap.notbacked)) %>% 
  # Account for cap.notbacked
  mutate(cap.notbacked = cap.notbacked - cap.backed.subsidies) %>%
  group_by(scenario, carbonprice) %>% 
  # Cumulative sums
  arrange(year.cap) %>% 
  mutate(cumcap.backed.demand = cumsum(cap.backed.demand),
         cumcap.backed.subsidies = cumsum(cap.backed.subsidies),
         cumcap.notbacked = cumsum(cap.notbacked))

# Plotting
data.plot <- data.cap.supported %>% 
  pivot_longer(cols = c("cumcap.certain", "cumcap.backed.demand", "cumcap.backed.subsidies", "cumcap.notbacked"),
               names_to = "variable",
               values_to = "value",
               names_prefix = "cumcap.") %>% 
  order.levels(scenario = c("Default", "Progressive", "Conservative"),
               variable = c("notbacked", "backed.subsidies", "backed.demand", "certain"))

colors.cap.type <- c(
  "notbacked" = "grey80",
  "backed.subsidies" = "#00A087FF",
  "backed.demand" = "#d8593b",
  "certain" = "#3C5488FF"
)

labels.cap.type <- c(
  "notbacked" = "Unsupported",
  "backed.subsidies" = "Supported by announced subsidies",
  "backed.demand" = "Supported by implemented demand-side policies",
  "certain" = "Operational/FID/Construction"
)

p <- ggplot() +
  geom_col(data = data.plot,
           mapping = aes(x = year.cap, y = value, fill = variable)) +
  scale_fill_manual(
    name = "Project announcements",
    values = colors.cap.type,
    labels = labels.cap.type) +
  facet_grid(scenario ~ carbonprice,
             labeller = labeller(carbonprice = c("TRUE" = "With ambitious carbon price",
                                                 "FALSE" = "Without carbon price"),
                                 scenario = c("Default" = "Central"))) +
  scale_x_continuous(name = "Year",
                     breaks = seq(2024, 2030, 2)) +
  ylab("Capacity (GW)") + 
  ggtitle("Global green hydrogen project announcements by policy support") +
  theme(strip.background =element_rect(fill = NA),
        axis.line = element_line(linewidth = 0),
        panel.border = element_rect(color = "black"))

# Save capacities backed by demand-side regulation or subsidies into csv
# This is read in in the beginning for the illustration in Figure 1
data.save <- data.plot %>%
  filter(variable != "notbacked",
         carbonprice == FALSE) %>% 
  group_by(scenario, carbonprice, year.cap) %>% 
  summarise(value = sum(value)) %>% 
  rename(year = year.cap)

write.csv(data.save, "data/fig1_capacity_backed_by_subsidies.csv", row.names = FALSE)

ggsave("figures/SUPP_FIGURE17_capacity_supported_by_policies.png",
       width = 120, height = 80, units = "mm", bg = "white", dpi = 600)
ggsave("figures/SUPP_FIGURE17_capacity_supported_by_policies.svg",
       width = 120, height = 80, units = "mm", bg = "white")

```

# SUPPLEMENTARY TABLES

## SUPPLEMENTARY TABLE 1: Number of validated projects

```{r supplementary table 1}

# IEA DB v2021
data.table.v2021 <- full_join(
  data.v2021.proj %>% 
    filter(status != "Decommissioned",
           reference != 0),
  data.v2021.orig.proj %>% 
    filter(status != "Decommissioned",
           reference != 0),
  by = c("reference" = "reference"))

# IEA DB v2021, projects announced for 2022
data.table.v2021 %>% 
  filter(year.y == 2022,
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

# IEA DB v2021, projects announced for 2023
data.table.v2021 %>% 
  filter(year.y == 2023,
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

# IEA DB v2021, projects announced for 2024-2030
data.table.v2021 %>% 
  filter(year.y %in% seq(2024, 2030),
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

# IEA DB v2022
data.table.v2022 <- full_join(
  data.v2022.proj %>% 
    filter(status != "Decommissioned",
           reference != 0),
  data.v2022.orig.proj %>% 
    filter(status != "Decommissioned",
           reference != 0),
  by = c("reference" = "reference"))

# IEA DB v2022, projects announced for 2022
data.table.v2022 %>% 
  filter(year.y == 2022,
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

# IEA DB v2022, projects announced for 2023
data.table.v2022 %>% 
  filter(year.y == 2023,
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

# IEA DB v2022, projects announced for 2024-2030
data.table.v2022 %>% 
  filter(year.y %in% seq(2024, 2030),
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

# IEA DB v2023
data.table.v2023 <- full_join(
  data.v2023.proj %>% 
    filter(status != "Decommissioned",
           reference != 0),
  data.v2023.orig.proj %>% 
    filter(status != "Decommissioned",
           reference != 0),
  by = c("reference" = "reference"))

# IEA DB v2023, projects announced for 2022
data.table.v2023 %>% 
  filter(year.y == 2022,
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

# IEA DB v2023, projects announced for 2023
data.table.v2023 %>% 
  filter(year.y == 2023,
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

# IEA DB v2023, projects announced for 2024-2030
data.table.v2023 %>% 
  filter(year.y %in% seq(2024, 2030),
         !is.na(capacity.y)) %>% 
  mutate(checked = replace_na(checked, TRUE)) %>% 
  summarise(
    projects.count = n(),
    projects.checked = sum(checked),
    projects.checked.capshare = sum(capacity.y * checked) / sum(capacity.y)) %>% 
  print()

```