Administrative Map of Vietnam + Spratly and Paracel Islands

Daily Graph Series

Introduction

Choropleth Maps là một trong những công cụ hình ảnh hóa dữ liệu hữu ích và được sử dụng phổ biến. Sử dụng công cụ hình ảnh này có truyền tải một cách nhanh chóng và ấn tượng, ví dụ, thị phần của hãng Coca Cola tại các tỉnh ở VN ra sao hoặc tỉ lệ đói nghèo theo khu vực địa lí bằng cách căn cứ vào độ đậm hay nhạt của màu sắc. Trong bài viết này chúng ta sẽ thực hiện một số lệnh thực hành tạo ra bản đồ hành chính cũng như Choropleth Map với dữ liệu cho Việt Nam.

Administrative Map of Vietnam

Để vẽ bản đồ hành chính của Việt Nam trước hết chúng ta cần có Geo-spatial Data. Hiện R có nhiều package hỗ trợ cho việc download và lấy loại dữ liệu này một cách thuận tiện ví dụ như raster, sf. Tuy nhiên các packages này, vì nhiều lí do, sẽ chưa có dữ liệu của hai quần đảo Trường Sa (Spratly Islands) và Hoàng Sa (Paracel Islands) thuộc Việt Nam. Do vậy chúng ta có thể sử dụng các nguồn dữ liệu Geo-spatial Data cho hai quần đảo này được cung cấp bởi các tổ chức và cá nhân ở Việt Nam như dưới đây:

# Clear workspace: 
rm(list = ls())

# Load tidyverse package: 

library(tidyverse)
 
# Get geo-spatial data for Vietnam provided: 

link2 <- "https://data.opendevelopmentmekong.net/dataset/999c96d8-fae0-4b82-9a2b-e481f6f50e12/resource/2818c2c5-e9c3-440b-a9b8-3029d7298065/download/diaphantinhenglish.geojson?fbclid=IwAR1coUVLkuEoJRsgaH81q6ocz1nVeGBirqpKRBN8WWxXQIJREUL1buFi1eE"

link1 <- "https://raw.githubusercontent.com/nguyenduy1133/data/main/Dia_phan_Tinh_cap_nhat.geojson?fbclid=IwAR1cPIIswZ9y8ZvsmW8ioqh57malxnlrBr2L3EjYRfMPiDjCqu4kCA2PUxQ"

vn_spatial <- sf::st_read(link1)

## Reading layer `Dia_phan_Tinh_cap_nhat' from data source `https://raw.githubusercontent.com/nguyenduy1133/data/main/Dia_phan_Tinh_cap_nhat.geojson?fbclid=IwAR1cPIIswZ9y8ZvsmW8ioqh57malxnlrBr2L3EjYRfMPiDjCqu4kCA2PUxQ' using driver `GeoJSON'
## Simple feature collection with 63 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 102.1421 ymin: 6.953306 xmax: 116.9473 ymax: 23.3939
## geographic CRS: WGS 84

# Rename for Dak Lak province: 

vn_spatial %>% mutate(Name = case_when(Name == "Dak Lak" ~ "Dac Lac", TRUE ~ Name)) -> vn_spatial

# Make a draft of map: 

ggplot() + 
  geom_sf(data = vn_spatial) + 
  labs(title = "Vietnam Map - Version 1")

Hoặc theo một cách khác:

ggplot() + 
  geom_sf(data = vn_spatial, aes(fill = Name), show.legend = FALSE) + 
  labs(title = "Vietnam Map - Version 2")

Trong nhiều tình huống chúng ta cần bản đồ mà thể hiện cả cấp vùng (Region) hoặc tiểu vùng (sub-region). Trước hết cần lấy dữ liệu về vùng/tiểu vùng (vốn rất nhiều nguồn) như sau:

# Scrap list of provinces by region and sub-region level: 

library(rvest)

# Collect region/sub-region data from Wiki: 

provinces <- "https://en.wikipedia.org/wiki/Provinces_of_Vietnam"

provinces %>% 
  read_html() %>% 
  html_nodes(xpath = '//*[@id="mw-content-text"]/div[1]/table[6]') %>% 
  html_table() %>% 
  .[[1]] -> provinces_vn

# Function extracts data in table: 

extract_table <- function(i) {
  
  provinces_vn %>% 
    slice(i) %>% 
    pull(3) %>% 
    str_split("\n", simplify = TRUE) %>% 
    str_replace_all("†", "") %>% 
    str_squish() %>% 
    as.vector() -> province_names
  
  provinces_vn %>% 
    slice(i) %>% 
    pull(1) %>% 
    str_split("\\(", simplify = TRUE) %>% 
    str_split("\\,", simplify = TRUE) %>% 
    str_replace_all("Vietnam", "") %>% 
    as.vector() %>% 
    str_squish() -> region
  
  provinces_vn %>% 
    slice(i) %>% 
    pull(2) %>% 
    str_split("\\(", simplify = TRUE) %>% 
    str_replace_all("\\)", "") %>% 
    str_replace_all("Vietnam", "") %>% 
    as.vector() %>% 
    str_squish() -> sub_region 
  
  tibble(province = province_names, region_vn = region[2], region_en = region[1], 
         sub_region_vn = sub_region[2], sub_region_en = sub_region[1]) -> df_final
  
  return(df_final)
  
}

# Use the function: 

lapply(1:nrow(provinces_vn), extract_table) -> province_region

do.call("bind_rows", province_region) -> province_region_vietnam

# Rename for some provinnces: 

library(stringi)

province_region_vietnam %>% 
  mutate(province_latin = stri_trans_general(province, "Latin-ASCII")) %>% 
  mutate(province_latin = case_when(province_latin == "Thua Thien-Hue" ~ "Thua Thien - Hue", 
                                    province_latin == "Ba Ria-Vung Tau" ~ "Ba Ria - Vung Tau", 
                                    province_latin == "Ho Chi Minh City" ~ "TP. Ho Chi Minh", 
                                    TRUE ~ province_latin)) -> province_region_vietnam

Sau đó join hai bộ dữ liệu (Geo-spatial and Region Datasets) lại và mapping:

# Combine the two datasets: 

full_join(vn_spatial, province_region_vietnam, by = c("Name" = "province_latin")) -> vn_spatial_region

# Administrative map of Vietnam: 

ggplot() + 
  geom_sf(data = vn_spatial_region, aes(fill = sub_region_en)) + 
  labs(title = "Administrative Regions of Viet Nam")

Hoặc làm cho bản đồ trở nên đẹp hơn:

# Elegant map of Vietnam: 

my_colors <- c('#e41a1c','#377eb8','#4daf4a','#984ea3','#ff7f00','#ffff33','#a65628','#f781bf')
bgr_color <- "azure"

library(extrafont)
my_font <- "Roboto Condensed"

ggplot() + 
  geom_sf(data = vn_spatial_region, aes(fill = sub_region_en)) + 
  scale_fill_manual(values = my_colors, name = "Region:") + 
  theme(legend.position = c(0.82, 0.83)) + 
  theme(text = element_text(family = my_font)) + 
  theme(legend.key.size = unit(0.5, "cm"), legend.key.width = unit(0.5, "cm") ) + 
  theme(legend.text = element_text(family = my_font, color = "grey20")) + 
  theme(panel.background = element_rect(fill = bgr_color)) +
  theme(plot.background = element_rect(fill = bgr_color)) + 
  theme(legend.background = element_rect(fill = bgr_color)) + 
  theme(plot.margin = unit(rep(0.7, 4), "cm")) + 
  theme(plot.title = element_text(size = 16)) + 
  labs(title = "Administrative Regions of Viet Nam", 
       caption = "Data Source: https://data.opendevelopmentmekong.net")

Choropleth Map

Trước hết lấy dữ liệu về mật độ dân số theo cả cấp Huyện và Xã cho các đơn vị hành chính của Việt Nam:

#================================
#  Stage 1: Extract all links
#================================

all_links <- "https://www.citypopulation.de/Vietnam.html"
pg <- read_html(all_links)
m <- html_nodes(pg, "a")
k <- html_attr(m, "href")
all_links_communes_level <-  str_c("https://www.citypopulation.de/en/vietnam/", k[-c(1:6)])

#============================================================
#  Stage 2: Function for getting popolation by commune level
#============================================================

get_pop_data <- function(link) {
  
  link %>% 
    read_html() %>% 
    html_nodes(xpath = '//*[@id="tl"]') %>% 
    html_table(fill = TRUE) %>% 
    .[[1]] -> df
  
  df %>%   
    select(1:3) -> df
  
  names(df) <- c("Name", "Status", "Pop")
  
  df %>% 
    mutate(Name_Latin = stringi::stri_trans_general(Name, "Latin-ASCII"), 
           Pop = str_replace_all(Pop, "[^0-9]", "") %>% as.integer()) %>% 
    return()
  
}


# Use the function: 

lapply(all_links_communes_level, get_pop_data) -> pop_cummunes_list

do.call("bind_rows", pop_cummunes_list) -> pop_cummunes_df

# Show some observations: 

library(knitr)

pop_cummunes_df %>% 
  select(-1) %>% 
  filter(!str_detect(Status, "District")) %>% 
  head() %>% 
  kable(caption = "Table 1: Population by Commune")

Table 1: Population by Commune

Status	Pop	Name_Latin
Ward	18499	An Thoi
Ward	18307	Binh Thuy
Ward	11745	Bui Huu Nghia
Ward	16450	Long Hoa
Ward	15232	Long Tuyen
Ward	10773	Thoi An Dong

Rồi thực hiện một số thao tác xử lí dữ liệu và lấy Geo-spatial Data bằng hàm raster::getData():

# Population by district level: 

commune_df_pop <- pop_cummunes_df %>% filter(str_detect(Status, "District"))

# Get geo-spatial data by district level: 

vietnam_dis <- raster::getData("GADM", country = "Vietnam", level = 2)

# Convert to data frame: 

vietnam_df <- vietnam_dis %>% fortify(region = "NAME_2")

vietnam_df %>% 
  mutate(dis_names = stri_trans_general(id, "Latin-ASCII")) -> vietnam_df

# Join the two data set: 

full_join(vietnam_df, commune_df_pop, by = c("dis_names" = "Name_Latin")) -> df_pop_dis_vn

df_pop_dis_vn %>% mutate(Pop = Pop / 1000) -> df_pop_dis_vn

Và vẽ:

# Choropleth map: 

special_provinces <- c("Khanh Hoa", "Da Nang")
library(viridis)

ggplot() +
  geom_sf(data = vn_spatial_region %>% filter(Name %in% special_provinces)) +  
  geom_polygon(data = df_pop_dis_vn, aes(fill = Pop, x = long, y = lat, group = group), color = "grey80") + 
  theme(legend.position = c(0.8, 0.9)) + 
  theme(text = element_text(family = my_font)) + 
  theme(legend.key.size = unit(0.5, "cm"), legend.key.width = unit(0.5, "cm") ) + 
  theme(legend.text = element_text(family = my_font, color = "grey20")) + 
  theme(legend.title = element_text(family = my_font, color = "black", size = 11)) + 
  theme(panel.background = element_rect(fill = bgr_color)) +
  theme(plot.background = element_rect(fill = bgr_color)) + 
  theme(legend.background = element_rect(fill = bgr_color)) + 
  theme(axis.text = element_blank()) + 
  theme(axis.title = element_blank()) + 
  theme(axis.ticks = element_blank()) + 
  theme(plot.margin = unit(rep(0.7, 4), "cm")) + 
  theme(plot.title = element_text(size = 16)) + 
  labs(title = "Population Density (in thousands) by District for Vietnam", caption = "Data Source: https://www.citypopulation.de/Vietnam.html") + 
  annotate("text", x = 112, y = 18, label = "Paracel Islands\n(Da Nang City)", family = my_font, size = 3, fontface = "italic") +  
  annotate("text", x = 115.5, y = 12.5, label = "Spratly Islands\n(Khanh Hoa)", family = my_font, size = 3, fontface = "italic") + 
  scale_fill_viridis(direction = -1, 
                     option = "C", 
                     name = "Density:", 
                     guide = guide_colourbar(direction = "horizontal",
                                             barheight = unit(3, units = "mm"),
                                             barwidth = unit(30, units = "mm"),
                                             title.hjust = 0.5,
                                             label.hjust = 0.5, 
                                             limits = c(0, 800), 
                                             title.position = "top"))

District and Commune Level

Dưới đây là R codes vẽ map cho tỉnh Nghệ An theo cấp huyện:

nghean <- vietnam_dis[vietnam_dis$NAME_1 == "Nghệ An", ]

nghean_df <- nghean %>% 
  fortify(region = "NAME_2") %>% 
  mutate(id = stri_trans_general(id, "Latin-ASCII"))

ggplot() +
  geom_polygon(data = nghean_df,
               aes(fill = id, x = long, y = lat, group = group), 
               color = "gray69") + 
  theme(text = element_text(family = my_font)) + 
  theme(plot.margin = unit(rep(0.7, 4), "cm")) + 
  labs(title = "Nghe An Province by District", 
       caption = "Data Source: http://diva-gis.org/gdata")

Theo cấp xã (rất nhiều xã nên bản đồ vẽ kiểu này không có nhiều ý nghĩa):

# Get geo-spatial data by commune level: 

vietnam_com <- raster::getData("GADM", country = "Vietnam", level = 3)

nghean_com <- vietnam_com[vietnam_com$NAME_1 == "Nghệ An", ]

nghean_df_com <- nghean_com %>% 
  fortify(region = "NAME_3") %>% 
  mutate(id = stri_trans_general(id, "Latin-ASCII"))

ggplot() +
  geom_polygon(data = nghean_df_com,
               aes(fill = id, x = long, y = lat, group = group), 
               color = "gray69", show.legend = FALSE) + 
  theme(text = element_text(family = my_font)) + 
  theme(plot.margin = unit(rep(0.7, 4), "cm")) + 
  labs(title = "Nghe An Province by Commune", 
       caption = "Data Source: http://diva-gis.org/gdata")

Conclusion

R là một công cụ mạnh và linh hoạt cho Data Visualization. Mặc dù không phải là một tool chuyên về Mapping và GIS nhưng R cũng đủ mạnh để có thể thực hiện hầu hết các tình huống về về mapping như trong post này đã chỉ ra.

References

Các cases trình bày trong post này nhất định là không thể cover hết tất cả các tình huống về mapping trong thực tế. Do vậy các nguồn thao khảo sau là hữu ích:

1. https://mikegruz.github.io/articles/2017-03/mapping-in-R
2. https://www.r-spatial.org/r/2018/10/25/ggplot2-sf-3.html
3. https://www.jessesadler.com/post/simple-feature-objects/
4. https://cengel.github.io/R-spatial/intro.html
5. https://dominicroye.github.io/en/2018/access-to-climate-reanalysis-data-from-r/
6. https://geocompr.robinlovelace.net/geometric-operations.html
7. https://www.r-spatial.org/r/2018/10/25/ggplot2-sf-2.html
8. https://github.com/ropensci/geojsonio
9. https://geocompr.robinlovelace.net/attr.html