How climate changed over the past decade in Vietnam ?

Data visualization series

Motivation

This post was initialized from my graduate dissertation about the impact of climate change on Vietnamese farmers’ income. At the time I was doing it, I could not explore the topic by visualization approach because of some limitations. Before digging into the relationship between farmers’ income and climate change, I would take a look of how climate in Vietnam has changed over the past decade at different levels of administrative boundaries by visualization.

All data used in this post is avaiable and open source. Monthly climte data was downloaded from WorldClim, variables were average minimum temperature (°C), average maximum temperature (°C) and total precipitation (mm) from 2010 to 2020. Geographic data of Vietnam was downloaded from GADM. I prefer to apply animation to the plots so as to get an intuitive sense of the change.

Disclaimer: The maps used on this blog are for illustrative purposes and do not express my implication or my opinion, concerning the legal status of any country or territory or concerning the delimitation of frontiers or boundaries.

Visualization of climate change in Vietnam

Temparature

R Codes for data processing and visualization

Mapping code part was inspired by this Milos’s tutorial. I was trying to employ the gganimate package to add transition of the time. Unfortunately, my laptop was not able to render such a long time, like a decade, to illustrate the time continuously. In the future, I will try with other approaches, it could be using a stronger computer or other transition types of gganimate. Therefore, the plot and codes here are still room for improvement.

# 1. Load packages -----
pacman::p_load(
   data.table, # data wrangling
   collapse, # data wragling
   sf, # geospatial data wrangling
   terra, # spatial data wrangling
   # raster,
   stringr,
   lubridate,
   # hrbrthemes,
   ggplot2, # visualization
   gganimate,
   # ggspatial,
   gifski,
   extrafont
   # showtext
)


# 2. Data processing -----
## 2.1 Geographic data -----

### 2.1.1 Import data -----
st_read("d:/R/rpubs_project/visualization/map_climate-change-VN/00_raw_data/geography/gadm41_VNM_0.shp") -> geom_dt

# Varname to lower case
setnames(geom_dt, new = str_to_lower(names(geom_dt)))


## 2.2 CLimate data -----

### 2.2.1 Average maximum temp -----

#### 2.2.1.1 import data -----
list.files(c("d:/R/rpubs_project/visualization/map_climate-change-VN/00_raw_data/climate/wc2.1_cruts4.06_2.5m_tmax_2010-2019/",
             "d:/R/rpubs_project/visualization/map_climate-change-VN/00_raw_data/climate/wc2.1_cruts4.06_2.5m_tmax_2020-2021/"), 
           full.names = T) |> 
   str_subset("(tmax_2021-\\d{2})(?=\\.tif)", negate = T) |> 
   str_subset("(tmax_2020-0[2-9])(?=\\.tif)", negate = T) |> 
   str_subset("(tmax_2020-1(0|1|2))(?=\\.tif)", negate = T) |> 
   (
      \(.) future.apply::future_lapply(., rast) |> setNames(nm = str_extract(., "(tmax_\\d{4}-\\d{2})(?=\\.)"))
   )() -> t_max_sr_10_20_list

mapply(\(x, y) `names<-`(y, x), 
       names(t_max_sr_10_20_list), t_max_sr_10_20_list) -> 
   t_max_sr_10_20_list

# whether raster and vector objects are identical
all.equal(st_crs(geom_dt) , st_crs(t_max_sr_10_20_list[[1]]))

# transform crs of geom_dt that matching t_max SpatRaster
geom_dt <- st_transform(geom_dt, t_max_sr_10_20_list[[1]] |> st_crs())

# quick ploting 1 layer of t_max spatraster list
# t_max_sr_10_20_list[[1]] |> 
#    as.data.frame(xy = T, na.rm = T) |> 
#    ggplot() +
#    geom_tile(aes(x = x, y = y, fill = wc2.1_2.5m_tmax_01)) +
#    coord_sf() +
#    scale_fill_viridis_c(option = "magma", direction = -1) +
#    theme_void() +
#    theme(legend.justification.right = c(0, .8),
#          legend.title = element_blank())

# cropping global raster to local raster by specified vector (reduced rectangular size)
future.apply::future_lapply(t_max_sr_10_20_list, 
       \(.) crop(., geom_dt) |> 
          # masking local raster by vector shap
          mask(geom_dt) |> 
          as.data.frame(xy = T, na.rm = T) |> 
          as.data.table() |> 
          # transform to long format
          melt(id.vars = c("x", "y"))) |> 
   rbindlist() -> t_max_sr_10_20_dt

t_max_sr_10_20_dt[, layer := variable][, variable := NULL][]

# create date from layer name
t_max_sr_10_20_dt[, date := str_extract(layer, "\\d{4}-\\d{2}") |> 
                     str_c("-01") |> as.IDate()][]

t_max_sr_10_20_dt[, month := date |> month(label = T)] |> 
   _[, year := date |> year()] |> 
   _[, month_year := paste0(month, ", ", year)][] -> 
   t_max_sr_10_20_dt

#### 2.2.1.3 breaks -----
min_value <- t_max_sr_10_20_dt[, min(value)]
min_value
max_value <- t_max_sr_10_20_dt[, max(value)]
max_value


breaks <- classInt::classIntervals(
   t_max_sr_10_20_dt$value,
   n = 14,
   style = "pretty"
)$brks

breaks

#### 2.2.1.5 color -----
cols <- colorRampPalette(rev(RColorBrewer::brewer.pal(11, "Spectral")))

#### 2.2.16 font -----

font_import()
loadfonts()


# setnames(fonts_dt, new = str_to_lower(names(fonts_dt)))
# 
# fonts_dt[, familyname |> unique()]



### 2.2.2 map -----
make_plot <- function(){
   t_max_sr_10_20_dt |> split(t_max_sr_10_20_dt$date) |>
      future.apply::future_lapply(\(data) {
         map_vn <-
            ggplot(data) +
            # ggplot(t_max_sr_10_20_dt[date == "2010-01-01"]) +
            geom_tile(aes(x = x, y = y, fill = value)) +
            # facet_wrap(~ date) +
            scale_fill_gradientn(
               name = "\u00B0C",
               colors = cols(14),
               limits = c(5, 40),
               breaks = seq(5, 40, 5),
               guide = guide_colorbar(
                  barheight = 6,
                  barwidth = .8,
                  nbin = 8,
                  direction = "vertical"
               )
            ) +
            coord_sf(crs = crs(t_max_sr_10_20_list[[1]])) +
            annotate("text", x = 108.1, y = 19, label = data$month_year |> unique(),
               size = 4, family = "Libre Baskerville") +
            theme_minimal() +
            theme(
               # text = element_text(family = "Josefin Sans"),
               axis.line = element_blank(),
               axis.title.x = element_blank(),
               axis.title.y = element_blank(),
               axis.text.x = element_blank(),
               axis.text.y = element_blank(),
               axis.ticks = element_blank(),
               legend.position = c(.2, .45),
               # legend.title = element_blank(),
               legend.text = element_text(
                  size = 8, color = "grey10",
                  family = "Merriweather Sans"
               ),
               plot.title = element_text(
                  size = 32, color = "gray10",
                  hjust = .5, vjust = 0, 
                  family = "Merriweather Sans Medium"
               ),
               plot.subtitle = element_text(
                  size = 16, color = "#c43c4e",
                  hjust = .5, vjust = -1.5,
                  family = "Merriweather Sans"
               ),
               plot.caption = element_text(
                  family = "Libre Baskerville"
               ),
               panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               plot.margin = unit(c(t= 1, r = 0, l = 0, b = 0), "lines")
            ) +
            labs(
               x = "",
               y = "",
               title = "Climate change in Vietnam",
               subtitle = "Monthly maximum temperature (2010 - 2020)",
               caption = "Author: Duy Khanh\nSource: WorldClim"
            )
         print(map_vn)
         }
      )
   }
# ggsave("d:/R/rpubs_project/visualization/map_climate-change-VN/03_plots/temp_vn.png", map_vn,
#        height = 8,
#        width = 7,
#        units = "in",
#        device = ragg::agg_png)
   
# timelapse_vn_map <- map_vn +
#    annotate("text", x = 108, y = 19, label = "{frame_time}",
#          size = 4.5, family = "Josefin Sans") +
#    transition_time(date) +
#    enter_fade() +
#    exit_fade() +
#    ease_aes("linear", interval = .2)
# 
# animated_time_map <- gganimate::animate(
#    timelapse_vn_map,
#    nframes = 121 * 2,
#    # duration = 20,
#    start_pause = 5,
#    end_pause = 30,
#    height = 8,
#    width = 7,
#    units = "in",
#    res = 200
#    # fps = 15,
#    # renderer = gifski_renderer(loop = T)
# )
# gganimate::anim_save(
#    "d:/R/rpubs_project/visualization/map_climate-change-VN/03_plots/vn_temp.gif",
#    animated_time_map
# )

save_gif(make_plot(), "d:/R/rpubs_project/visualization/map_climate-change-VN/03_plots/vn_temp.gif", 
         height = 1100, width = 1100, delay = .5, res = 150)

Precipitation

R codes for data processing and visualization

rm(list = ls())

# 2. Data processing -----
## 2.1 Geographic data -----

### 2.1.1 Import data -----
# st_read("d:/R/rpubs_project/visualization/map_climate-change-VN/00_raw_data/geography/gadm41_VNM_0.shp") -> geom_dt
# 
# # Varname to lower case
# setnames(geom_dt, new = str_to_lower(names(geom_dt)))


## 2.2 CLimate data -----

### 2.2.1 Total precipitation -----

#### 2.2.1.1 import data -----
list.files(c("d:/R/rpubs_project/visualization/map_climate-change-VN/00_raw_data/climate/wc2.1_cruts4.06_2.5m_prec_2010-2019/",
             "d:/R/rpubs_project/visualization/map_climate-change-VN/00_raw_data/climate/wc2.1_cruts4.06_2.5m_prec_2020-2021/"), 
           full.names = T) |> 
   str_subset("(prec_2021-\\d{2})(?=\\.tif)", negate = T) |> 
   str_subset("(prec_2020-0[2-9])(?=\\.tif)", negate = T) |> 
   str_subset("(prec_2020-1(0|1|2))(?=\\.tif)", negate = T) |> 
   (
      \(.) future.apply::future_lapply(., rast) |> setNames(nm = str_extract(., "(prec_\\d{4}-\\d{2})(?=\\.)"))
   )() -> prec_sr_10_20_list

mapply(\(x, y) `names<-`(y, x), 
       names(prec_sr_10_20_list), prec_sr_10_20_list) -> 
   prec_sr_10_20_list

# whether raster and vector objects are identical
all.equal(st_crs(geom_dt) , st_crs(prec_sr_10_20_list[[1]]))

# transform crs of geom_dt that matching prec SpatRaster
geom_dt <- st_transform(geom_dt, prec_sr_10_20_list[[1]] |> st_crs())

# quick ploting 1 layer of prec spatraster list
# prec_sr_10_20_list[[1]] |> 
#    as.data.frame(xy = T, na.rm = T) |> 
#    ggplot() +
#    geom_tile(aes(x = x, y = y, fill = wc2.1_2.5m_prec_01)) +
#    coord_sf() +
#    scale_fill_viridis_c(option = "magma", direction = -1) +
#    theme_void() +
#    theme(legend.justification.right = c(0, .8),
#          legend.title = element_blank())

# cropping global raster to local raster by specified vector (reduced rectangular size)
future.apply::future_lapply(prec_sr_10_20_list, 
       \(.) crop(., geom_dt) |> 
          # masking local raster by vector shap
          mask(geom_dt) |> 
          as.data.frame(xy = T, na.rm = T) |> 
          as.data.table() |> 
          # transform to long format
          melt(id.vars = c("x", "y"))) |> 
   rbindlist() -> prec_sr_10_20_dt

prec_sr_10_20_dt[, layer := variable][, variable := NULL][]

# create date from layer name
prec_sr_10_20_dt[, date := str_extract(layer, "\\d{4}-\\d{2}") |> 
                     str_c("-01") |> as.IDate()][]

prec_sr_10_20_dt[, month := date |> month(label = T)] |> 
   _[, year := date |> year()] |> 
   _[, month_year := paste0(month, ", ", year)][] -> 
   prec_sr_10_20_dt

#### 2.2.1.3 breaks -----
min_value <- prec_sr_10_20_dt[, min(value)]
min_value
max_value <- prec_sr_10_20_dt[, max(value)]
max_value


breaks <- classInt::classIntervals(
   prec_sr_10_20_dt$value,
   n = 14,
   style = "pretty"
)$brks

breaks

#### 2.2.1.5 color -----
cols <- colorRampPalette(RColorBrewer::brewer.pal(9, "Blues")[3:9])

#### 2.2.16 font -----

# font_import()
# loadfonts()


# setnames(fonts_dt, new = str_to_lower(names(fonts_dt)))
# 
# fonts_dt[, familyname |> unique()]



### 2.2.2 map -----
make_plot <- function(){
   prec_sr_10_20_dt |> split(prec_sr_10_20_dt$date) |>
      future.apply::future_lapply(\(data) {
         map_vn <-
            ggplot(data) +
            # ggplot(prec_sr_10_20_dt[date == "2010-01-01"]) +
            geom_tile(aes(x = x, y = y, fill = value)) +
            # facet_wrap(~ date) +
            scale_fill_gradientn(
               name = "mm",
               colors = cols(16),
               limits = c(0, 2200),
               # breaks = breaks,
               guide = guide_colorbar(
                  barheight = 6,
                  barwidth = .8,
                  nbin = 8,
                  direction = "vertical"
               )
            ) +
            coord_sf(crs = crs(prec_sr_10_20_list[[1]])) +
            annotate("text", x = 108.1, y = 19, 
                     label = 
                        # "Jan, 2010",
                        data$month_year |> unique(),
               size = 4, family = "Libre Baskerville") +
            theme_minimal() +
            theme(
               # text = element_text(family = "Josefin Sans"),
               axis.line = element_blank(),
               axis.title.x = element_blank(),
               axis.title.y = element_blank(),
               axis.text.x = element_blank(),
               axis.text.y = element_blank(),
               axis.ticks = element_blank(),
               legend.position = c(.2, .45),
               # legend.title = element_blank(),
               legend.text = element_text(
                  size = 8, color = "grey10",
                  family = "Merriweather Sans"
               ),
               plot.title = element_text(
                  size = 32, color = "gray10",
                  hjust = .5, vjust = 0, 
                  family = "Merriweather Sans Medium"
               ),
               plot.subtitle = element_text(
                  size = 16, color = "#08306B",
                  hjust = .5, vjust = -1.5,
                  family = "Merriweather Sans"
               ),
               plot.caption = element_text(
                  family = "Libre Baskerville"
               ),
               panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               plot.margin = unit(c(t= 1, r = 0, l = 0, b = 0), "lines")
            ) +
            labs(
               x = "",
               y = "",
               title = "Climate change in Vietnam",
               subtitle = "Monthly total precipitation (2010 - 2020)",
               caption = "Author: Duy Khanh\nSource: WorldClim"
            )
         print(map_vn)
         }
      )
}

save_gif(make_plot(), "d:/R/rpubs_project/visualization/map_climate-change-VN/03_plots/vn_prec.gif", 
         height = 1100, width = 1100, delay = .5, res = 150)