Dutch Gas Transition Analysis (1990–2023)

Yaroslava Soroka

Date: 27.05.2026

Context and Relevance to Sustainable Development

The development of the Dutch gas market, which began with the discovery of the giant Groningen field in 1959, made the country a leader in gas production and exports in the EU, making blue fuel the basis of national energy policy and electricity generation (Corbeau, 2018). But after decades of dominance as Europe’s main energy hub, the Dutch gas market has reached a turning point: the cessation of production at the Groningen field has created fundamental challenges, ultimately leading to a large-scale transformation of the entire gas sector. (Riemersma et al., 2020). Additionally, the Netherlands is demonstrating an ambitious climate strategy aimed at achieving neutrality by 2050. The country is already outperforming the EU average, having reduced net emissions by 32.3% and setting an even more ambitious domestic target of a 60% reduction by 2030 (Erbach & Dewulf, 2024).

Therefore, the focus of this analysis lies on the recent transformation of the Dutch energy system, specifically the period marked by a rising share of renewable energy and a simultaneous decline in CO₂ emissions per capita. In this context, the central research question is: How did the share of renewable energy (X) and CO₂ emissions per capita (Y) change in the Netherlands between 2018 and 2023, and what factors explain this transition?

Figure: Renewable energy share vs. CO₂​ emissions per capita (2018–2023)

The data reveals a clear upward trend in the share of renewable energy alongside a simultaneous decline in CO₂ emissions per capita. This inverse correlation suggests that the Dutch energy system is undergoing a fundamental structural transformation toward decarbonization

1.1. Chapter …

Description

library(tidyverse)
library(eurostat)
library(dplyr)

# 2. Завантажуємо дані
gas_raw <- get_eurostat("nrg_cb_gas")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/df45c61508f0e62885a3089aceaccc5f.rds
Table  nrg_cb_gas  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/df45c61508f0e62885a3089aceaccc5f.rds

# 3. Обробка через пайп (%>%)
prod <- gas_raw %>%
  filter(geo == "NL") %>%                    # Тільки Нідерланди
  filter(unit == "TJ_GCV") %>%               # Тільки певна одиниця виміру
  filter(nrg_bal %in% c("IPRD")) %>%         # Тільки Виробництво 
  filter(siec == "G3000")  

exp<- gas_raw %>%
  filter(geo == "NL") %>%                    # Тільки Нідерланди
  filter(unit == "TJ_GCV") %>%               # Тільки певна одиниця виміру
  filter(nrg_bal %in% c( "EXP")) %>%        # Тільки Експорт
  filter(siec == "G3000")   


plot(prod$TIME_PERIOD, prod$values,    
     type = "l", col = "blue",lwd = 2,
     xaxt = "n", 
     xlab = "Year", ylab = "Gas Volume", 
     main = "Gas Production vs Exports in the Netherlands")
caption = "Source: Eurostat [nrg_cb_gas]"


lines(exp$TIME_PERIOD, exp$values, col = "red", lwd = 2)

axis.Date(1, at = seq(min(prod$TIME_PERIOD), max(prod$TIME_PERIOD), by = "5 years"), format = "%Y")

legend("topright", 
       legend = c("Production", "Exports"), 
       col = c("blue", "red"), 
       lty = 1, lwd = 2, cex = 0.7)

grid()

2.1. Name

Description

library(eurostat)
library(dplyr)


ren <- get_eurostat("nrg_ind_ren")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/78ad026d0e415a56d0f1b89d1e337078.rds
Table  nrg_ind_ren  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/78ad026d0e415a56d0f1b89d1e337078.rds

co2 <- get_eurostat("env_air_gge")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/58a9992b00c628ab521a3b92df80b45b.rds
Table  env_air_gge  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/58a9992b00c628ab521a3b92df80b45b.rds

co2_smallest_order <- co2 %>%
  pull(values) %>%
  order()
co2_smallest <- co2_smallest_order[1]

ren_nl = ren %>%
  filter(geo == "NL", TIME_PERIOD >= "2018-01-01",TIME_PERIOD <= "2023-12-31",nrg_bal=="REN" )

co2_nl = co2 %>%
  filter(geo == "NL", airpol == "CO2", TIME_PERIOD >= "2018-01-01", TIME_PERIOD <= "2023-12-31", src_crf == "TOTXMEMO", unit=="MIO_T")

ren_nl <- ren_nl %>%
  rename(renewable = values)
co2_nl <- co2_nl %>%
  rename(co2 = values)

data_final <- ren_nl %>%
  inner_join(co2_nl, by = "TIME_PERIOD")

cor(data_final$renewable, data_final$co2)

[1] -0.9897367

model <- lm(co2 ~ renewable, data = data_final)
summary(model)

Call:
lm(formula = co2 ~ renewable, data = data_final)

Residuals:
     1      2      3      4      5      6 
-1.026 -1.036  2.994  2.565 -1.192 -2.305 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 192.2778     3.7595   51.14 8.75e-07 ***
renewable    -3.9464     0.2849  -13.85 0.000157 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.47 on 4 degrees of freedom
Multiple R-squared:  0.9796,    Adjusted R-squared:  0.9745 
F-statistic: 191.9 on 1 and 4 DF,  p-value: 0.0001575

plot(data_final$renewable, data_final$co2, col = "blue",lwd = 2,
     xlab = "Renewable energy (%)",
     ylab = "CO2 per capita (MIO_T)",
     main = "Relationship between Renewable Energy and CO2")
model <- lm(co2 ~ renewable, data = data_final)
abline(model, col = "red",lwd = 2)

legend("topright", 
       legend = c("Experimental data", "Model data"), 
       col = c("blue", "red"), 
       lty = 1, lwd = 2, cex = 0.7)

grid()

2.2. Name

Description

library(tidyverse)
library(eurostat)

# 1. Завантажуємо дані (якщо ще не завантажені)
raw_energy <- get_eurostat("nrg_bal_s")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/d3e73add764e8b0bdd8298c9ee0a57ec.rds
Table  nrg_bal_s  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/d3e73add764e8b0bdd8298c9ee0a57ec.rds

raw_gdp    <- get_eurostat("nama_10_gdp")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/5590e4702bb787ef42ec6e58b1b97c07.rds
Table  nama_10_gdp  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/5590e4702bb787ef42ec6e58b1b97c07.rds

energy_gdp_data <- raw_gdp %>%
  # Фільтруємо ВВП
  filter(geo == "NL", unit == "CLV15_MEUR", na_item == "B1GQ") %>%
  select(TIME_PERIOD, gdp_val = values) %>%
  # Приєднуємо дані по енергії
  inner_join(
    raw_energy %>%
      filter(geo == "NL", unit == "GWH", nrg_bal == "AFC") %>%
      group_by(TIME_PERIOD) %>%
      summarise(energy_val = sum(values, na.rm = TRUE)),
    by = "TIME_PERIOD"
  ) %>%
  # Рахуємо індекси (перший рік = 100)
  mutate(gdp_index = (gdp_val / first(gdp_val)) * 100,
         energy_index = (energy_val / first(energy_val)) * 100)
str(energy_gdp_data)

tibble [30 × 5] (S3: tbl_df/tbl/data.frame)
 $ TIME_PERIOD : Date[1:30], format: "1995-01-01" ...
 $ gdp_val     : num [1:30] 472563 488728 509492 533171 560061 ...
 $ energy_val  : num [1:30] 1315474 1393935 1345128 1348896 1348048 ...
 $ gdp_index   : num [1:30] 100 103 108 113 119 ...
 $ energy_index: num [1:30] 100 106 102 103 102 ...

y_min = min(energy_gdp_data$energy_index)
y_max = max(energy_gdp_data$gdp_index)

plot(energy_gdp_data$TIME_PERIOD, energy_gdp_data$gdp_index,
     ylim = c(y_min,y_max),
     type = "l", col = "blue",lwd = 2,
     xaxt = "n", 
     xlab = "Year", ylab = "Index Value, 1995 = 100%", 
     main = "Economic Growth vs Energy Use (Netherlands)",
     )

mtext("Source: Eurostat", side = 1, line = 4, adj = 1)

lines(energy_gdp_data$TIME_PERIOD, energy_gdp_data$energy_index, col = "red", lwd = 2)

axis.Date(1, at = seq(min(energy_gdp_data$TIME_PERIOD), max(energy_gdp_data$TIME_PERIOD), by = "4 years"), format = "%Y")

legend("topleft", 
       legend = c("GDP Growth", "Energy Use"), 
       col = c("blue", "red"), 
       lty = 1, lwd = 2, cex = 0.7)

grid()

3.1. Name

Description

library(eurostat)
library(dplyr)
library(tidyr)

# 2. ЗАВАНТАЖЕННЯ ДАНИХ 
# nrg_inf_epc — це код бази даних Eurostat для потужностей електростанцій
capacity_raw <- get_eurostat("nrg_inf_epc", time_format = "date")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/0bbff8a7e565f26af137d75b3ab7b781.rds
Table  nrg_inf_epc  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/0bbff8a7e565f26af137d75b3ab7b781.rds

# 3. ОБРОБКА ДАНИХ
capacity_full <- capacity_raw %>%
  filter(geo == "NL", 
         siec %in% c("RA420", "RA300", "N9000", "RA100")) %>% 
  mutate(Year = as.numeric(format(TIME_PERIOD, "%Y")),
         tech = case_when(
           siec == "RA420" ~ "Solar",
           siec == "RA300" ~ "Wind",
           siec == "N9000" ~ "Nuclear",
           siec == "RA100" ~ "Hydro",
           TRUE ~ "Other"
         )) %>%
  group_by(Year, tech) %>%
  summarise(values = sum(values, na.rm = TRUE), .groups = 'drop') %>%
  pivot_wider(names_from = tech, values_from = values)

# 4. ПЕРЕВІРКА РЕЗУЛЬТАТУ

print("Оброблена таблиця capacity_full:")

[1] "Оброблена таблиця capacity_full:"

print(head(capacity_full))

#gas_prices_raw <- get_eurostat("nrg_pc_202", time_format = "date")

gas_prices_clean <- gas_prices_raw %>%
  filter(geo == "NL", 
         unit == "KWH", 
         tax == "I_TAX", 
         currency == "EUR") %>%
  mutate(Year = as.numeric(format(TIME_PERIOD, "%Y"))) %>%
  group_by(Year) %>%
  summarise(gas_price = mean(values, na.rm = TRUE), .groups = 'drop')

# Об'єднуємо за роком
geo_data_final <- gas_prices_clean %>%
  inner_join(capacity_full, by = "Year") %>%
  filter(Year >= 2018)

# Дивимось результат
print(geo_data_final)

# 1. Перераховуємо все у відсотки (Індексація: 2018 рік = 100)
geo_data_indexed <- geo_data_final %>%
  mutate(
    Solar_idx = (Solar / Solar[Year == 2018]) * 100,
    Wind_idx = (Wind / Wind[Year == 2018]) * 100,
    Nuclear_idx = (Nuclear / Nuclear[Year == 2018]) * 100,
    Hydro_idx = (Hydro / Hydro[Year == 2018]) * 100,
    Gas_idx = (gas_price / gas_price[Year == 2018]) * 100
  )

y_min = min(geo_data_indexed$Solar_idx)
y_max = max(geo_data_indexed$Solar_idx)

plot(geo_data_indexed$Year, geo_data_indexed$Solar_idx,
     ylim = c(y_min,y_max),
     type = "l", col = "red",lwd = 2,
     xaxt = "n", 
     xlab = "Year", ylab = "Growth relative to 2018, 2018 = 100%", 
     main = "Netherlands Energy Landscape (2018-2024)",)

grid()

mtext("Source: Eurostat Data Analysis (2026)", side = 1, line = 4, adj = 1)

lines(geo_data_indexed$Year, geo_data_indexed$Wind_idx, col = "blue", lwd = 2)

lines(geo_data_indexed$Year, geo_data_indexed$Gas_idx, col = "lightblue", lwd = 2)
lines(geo_data_indexed$Year, geo_data_indexed$Nuclear_idx, col = "black", lwd = 2)

lines(geo_data_indexed$Year, geo_data_indexed$Hydro_idx, col = "darkblue", lwd = 2)


axis(1, at = seq(min(geo_data_indexed$Year), max(geo_data_indexed$Year)))

legend("topleft", 
       legend = c("Solar", "Wind","Gas Price","Nuclear","Hydro"), 
       col = c("red", "blue","lightblue","black","darkblue"), 
       lty = 1, lwd = 2, cex = 0.7)

3.2. Name

Description

library(eurostat)
library(dplyr)
library(ggplot2)

# 1. Отримання даних з правильними ID
tax_raw <- get_eurostat("env_ac_tax")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/7ac6eba40fc48406f3b5a838b1bc7300.rds
Table  env_ac_tax  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/7ac6eba40fc48406f3b5a838b1bc7300.rds

# env_ac_epneis1 — видатки (замість env_ac_exp1, який видавав помилку)
exp_raw <- get_eurostat("env_ac_epneis1")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/9c36de815110a194639e1ba9b01144c6.rds
Table  env_ac_epneis1  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/9c36de815110a194639e1ba9b01144c6.rds

# 1. Processing Environmental Tax Revenues
tax_clean <- tax_raw %>%
  filter(geo == "NL", unit == "MIO_EUR") %>%
  mutate(Year = as.numeric(format(TIME_PERIOD, "%Y"))) %>%
  group_by(Year) %>%
  summarise(tax_revenue = sum(values, na.rm = TRUE), .groups = 'drop')

# 2. Processing Environmental Protection Expenditure
# Here we rename the variable to the full term you requested
exp_clean <- exp_raw %>%
  filter(geo == "NL", unit == "MIO_EUR") %>%
  mutate(Year = as.numeric(format(TIME_PERIOD, "%Y"))) %>%
  group_by(Year) %>%
  summarise(environmental_protection_expenditure = sum(values, na.rm = TRUE), .groups = 'drop')

tax_clean <-tax_clean %>%
  filter(Year >= 2018,Year <= 2022)

y_min = min(exp_clean$environmental_protection_expenditure)
y_max = max(tax_clean$tax_revenue)+10000


plot(tax_clean$Year, tax_clean$tax_revenue,
     ylim = c(y_min,y_max),
     xlim = c(2017.7,2022.3),
     type = "l", col = "red",lwd = 2,
     xaxt = "n", 
     xlab = "Year", ylab = "Amount (Million EUR)", 
     main = "Netherlands: Environmental Fiscal Trends",)

text(tax_clean$Year, tax_clean$tax_revenue, labels = tax_clean$tax_revenue, pos=1)

grid()

axis(1, at = seq(min(tax_clean), max(tax_clean$Year)))

legend("topleft", 
       legend = c("Tax Revenue", "Environmental Protection Expenditure"), 
       col = c("red", "blue"), 
       lty = 1, lwd = 2, cex = 0.7)


mtext("Source: Eurostat Data Analysis (2026)", side = 1, line = 4, adj = 1)

lines(exp_clean$Year, exp_clean$environmental_protection_expenditure, col = "blue", lwd = 2)
text(exp_clean$Year, exp_clean$environmental_protection_expenditure, labels = round(exp_clean$environmental_protection_expenditure), pos=3)

3.3. Name

Description

library(eurostat)
library(dplyr)

rd_e_gerdtot_raw <- get_eurostat("rd_e_gerdtot")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/88347b7c6c418358a0d682033c7044a1.rds
Table  rd_e_gerdtot  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/88347b7c6c418358a0d682033c7044a1.rds

sdg_09_40_raw <- get_eurostat("sdg_09_40")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/15abb4dea79736a34ef5ff3e0c5369e5.rds
Table  sdg_09_40  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/15abb4dea79736a34ef5ff3e0c5369e5.rds

# 1. Get R&D Data (X-axis equivalent: Time)
rd_recent <- rd_e_gerdtot_raw %>%
  mutate(Year = as.numeric(format(TIME_PERIOD, "%Y"))) %>%
  filter(geo == "NL", sectperf == "TOTAL", unit == "EUR_HAB", Year >= 2018,Year<=2024) %>%
  select(Year, values)

# 2. Get Patent Data
pat_recent <- sdg_09_40_raw %>%
  mutate(Year = as.numeric(format(TIME_PERIOD, "%Y"))) %>%
  filter(geo == "NL", unit == "NR", Year >= 2018,Year<=2024 ,coop_ptn == "APPL" ) %>%
  select(Year,  values)


y_min = min(rd_recent$values)
y_max = max(pat_recent$values)+2000

plot(pat_recent$Year, pat_recent$values,
     ylim = c(y_min,y_max),
     xlim = c(2017.8,2024.2),
     type = "l", col = "red",lwd = 2,
     xaxt = "n", 
     xlab = "Year", ylab = "Value (Units vary by indicator)", 
     main = "Netherlands: Innovation Indicators (2018-2023)\n Yearly trends in R&D investment\n and Climate Technology patents",)

text(pat_recent$Year, pat_recent$values, labels = pat_recent$values, pos=1)

grid()

axis(1, at = seq(2018, 2024))

legend("topleft", 
       legend = c("Climate Tech Patents (Count)", "R&D Expenditure (EUR/hab)"), 
       col = c("red", "blue"), 
       lty = 1, lwd = 2, cex = 0.7)


mtext("Source: Eurostat Data Analysis (2026)", side = 1, line = 4, adj = 1)

lines(rd_recent$Year, rd_recent$values, col = "blue", lwd = 2)
text(rd_recent$Year, rd_recent$values, labels = round(rd_recent$values), pos=3)

4.1. Name

Description

library(eurostat)
library(dplyr)
library(ggplot2)

# 1. Завантажуємо дані
gge_raw <- get_eurostat("env_air_gge")

Dataset query already saved in cache_list.json...
Reading cache file C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/58a9992b00c628ab521a3b92df80b45b.rds
Table  env_air_gge  read from cache file:  C:\Users\Admin\AppData\Local\Temp\RtmpMrh3Mh/eurostat/58a9992b00c628ab521a3b92df80b45b.rds

# 2. Фільтруємо дані 
gge_2023 <- gge_raw %>%
  filter(geo == "NL") %>%
  mutate(Year = as.numeric(format(TIME_PERIOD, "%Y"))) %>%
  # Беремо період 2018-2023
  filter(Year >= 2018 & Year <= 2023) %>%
  # Фільтруємо парникові гази та основні сектори
  filter(airpol == "GHG") %>%
  filter(src_crf %in% c("CRF1A1", "CRF1A2", "CRF1A3", "CRF1A4", "CRF3")) %>%
  # Групуємо, щоб прибрати дублікати
  group_by(Year, src_crf) %>%
  summarise(values = mean(values, na.rm = TRUE), .groups = 'drop')

# 3. Додаємо назви секторів
sector_labels <- c(
  "CRF1A1" = "Energy",
  "CRF1A2" = "Industry",
  "CRF1A3" = "Transport",
  "CRF1A4" = "Residential",
  "CRF3"   = "Agriculture"
)
gge_2023$Sector <- sector_labels[gge_2023$src_crf]

gge_2023_Energy <-gge_2023 %>% filter(Sector=="Energy")
gge_2023_Industry <-gge_2023 %>% filter(Sector=="Industry")
gge_2023_Transport <-gge_2023 %>% filter(Sector=="Transport")
gge_2023_Residential <-gge_2023 %>% filter(Sector=="Residential")
gge_2023_Agriculture <-gge_2023 %>% filter(Sector=="Agriculture")

# 4. МАЛЮЄМО

y_min = min(gge_2023$values)-1000
y_max = max(gge_2023$values)+3000

plot(gge_2023_Energy$Year, gge_2023_Energy$values,
     ylim = c(y_min,y_max),
     xlim = c(2017.8,2023.2),
     type = "l", col = "red",lwd = 2,
     xaxt = "n", 
     xlab = "Year", ylab = "Thousand Tonnes", 
     main = "Netherlands GHG Emissions (2018-2023)\n Progress towards climate goals by sector",)

text(gge_2023_Energy$Year, gge_2023_Energy$values, labels = round(gge_2023_Energy$values), pos=1, col = "red",cex=0.7)

grid()

axis(1, at = seq(2018, 2023))

legend("topright", 
       legend = c("Energy", "Industry","Transport","Residential","Agriculture"), 
       col = c("red", "blue","violet","orange","green"), 
       lty = 1, lwd = 2, cex = 0.7)


mtext("Source: Eurostat Data Analysis (2026)", side = 1, line = 4, adj = 1)

lines(gge_2023_Industry$Year, gge_2023_Industry$values, col = "blue", lwd = 2)
text(gge_2023_Industry$Year, gge_2023_Industry$values, labels = round(gge_2023_Industry$values), pos=3,cex = 0.7, col = "blue")

lines(gge_2023_Transport$Year, gge_2023_Transport$values, col = "violet", lwd = 2)
text(gge_2023_Transport$Year, gge_2023_Transport$values, labels = round(gge_2023_Transport$values), pos=1,cex = 0.7, col = "violet")

lines(gge_2023_Residential$Year, gge_2023_Residential$values, col = "orange", lwd = 2)
text(gge_2023_Residential$Year, gge_2023_Residential$values, labels = round(gge_2023_Residential$values), pos=3 ,cex = 0.7, col = "orange")

lines(gge_2023_Agriculture$Year, gge_2023_Agriculture$values, col = "green", lwd = 2)
text(gge_2023_Agriculture$Year, gge_2023_Agriculture$values, labels = round(gge_2023_Agriculture$values), pos=1,cex = 0.7,col = "green")

NA
NA

5.1. Name

Description

library(eurostat)
library(dplyr)

# ЗАВАНТАЖЕННЯ 
#ge_raw <- get_eurostat("env_air_gge")

ge_clean <- ge_raw %>%
  filter(geo == "NL", airpol == "GHG", unit == "THS_T") %>%
  mutate(Year = as.numeric(format(TIME_PERIOD, "%Y"))) %>%
  filter(Year %in% 2018:2023)

# Створення категорій та групування
compare_df <- ge_clean %>%
  filter(grepl("MEMO", src_crf)) %>%
  mutate(Indicator = case_when(
    grepl("X4|XL", src_crf) ~ "Excl",
    TRUE ~ "Incl"
  )) %>%
  group_by(Year, Indicator) %>%
  summarise(values = mean(values, na.rm = TRUE), .groups = 'drop')
# Рахуємо математику перед графіком
analysis_data <- compare_df %>% 
  pivot_wider(names_from = Indicator, values_from = values)

correlation <- cor(analysis_data$Excl, analysis_data$Incl)
reg_summary <- summary(lm(Incl ~ Excl, data = analysis_data))


# Розділення даних
excl <- subset(compare_df, Indicator == "Excl")
incl <- subset(compare_df, Indicator == "Incl")

# Відступи
x_pad <- 0.3
y_pad <- diff(range(compare_df$values)) * 0.05

# Поля для legend/caption
par(mar = c(7,4,5,2))


y_min = min(excl$values)-1000
y_max = max(excl$values)+3000

# Основний графік
plot(excl$Year, excl$values,
     xlim = c(2018,2023),
     ylim = c(y_min,y_max),
     type = "o", pch = 16, lwd = 2,
     col = "red",
     main = "Comparison of Emissions: Total vs. Total with LULUCF",
     sub = paste0("Correlation: ",round(correlation, 4)," | R²: ",round(reg_summary$r.squared, 4)     ),
     xlab = "Year",  ylab = "Thousand Tonnes (CO2 eq.)",
     col.sub = "grey30")


grid()

# Власна вісь X
axis(1, at = seq(2018, 2023))

# Друга лінія
lines(incl$Year, incl$values, type = "o", pch = 16,lwd = 2, col = "blue")

# Лінія тренду для Excl
abline(lm(values ~ Year, data = excl),col = "red",lty = 2,  lwd = 1)

# Лінія тренду для Incl
abline(lm(values ~ Year, data = incl),col = "blue",  lty = 2,  lwd = 1)

# Легенда
legend("topright", legend = c("Total (Excl. LULUCF)","Total (Incl. LULUCF)"),
       col = c("red", "blue"),
       lty = 1, lwd = 2, pch = 16, cex = 0.7)

# Caption
mtext(  "Source: Eurostat [env_air_gge] | Dashed lines represent linear trends",
  side = 1,  line = 5,  cex = 0.8,  col = "grey40"
)

NA
NA