#--------------------------------------------
########## Figure 3 ##########
#--------------------------------------------

rm(list = ls())
setwd("./Replication_AQ_EV/Map/")

library(showtext) 
showtext_auto(enable = TRUE) 
font_add("cnfont", regular = "Tiejili Regular.ttf")
cnfont <- "cnfont"
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
library(RColorBrewer)
library(sf)
library(tidyverse)
library(ggspatial)
library(showtext) 
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
setwd("./Replication_AQ_EV/")

city_map <- read_sf("Map/Location/ChinaCity.geojson") %>% 
  st_transform(mycrs) %>% 
  st_make_valid()
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
read_sf("Map/Haianxian/Haianxian.shp") %>% 
  st_transform(mycrs) -> hax
read_sf("Map/Jiuduanxian/Jiuduanxian.shp") %>% 
  st_transform(mycrs) -> jdx 
read_sf("Map/Qinhuaixian.geojson") %>% 
  st_transform(mycrs) -> qh 
read_sf("Map/Location/County.shp") %>% 
  st_transform(mycrs) %>% 
  rename(countycode = PAC, county = NAME) %>% 
  st_make_valid() -> county

dt_data <- read.csv("Data/EVCAperStation.csv")

dt_data$ChargePs_mean[dt_data$ChargePs_mean == 0] <- NA
dt_data$EVCA_PM25 <- (dt_data$ChargePs_mean / 141.9836) * 0.00147664 * 1000
dt_data$EVCA_NO2 <- (dt_data$ChargePs_mean / 141.9836) * 0.00208288 * 1000
dt_data$EVCA_PM25_CIup <- (dt_data$ChargePs_mean / 141.9836) * 0.002141128 * 1000
dt_data$EVCA_PM25_CIdow <- (dt_data$ChargePs_mean / 141.9836) * 0.000812152 * 1000
dt_data$EVCA_NO2_CIup <- (dt_data$ChargePs_mean / 141.9836) * 0.002421348 * 1000
dt_data$EVCA_NO2_CIdow <- (dt_data$ChargePs_mean / 141.9836) * 0.001484052 * 1000

county %>% 
  left_join(dt_data, by = c("市" = "city")) -> dfcounty2

###### Map: PM2.5 ###### 
dfcounty2$PM25_binned <- cut(dfcounty2$EVCA_PM25, 
                             breaks = c(-Inf, 0.1, 0.5, 2, 4, Inf),
                             labels = c("≤0.0, >-0.1", "≤-0.1, >-0.5", "≤-0.5, >-2.0", 
                                        "≤-2.0, >-4.0", "≤-4.0"),
                             include.lowest = TRUE)
dfcounty2$PM25_binned <- factor(dfcounty2$PM25_binned, 
                                levels = c("≤-4.0", "≤-2.0, >-4.0", "≤-0.5, >-2.0",
                                           "≤-0.1, >-0.5", "≤0.0, >-0.1"))
pm25_map <- ggplot() + 
  geom_sf(data = city_map, color = NA) +
  geom_sf(data = hax, color = "#0055AA", linewidth = 0.3) +
  geom_sf(data = jdx, color = "black", linewidth = 0.6) + 
  geom_sf(data = dfcounty2, aes(fill = PM25_binned), color = "grey70", shape = 21) +
  geom_sf(data = qh, color = "#EFC000FF", linewidth = 2) + 
  scale_fill_manual(
    values = c("#356B5C", "aquamarine4", "aquamarine3", "aquamarine", "#B3FFE0"),
    name = expression("PM"["2.5"]*" (μg/m³)"),
    na.value = "grey97" 
  ) +
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = c(0.06, 0.06), 
    legend.justification = c("left", "bottom"), 
    legend.key.size = unit(7, "mm"),
    legend.text = element_text(size = 22),
    legend.title = element_text(size = 22), 
    axis.title = element_blank(), 
    axis.text = element_blank(), 
    panel.grid = element_blank()
  ) +
  annotation_scale(location = "bl", 
                   width_hint = 0.3, 
                   text_family = cnfont) + 
  labs(title = "",
       caption = "",
       fill = "") + 
  annotation_north_arrow(
    location = "tl", 
    which_north = "false",
    style = north_arrow_fancy_orienteering(text_family = cnfont),
    pad_x = unit(0.01, "in"),
    pad_y = unit(0.75, "in") 
  )


###### Map: NO2 ###### 
dfcounty2$NO2_binned <- cut(dfcounty2$EVCA_NO2, 
                             breaks = c(-Inf, 0.5, 1, 3, 6, Inf),
                             labels = c("≤0, >-0.5", "≤-0.5, >-1.0", "≤-1.0, >-3.0", 
                                        "≤-3.0, >-6.0", "≤-6.0"),
                             include.lowest = TRUE)
dfcounty2$NO2_binned <- factor(dfcounty2$NO2_binned, 
                                levels = c("≤-6.0", "≤-3.0, >-6.0", "≤-1.0, >-3.0",
                                           "≤-0.5, >-1.0","≤0, >-0.5"))
no2_map <- ggplot() + 
  geom_sf(data = city_map, color = NA) +
  geom_sf(data = hax, color = "#0055AA", linewidth = 0.3) +
  geom_sf(data = jdx, color = "black", linewidth = 0.6) + 
  geom_sf(data = dfcounty2, aes(fill = NO2_binned), color = "grey70", shape = 21) +
  geom_sf(data = qh, color = "#EFC000FF", linewidth = 2) + 
  scale_fill_manual(
    values = c("indianred4", "indianred", "indianred2", "#FF9999", "#FFCCCC"),
    name = expression("NO"["2"]*" (μg/m³)"),
    na.value = "grey97"
  ) +
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = c(0.06, 0.06), 
    legend.justification = c("left", "bottom"),
    legend.key.size = unit(7, "mm"), 
    legend.text = element_text(size = 22), 
    legend.title = element_text(size = 22),
    axis.title = element_blank(),
    axis.text = element_blank(), 
    panel.grid = element_blank()
  ) +
  annotation_scale(location = "bl", 
                   width_hint = 0.3, 
                   text_family = cnfont) + 
  labs(title = "",
       caption = "",
       fill = "") + 
  annotation_north_arrow(
    location = "tl", 
    which_north = "false",
    style = north_arrow_fancy_orienteering(text_family = cnfont),
    pad_x = unit(0.01, "in"),
    pad_y = unit(0.75, "in")
  )



###### EVCA ###### 
df_mean <- dt_data %>%
  group_by(region) %>%
  summarize(
    EVCA_PM25_mean = mean(EVCA_PM25, na.rm = TRUE),
    EVCA_NO2_mean = mean(EVCA_NO2, na.rm = TRUE),
    EVCA_PM25_CIup_mean = mean(EVCA_PM25_CIup, na.rm = TRUE),
    EVCA_PM25_CIdow_mean = mean(EVCA_PM25_CIdow, na.rm = TRUE),
    EVCA_NO2_CIup_mean = mean(EVCA_NO2_CIup, na.rm = TRUE),
    EVCA_NO2_CIdow_mean = mean(EVCA_NO2_CIdow, na.rm = TRUE)
  )

df_long <- df_mean %>%
  pivot_longer(
    cols = c(EVCA_PM25_mean, EVCA_NO2_mean),
    names_to = "pollutant",
    values_to = "mean"
  ) %>%
  mutate(
    ymin = ifelse(pollutant == "EVCA_PM25_mean", EVCA_PM25_CIdow_mean, EVCA_NO2_CIdow_mean),
    ymax = ifelse(pollutant == "EVCA_PM25_mean", EVCA_PM25_CIup_mean, EVCA_NO2_CIup_mean),
    type = ifelse(pollutant == "EVCA_PM25_mean", "PM2.5", "NO2")
  )

evca_plot <- ggplot(df_long, aes(x = type, y = mean, fill = region)) +
  geom_bar(stat = "identity", position = "dodge", alpha = 0.6) +
  geom_errorbar(aes(ymin = ymin, ymax = ymax), position = position_dodge(width = 0.9), width = 0.2) +
  labs(y = "Estimated reduction (μg/m³)", x = "Pollutant Type") +
  scale_fill_manual(values = c("North" = "#EDB176", "South" = "#6959cd"), name = "Region") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, size = 22),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 22),
        axis.text.y = element_text(size = 22),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 22),
        plot.title = element_blank(),
        axis.line.x = element_line(size = 0.8, color = "black"), 
        axis.line.y = element_line(size = 0.8, color = "black"),
        axis.ticks = element_line(size = 0.8, color = "black"))



dt_data <- read.csv("Data/PollutantsbyRegion.csv")

dt_data$EVCA_NO2 <- (dt_data$sum_NO2 * 0.00008) * (100000 / 141.9836)
dt_data$EVCA_PM25 <- (dt_data$sum_PM25 * 0.00004) * (100000 / 141.9836)
dt_data$EVCA_PM25_CIup <- (dt_data$sum_PM25 * 0.000058) * (100000 / 141.9836)
dt_data$EVCA_PM25_CIdow <- (dt_data$sum_PM25 * 0.000022) * (100000 / 141.9836)
dt_data$EVCA_NO2_CIup <- (dt_data$sum_NO2 * 0.000093) * (100000 / 141.9836)
dt_data$EVCA_NO2_CIdow <- (dt_data$sum_NO2 * 0.000059) * (100000 / 141.9836)

county %>% 
  left_join(dt_data, by = c("市" = "city")) -> dfcounty2

dfcounty2$PM25_binned <- cut(dfcounty2$EVCA_PM25, 
                             breaks = c(-Inf, 0.5, 1.0, 1.5, 1.8, Inf),
                             labels = c("≤0.0, >-0.5", "≤-0.5, >-1.0", "≤-1.0, >-1.5", 
                                        "≤-1.5, >-1.8", "≤-1.8"),
                             include.lowest = TRUE)
dfcounty2$PM25_binned <- factor(dfcounty2$PM25_binned, 
                                levels = c("≤-1.8", "≤-1.5, >-1.8", "≤-1.0, >-1.5",
                                           "≤-0.5, >-1.0", "≤0.0, >-0.5"))

pm25_map_ca <- ggplot() + 
  geom_sf(data = city_map, color = NA) +
  geom_sf(data = hax, color = "#0055AA", linewidth = 0.3) +
  geom_sf(data = jdx, color = "black", linewidth = 0.6) + 
  geom_sf(data = dfcounty2, aes(fill = PM25_binned), color = "grey70", shape = 21) +
  geom_sf(data = qh, color = "#EFC000FF", linewidth = 2) + 
  scale_fill_manual(
    values = c("#356B5C", "aquamarine4", "aquamarine3", "aquamarine", "#B3FFE0"),
    name = expression("PM"["2.5"]*" (μg/m³)"),
    na.value = "grey97" 
  ) +
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = c(0.06, 0.06), 
    legend.justification = c("left", "bottom"), 
    legend.key.size = unit(7, "mm"), 
    legend.text = element_text(size = 26), 
    legend.title = element_text(size = 26),
    axis.title = element_blank(), 
    axis.text = element_blank(),
    panel.grid = element_blank()
  ) +
  annotation_scale(location = "bl", 
                   width_hint = 0.3, 
                   text_family = cnfont) + 
  labs(title = "",
       caption = "",
       fill = "") + 
  annotation_north_arrow(
    location = "tl", 
    which_north = "false",
    style = north_arrow_fancy_orienteering(text_family = cnfont),
    pad_x = unit(0.01, "in"),
    pad_y = unit(0.75, "in") 
  )


dfcounty2$NO2_binned <- cut(dfcounty2$EVCA_NO2, 
                            breaks = c(-Inf, 0.5, 1.0, 1.5, 2.5, Inf),
                            labels = c("≤0.0, >-0.5", "≤-0.5, >-1.0", "≤-1.0, >-1.5", 
                                       "≤-1.5, >-2.5", "≤-2.5"),
                            include.lowest = TRUE)

dfcounty2$NO2_binned <- factor(dfcounty2$NO2_binned, 
                               levels = c("≤-2.5", "≤-1.5, >-2.5", "≤-1.0, >-1.5",
                                          "≤-0.5, >-1.0", "≤0.0, >-0.5"))

no2_map_ca <- ggplot() + 
  geom_sf(data = city_map, color = NA) +
  geom_sf(data = hax, color = "#0055AA", linewidth = 0.3) +
  geom_sf(data = jdx, color = "black", linewidth = 0.6) + 
  geom_sf(data = dfcounty2, aes(fill = NO2_binned), color = "grey70", shape = 21) +
  geom_sf(data = qh, color = "#EFC000FF", linewidth = 2) + 
  scale_fill_manual(
    values = c("indianred4", "indianred", "indianred2", "#FF9999", "#FFCCCC"),
    name = expression("NO"["2"]*" (μg/m³)"),
    na.value = "grey97"
  ) +
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = c(0.06, 0.06), 
    legend.justification = c("left", "bottom"),
    legend.key.size = unit(7, "mm"), 
    legend.text = element_text(size = 26),
    legend.title = element_text(size = 26), 
    axis.title = element_blank(),
    axis.text = element_blank(),
    panel.grid = element_blank()
  ) +
  annotation_scale(location = "bl", 
                   width_hint = 0.3, 
                   text_family = cnfont) + 
  labs(title = "",
       caption = "",
       fill = "") + 
  annotation_north_arrow(
    location = "tl", 
    which_north = "false",
    style = north_arrow_fancy_orienteering(text_family = cnfont),
    pad_x = unit(0.01, "in"), 
    pad_y = unit(0.75, "in")
  )

plot_E <- no2_map_ca + labs(tag = expression(bold("D. estimated NO"["2"]*" reduction due to 100 MWh"))) + theme(plot.tag = element_text(size = 26),
                                                                                                        plot.tag.position = c(0.5, 0.96))
plot_F <- pm25_map_ca + labs(tag = expression(bold("E. estimated PM"["2.5"]*" reduction due to 100 MWh"))) + theme(plot.tag = element_text(size = 26),
                                                                                                           plot.tag.position = c(0.5, 0.96))


###### EVCA ###### 
df_mean <- dt_data %>%
  group_by(region) %>%
  summarize(
    EVCA_PM25_mean = mean(EVCA_PM25, na.rm = TRUE),
    EVCA_NO2_mean = mean(EVCA_NO2, na.rm = TRUE),
    EVCA_PM25_CIup_mean = mean(EVCA_PM25_CIup, na.rm = TRUE),
    EVCA_PM25_CIdow_mean = mean(EVCA_PM25_CIdow, na.rm = TRUE),
    EVCA_NO2_CIup_mean = mean(EVCA_NO2_CIup, na.rm = TRUE),
    EVCA_NO2_CIdow_mean = mean(EVCA_NO2_CIdow, na.rm = TRUE)
  )

df_long <- df_mean %>%
  pivot_longer(
    cols = c(EVCA_PM25_mean, EVCA_NO2_mean),
    names_to = "pollutant",
    values_to = "mean"
  ) %>%
  mutate(
    ymin = ifelse(pollutant == "EVCA_PM25_mean", EVCA_PM25_CIdow_mean, EVCA_NO2_CIdow_mean),
    ymax = ifelse(pollutant == "EVCA_PM25_mean", EVCA_PM25_CIup_mean, EVCA_NO2_CIup_mean),
    type = ifelse(pollutant == "EVCA_PM25_mean", "PM2.5", "NO2")
  )

evelec_plot <- ggplot(df_long, aes(x = type, y = mean, fill = region)) +
  geom_bar(stat = "identity", position = "dodge", alpha = 0.6) +
  geom_errorbar(aes(ymin = ymin, ymax = ymax), position = position_dodge(width = 0.9), width = 0.2) +
  labs(y = "Estimated reduction (μg/m³)", x = "Pollutant Type") +
  scale_fill_manual(values = c("North" = "#EDB176", "South" = "#6959cd"), name = "Region") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, size = 22),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 22),
        axis.text.y = element_text(size = 22),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.title = element_blank(),
        legend.text = element_text(size = 22),
        plot.title = element_blank(),
        axis.line.x = element_line(size = 0.8, color = "black"), 
        axis.line.y = element_line(size = 0.8, color = "black"),
        axis.ticks = element_line(size = 0.8, color = "black"))


EV_Benefits <- read.csv("Data/EV_Benefits.csv")
EV_Benefits$diff <- EV_Benefits$NO2_EV - EV_Benefits$NO2_mean 
EV_Benefits$mid_point <- (EV_Benefits$NO2_mean + EV_Benefits$NO2_EV) / 2 
EV_Benefits$diff_pm25 <- EV_Benefits$PM25_EV - EV_Benefits$PM25_mean 
EV_Benefits$mid_point_pm25 <- (EV_Benefits$PM25_mean + EV_Benefits$PM25_EV) / 2 


#######NO2_Benefits_plot#######
NO2_Benefits_plot <- ggplot(EV_Benefits, aes(x = year_only)) +
  geom_line(aes(y = NO2_mean, color = "NO2_mean"), size = 2) +
  geom_line(aes(y = NO2_EV, color = "NO2_EV"), linetype = "solid", size = 2) +
  geom_point(aes(y = NO2_mean, color = "NO2_mean"), size = 5, shape = 21, fill = "#BA3E45") +
  geom_point(aes(y = NO2_EV, color = "NO2_EV"), size = 5, shape = 21, fill = "#FF9999") +
  geom_ribbon(aes(ymin = NO2_EV_lower, ymax = NO2_EV_upper), color = "#D69D98", alpha = 0.2) +
  annotate(
    "text", x = 2016.5-1.0, y = min(EV_Benefits$NO2_mean) + 0.8,
    label = "Average reduction by 3.28 μg/m³ (2016-2023)
              40.9% of the total decrease",
    hjust = 0, color = "black", size = 7.8
  ) +
  scale_color_manual(values = c("NO2_mean" = "#BA3E45", "NO2_EV" = "#FF9999"),
                     labels = c("NO2_mean" = expression("Predicted NO"["2"]*" level in the absence of EVs"), "NO2_EV" = expression("Actual NO"["2"]*" level")) # 自定义图例标签
  ) +
  guides(
    color = guide_legend(
      override.aes = list(
        linetype = c("solid", "solid"),
        shape = NA,                      
        size = 1.5,                    
        color = c("#BA3E45", "#FF9999")       
      ),
      keywidth = 5.5, 
      keyheight = 0.5, 
      title = NULL
    )
  ) +
  labs(
    y = expression("Average NO"["2"]*" concentration (μg/m³)")
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 26), 
    plot.title = element_text(hjust = 0.5, size = 26), 
    legend.position = c(1.03, 1.03),
    legend.justification = c(1, 1),
    legend.text.align = 0,
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(color = "black"),
    axis.text.x = element_text(color = "black", size = 26),
    axis.text.y = element_text(color = "black", size = 26),
    axis.ticks = element_line(size = 0.8, color = "black"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size = 26)
  ) +
  scale_x_continuous(breaks = seq(min(EV_Benefits$year_only), max(EV_Benefits$year_only), by = 1))



#######PM25_Benefits_plot#######
PM25_Benefits_plot <- ggplot(EV_Benefits, aes(x = year_only)) +
  geom_line(aes(y = PM25_mean, color = "PM25_mean"), size = 2) +
  geom_line(aes(y = PM25_EV, color = "PM25_EV"), linetype = "solid", size = 2) +
  geom_point(aes(y = PM25_mean, color = "PM25_mean"), size = 5, shape = 21, fill = "#00688b") +
  geom_point(aes(y = PM25_EV, color = "PM25_EV"), size = 5, shape = 21, fill = "aquamarine") +
  geom_ribbon(aes(ymin = PM25_EV_lower, ymax = PM25_EV_upper), color = "#9FBBD5", alpha = 0.2) +
  annotate(
    "text", x = 2016.5-1.5, y = min(EV_Benefits$PM25_mean) + 0.8,
    label = "Average reduction by 2.4 μg/m³ (2016-2023)
              16.7% of the total decrease",
    hjust = 0, color = "black", size = 7.8
  ) +
  scale_color_manual(values = c("PM25_mean" = "#00688b", "PM25_EV" = "aquamarine"),
                     labels = c("PM25_mean" = expression("Predicted PM"["2.5"]*" level in the absence of EVs"), "PM25_EV" = expression("Actual PM"["2.5"]*" level")) # 自定义图例标签
  ) +
  guides(
    color = guide_legend(
      override.aes = list(
        linetype = c("solid", "solid"),
        shape = NA,   
        size = 1.5,              
        color = c("#00688b", "aquamarine")    
      ),
      keywidth = 5.5, 
      keyheight = 0.5,
      title = NULL
    )
  ) +
  labs(
    y = expression("Average PM"["2.5"]*" concentration (μg/m³)")
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 26), 
    plot.title = element_text(hjust = 0.5, size = 26), 
    legend.position = c(1.02, 1.02), 
    legend.justification = c(1, 1), 
    legend.text.align = 0,
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.line = element_line(color = "black"), 
    axis.text.x = element_text(color = "black", size = 26),
    axis.text.y = element_text(color = "black", size = 26),
    axis.ticks = element_line(size = 0.8, color = "black"), 
    axis.title.x = element_blank(),
    axis.title.y = element_text(size = 26)
  ) +
  scale_x_continuous(breaks = seq(min(EV_Benefits$year_only), max(EV_Benefits$year_only), by = 1)) 





######  Combine All #######
library(cowplot)
plot_A_EVCS <- no2_map + labs(tag = expression(bold("A. estimated NO"[2]*" reduction due to 1,000 EVCS"))) + theme(plot.tag = element_text(size = 26),
                                                                                                  plot.tag.position = c(0.49, 0.96))
plot_B_EVCS <- pm25_map + labs(tag = expression(bold("B. estimated PM"["2.5"]*" reduction due to 1,000 EVCS"))) + theme(plot.tag = element_text(size = 26),
                                                    plot.tag.position = c(0.52, 0.96))

plot_evca_plot <- evca_plot + labs(tag = expression(atop(bold("C. estimated air quality"), bold(" benefits due to 1,000 EVCS")))) + theme(plot.tag = element_text(size = 26),
                                                                                                          plot.tag.position = c(0.52, 1.12))
plot_evelec_plot <- evelec_plot + labs(tag = expression(atop(bold("F. estimated air quality"), bold("benefits due to 100 MWh")))) + theme(plot.tag = element_text(size = 26),
                                                                                                          plot.tag.position = c(0.53, 1.12))

plot_A <- NO2_Benefits_plot + labs(tag = expression(atop(bold("G. annual trends in predicted NO"["2"]*" level"), bold("in the absence of EV deployment")))) + theme(plot.tag = element_text(size = 26),
                                                                                                                   plot.tag.position = c(0.4, 1.1))
plot_B <- PM25_Benefits_plot + labs(tag = expression(atop(bold("H. annual trends in predicted PM"["2.5"]*" level"), bold("in the absence of EV deployment")))) + theme(plot.tag = element_text(size = 26),
                                                                                                                      plot.tag.position = c(0.4, 1.1))
combined_plot <- ggdraw() +
  draw_plot(plot_A_EVCS, x = 0.01, y = 0.53, width = 0.35, height = 0.6) +
  draw_plot(plot_B_EVCS, x = 0.38, y = 0.53, width = 0.35, height = 0.6) +
  draw_plot(plot_E, x = 0.01, y = 0.19, width = 0.35, height = 0.6) +
  draw_plot(plot_F, x = 0.38, y = 0.19, width = 0.35, height = 0.6) +
  draw_plot(plot_A, x = 0.01, y = 0.01, width = 0.465, height = 0.26) +
  draw_plot(plot_B, x = 0.5, y = 0.01, width = 0.465, height = 0.26) +
  draw_plot(plot_evca_plot, x = 0.745, y = 0.69, width = 0.26, height = 0.26) +
  draw_plot(plot_evelec_plot, x = 0.74, y = 0.35, width = 0.26, height = 0.26) 

ggsave("Figure3.svg", plot = combined_plot, width = 22, height = 27, device = "svg")