CO2 emission in Vietnam

Visualization

Data processing

# 1. load packages -----
pacman::p_load(
   data.table,
   collapse,
   sf,
   terra,
   ggplot2,
   ggrepel,
   ggpp,
   ggpubr,
   grid,
   gridExtra
)

# 2. Data processing -----
## 2.1 Import data -----
## 2.1.1 Vietnam provinces data -----
"d:/R/rpubs_project/visualization/map_co2-emission_VN_2022/00_raw_data/gadm/gadm41_VNM_1.shx" |> 
   st_read() |> as.data.table() -> province

## 2.1.2 Vietnam population by province -----
"d:/R/rpubs_project/visualization/map_co2-emission_VN_2022/00_raw_data/pop/pop_22.xlsx" |> 
   readxl::read_xlsx() |> 
   as.data.table() |> 
   province[i = _, on = .(VARNAME_1 == province)] -> 
   province

## 2.1.3 CO2 emission -----
"d:/R/rpubs_project/visualization/map_co2-emission_VN_2022/00_raw_data/co2/v8.0_FT2022_GHG_CO2_2022_TOTALS_emi.nc" |> 
   rast() -> co2_sp

## 2.2 Finalize dataset -----
# set crs of co2 matching province
st_crs(province$geometry) # EPSG: 4326
crs(co2_sp) <- "epsg:4326"

# Extract co2 data for each province by summation
exactextractr::exact_extract(co2_sp, province$geometry, "sum") -> province$co2_sum

# Calculate CO2 per capita
province[, co2_pc := co2_sum / (pop * 1000)][]

show_save <- function(plot){
   ggsave("d:/R/rpubs_project/visualization/map_co2-emission_VN_2022/03_plots/co2_map.png", plot = plot, height = 9, width = 8, dpi = 300,
          device = ragg::agg_png())
   file.show("d:/R/rpubs_project/visualization/map_co2-emission_VN_2022/03_plots/co2_map.png")
}

# 3. Visualization -----
## 3.1 Setup -----

# Color palette
cols <- hcl.colors(6, "Inferno", rev = T)[1:4]
pal <- colorRampPalette(cols)(64)

province |> st_as_sf() |> 
   sbt(co2_pc > 8, VARNAME_1) |> 
   st_centroid() |> 
   (
      \(.) data.table(.$VARNAME_1, st_coordinates(.))
   )() |> setnames(new = .c(province, x, y)) -> top5_co2_max

province |> as.data.table() |> 
   _[, top_co2_pc_max := 
        fcase(co2_pc > 8, VARNAME_1,
              co2_pc < 8, "")][] -> 
   province

province |> as.data.table() |> 
   _[, top_co2_sum_max := 
        fcase(co2_sum > 14 * 10^6, VARNAME_1,
              co2_sum < 14 * 10^6, "")][] -> 
   province

province |> tfm(co2_sum_mil = co2_sum / 10^6) -> province

province |> st_as_sf() -> province

## 3.2 Mapping -----

### 3.2.1 CO2 emission per capita by province -----



province |> 
   ggplot() +
      geom_sf(
         aes(fill = co2_pc),
         color = "white",
         size = .15
      ) +
      scale_fill_gradientn(
         name = "",
         colors = pal
         # limits = c(0, 35),
         # breaks = seq(5, 35, 5),
         # labels = c("<5", "5-10", "10-15", "15-20", "20-25", "25-30", ">30"),
         # guide = guide_colorbar(
         #    barheight = 9,
         #    barwidth = 1
         # )
      ) +
      geom_text_repel(aes(label = top_co2_pc_max,
                          geometry = geometry),
                      stat = "sf_coordinates",
                      min.segment.length = 0,
                      nudge_x = 1.5,
                      color = "gray20",     # text color
                      segment.color = "gray30",
                      family = "Merriweather Sans",
                      size = 3.3
                      ) +
      # labs(title = "Per capita CO₂ emissions by province, 2022") +
      labs(caption = "Per capita by province (tons)") +
      theme_void() +
      theme(
         legend.position = c(.2, .5),
         legend.text = element_text(
            color = "gray20",
            vjust = 1.5, size = 10,
            family = "Work Sans"
         ),
         plot.margin = unit(c(
            t = 1, b = 0, l = 0, r = 0
         ), "lines"),
         plot.caption = element_text(
            family =  "Merriweather Sans",
            size = 16, hjust = .5
         )
      ) -> p1

### 3.2.2 CO2 emission by province -----
province |> 
      ggplot() +
      geom_sf(
         aes(fill = co2_sum_mil),
         color = "white",
         size = .15
      ) +
      scale_fill_gradientn(
         name = "",
         colors = pal,
         limits = c(0, 45),
         breaks = seq(5, 45, 5),
         # labels = c("<5", "5-10", "10-15", "15-20", "20-25", "25-30", "30-35", "35-40", ">40"),
         guide = guide_colorbar(
            barheight = .5,
            barwidth = 38
         )
      ) +
      geom_text_repel(aes(label = top_co2_sum_max,
                          geometry = geometry),
                      stat = "sf_coordinates",
                      min.segment.length = 0,
                      nudge_x = 1.5,
                      color = "gray20",     # text color
                      segment.color = "gray30",
                      family = "Merriweather Sans",
                      size = 3.3
      ) +
      labs(captions = "Total by province (mil. tons)") +
      theme_void() +
      theme(
         legend.position = "bottom",
         legend.text = element_text(
            color = "gray20",
            size = 13,
            family = "Work Sans",
            hjust = 1
         ), legend.ticks = element_line(
            linewidth = 2,
            colour = "white"),
         plot.margin = unit(c(
            t = 1, b = 0, l = 0, r = 0
         ), "lines"),
         plot.caption = element_text(
            family =  "Merriweather Sans",
            size = 16, hjust =.5
         )
      ) -> p2

share_leg <- get_legend(p2)



ggarrange(
   p2 + theme(legend.position = "none",
              plot.title = element_blank()),
   p1 + theme(legend.position = "none",
              plot.title = element_blank()),
   nrow = 1
   # legend.grob = share_leg,
   # legend = 'bottom'
) |> ggarrange(
      ..1 = text_grob(
         label = "CO2 emission in Vietnam, 2022",
         family = "Merriweather Sans ExtraBold",
         face = "bold",
         size = 30, just = "left", x = .04,
         color = "gray20") |> 
         as_ggplot(),
      ..2 = text_grob(
         label = "Top 5 provinces that produced the most CO2 indicated at each plot", 
         size = 16, just = "left", x = .04, 
         color = "gray30", y = .6,
         family = "Merriweather Sans"
      ),
      ..3 = _, 
      ..4 = as_ggplot(share_leg),
      nrow = 4,
      heights = c(.075, .035, .83, .07)
   ) |> grid.arrange(
      bottom = textGrob(
         "By: Duy Khanh",
         gp = gpar(family = "Libre Baskerville", 
                   fontsize = 12),
         x = .98,
         y = .5,
         just = "right"
      )
   ) |> as_ggplot() +
   annotation_custom(
      rectGrob(gp = gpar(fill = "firebrick")),
      xmin = 0, xmax = .02,
      ymin = .905, ymax = .98
   ) -> my_map
show_save(my_map)