Motivations

Hình ảnh hóa dữ liệu (Data Visualization) nói chung và hình ảnh hóa dữ liệu không gian (Spatial Data Visualization) thường là phần không thể thiếu trong các báo cáo cũng như nghiên cứu liên quan đến các vấn đề kinh tế - xã hội. Dưới đây là một vài ví dụ điển hình:

Có nhiều tools có thể được sử dụng cho hình ảnh hóa dữ liệu không gian. Nhiều trong số đó là những phần mềm chuyên biệt có bản quyền và rất đắt. Khóa học này được thiết kế để chuẩn bị cho người học nắm được và vận dụng R cho việc hình ảnh hóa dữ liệu không gian.

Objectives

Khóa học này trang bị và hướng dẫn cho người học vẽ được hầu hết các kiểu bản đồ thường xuất hiện trong các báo cáo và nghiên cứu liên quan đến các vấn đề kinh tế - xã hội bằng cách chỉ sử dụng ggplot2.

Prerequisite

Khóa học này phù hợp với những người học mà: (1) sử dụng R ở mức cơ bản trở lên, (2) tương đối quen thuộc với ggplot2.

Preparation

Trước hết cần cài đặt:

Sau khi cài đặt xong cần đảm bảo rằng những thứ cần cài đạt đã xong bằng cách test đoạn mã sau:

vietnam_province <- raster::getData("GADM", country = "Vietnam", level = 1)

Nếu chạy đoạn mã trên mà báo lỗi - chẳng hạn lỗi cài đặt raster, thì cài đặt lại package này.

The Grammar of Graphics

Phần này nhắc lại ngữ pháp hình ảnh (The Grammar of Graphics) của ggplot2. Trước hết tạo một bộ dữ liệu cho mục đích minh họa:

# Clear R environment: 

rm(list = ls())

# Create data for our illutration: 

df1 <- data.frame(point_name = c("A", "B", "C", "D"), 
                  x_value = c(1, 2, 3, 1), 
                  y_value = c(1, 2, 1, 1))

# Show data: 

df1
##   point_name x_value y_value
## 1          A       1       1
## 2          B       2       2
## 3          C       3       1
## 4          D       1       1

Vẽ một tam giác màu xanh với đường viền màu đỏ:

library(dplyr) # For data cleansing/wrangling. 
library(ggplot2) # For data visualization. 

# Draw a triangle (method 1): 

df1 %>%
  ggplot(aes(x = x_value, y = y_value)) + 
  geom_polygon(fill = "green") + 
  geom_path(color = "red", size = 2)

# Method 2: 

ggplot() + 
  geom_polygon(data = df1, aes(x = x_value, y = y_value), fill = "green") + 
  geom_path(data = df1, aes(x = x_value, y = y_value), color = "red", size = 2)

Draw a Basic Map

Để vẽ bản đồ của Việt Nam ở cấp quốc gia chẳng hạn, trước hết chúng ta cần lấy dữ liệu không gian (Spatial Data) tương ứng. Việc này có thể được thực hiện dễ dàng bằng hàm getData() như sau:

# Get geospatial data (nation level) for Viet Nam: 

geoData_vn_nation <- raster::getData("GADM", country = "Vietnam", level = 0)

# Convert to data frame: 

vn_nation_df <- geoData_vn_nation %>% fortify(region = "NAME_0")

Với Spatial Data đã có chúng ta có thể vẽ bản đồ cho Việt Nam bằng hàm geom_polygon() hoặc geom_path() theo cú pháp sau:

# Use geom_polygon(): 

vn_nation_df %>% 
  ggplot(aes(x = long, y = lat, group = group)) + 
  geom_polygon()

# Use geom_path(): 

vn_nation_df %>% 
  ggplot(aes(x = long, y = lat, group = group)) + 
  geom_path()

Hoặc kết hợp cả hai hàm này:

vn_nation_df %>% 
  ggplot(aes(x = long, y = lat, group = group)) + 
  geom_polygon(fill = "green") + 
  geom_path(color = "red")

Giả sử muốn thể hiện vị trí của Hà Nội chẳng hạn, trước hết lấy Spatial Data cho Hà Nội:

# Spatial data for Hanoi from https://www.latlong.net/place/hanoi-vietnam-14947.html

geo_hanoi <- data.frame(long = 105.804817, lat = 21.028511, group = "Vietnam.1")

Rồi hiển thị vị trí của Hà Nội bằng một điểm màu cam:

# Mapping Vietnam + Hanoi: 

vn_nation_df %>% 
  ggplot(aes(x = long, y = lat, group = group)) + 
  geom_path() + 
  geom_point(data = geo_hanoi, size = 3, color = "orange")

Còn để vẽ bản đồ Việt Nam ở cấp tỉnh, trước hết lấy Spatial Data ở cấp tương ứng:

# Get geospatial data (nation level) for Viet Nam: 

geoData_vn_province <- raster::getData("GADM", country = "Vietnam", level = 1)

# Convert to data frame: 

vn_province_df <- geoData_vn_province %>% fortify(region = "NAME_1")

Rồi vẽ bản đồ cho Việt Nam:

# Map of Vietnam, province level: 
vn_province_df %>% 
  ggplot(aes(x = long, y = lat, group = group, fill = id)) + 
  geom_polygon()

# Alternative: 

vn_province_df %>% 
  ggplot(aes(x = long, y = lat, group = group, fill = id)) + 
  geom_polygon(show.legend = FALSE)

Detailed Mapping

Trong nhiều trường hợp chúng ta có thể chỉ cần bản đồ cho, ví dụ, một số quận nội thành của Hà Nội. Trước hết lấy Spatial Data đến cấp huyện (vì đơn vị hành chính cấp Quận tương đương với Huyện) rồi chỉ lấy data tương ứng cho Hà Nội mà thôi:

library(stringi) # For processing text data. 

# Get geo-spatial data by district level for all provinces: 

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

# Only select Hanoi: 

hanoi <- vietnam_dis[vietnam_dis$NAME_1 == "Hà Nội", ]
# Convert to data frame and latin text for id column: 

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

# Select some districts: 

some_districts <- c("Hoan Kiem", "Hai Ba Trung", "Ba Dinh", "Tay Ho", 
                    "Hoang Mai", "Thanh Xuan", "Nam Tu Liem", "Bac Tu Liem",
                    "Cau Giay", "Dong Da", "Long Bien")

# Recode label: 

hanoi_df %>% mutate(id = case_when(!id %in% some_districts ~ "Others", TRUE ~ id)) -> hanoi_df_relabeled

# Only select inner districts: 

hanoi_df_relabeled %>% filter(id != "Others") -> geo_inner_hanoi

Với dữ liệu đã chuẩn bị chúng ta có thể vẽ bản đồ các quận nội thành (Version 1):

# Draw Map 1 (Only for inner districts): 

geo_inner_hanoi %>% 
  ggplot(aes(x = long, y = lat, group = group, fill = id)) +  
  geom_polygon(color = "white") + 
  labs(title = "Hanoi Map (Version 1) - Inner Districts", 
       x = "Longitude", 
       y = "Latitude", 
       caption = "Data Source: http://diva-gis.org/gdata")

Hoặc muốn nhấn mạnh các quận nội thành trong tổng thể toàn thành phố (Version 2):

hanoi_df_relabeled %>% filter(id == "Others") -> geo_sub_hanoi

# Version 2: 

geo_sub_hanoi %>% 
  ggplot(aes(x = long, y = lat, group = group), color = "gray69") + 
  geom_path() + 
  geom_polygon(data = geo_inner_hanoi, aes(fill = id), color = "black") + 
  labs(title = "Hanoi Map (Version 2)", 
       caption = "Data Source: http://diva-gis.org/gdata")

Draw a Sophisticated Map

Một kiểu mapping thể hiện nhiều thông tin hơn là Choropleth Map. Một ví dụ là thể hiện tỉ lệ nhiễm Covid 19.

Trước hết mô phỏng tỉ lệ nhiễm Covid 19 theo tỉnh:

#-----------------------------------
# Stage 1: Simulate Covid 19 data
#-----------------------------------

# All provinces: 

vn_provinces <- vn_province_df$id %>% unique()

# Number of provinces: 

n_prov <- length(vn_provinces)

# Covid data: 

covid_data <- runif(n = n_prov, min = 0, max = 0.9)

# Create a data frame: 

covid_df <- data.frame(id = vn_provinces, covid_rate = covid_data)

#-----------------------------------
# Stage 2: Join the two data sets
#-----------------------------------

# Join data sets: 

full_join(vn_province_df, covid_df, by = "id") -> vn_covid_geo_df

Với dữ liệu đã chuẩn bị như trên chúng ta có thể vẽ Choropleth Map:

vn_covid_geo_df %>% 
  ggplot(aes(x = long, y = lat, group = group, fill = covid_rate)) + 
  geom_polygon()

Hoặc một Choropleth Map có đường biên giới giữa các tỉnh:

vn_covid_geo_df %>% 
  ggplot(aes(x = long, y = lat, group = group, fill = covid_rate)) + 
  geom_polygon() + 
  geom_path()

A Real Case 1: Mapping PCI 2018 Indicator

Chúng ta sẽ tái lập lại bản đồ chỉ số PCI 2018 được in tại trang 26 của báo cáo Chỉ số năng lực cạnh tranh cấp tỉnh 2018. Trước hết download dữ liệu về chỉ số PCI 2018 từ đây rồi load dữ liệu này:

# Load PCI data: 

readxl::read_excel("_data_pci_file_general_file_1588582407-Indicator_PCI2018_VN.xlsx", sheet = 1) -> pci2018

# Select some column: 

pci2018 %>% select(1, 2, 3) -> pci2018_mini

# Remove some observations: 

pci2018_mini %>% slice(1:63) -> pci2018_63

# Rename for all columns: 

names(pci2018_63) <- c("id", "rank", "pci_score")

# Extract 63 provinces from PCI data: 

pci2018_63$id -> vn_provinces_pci

Trước khi join các bộ dữ liệu chúng ta phải kiểm tra sự đồng nhất của tên tỉnh. Dễ dàng thấy rằng bộ dữ liệu PCI có tỉnh là “BRVT” nhưng bộ dữ liệu Spatial thì lại là “Bà Rịa - Vũng Tàu”. Do vậy để join các dữ liệu lại cần phải đồng nhất tên tỉnh.

Vấn đề thứ hai liên quan đến encoding. Tỉnh “Yên Bái” từ bộ dữ liệu Spatial Data:

vn_provinces[63]
## [1] "Yên Bái"

Còn từ bộ PCI:

vn_provinces_pci[63]
## [1] "Yên Bái"

Nhìn bằng mắt thì cả hai đều là “Yên Bái” nhưng R sẽ không xem chúng là đồng nhất:

vn_provinces[63] == vn_provinces_pci[63]
## [1] FALSE

Một trong những giải pháp đơn giản là convert tên tỉnh ở hai bộ dữ liệu về kí tự Latin bằng hàm stri_trans_general() của thư viện stringi:

library(stringi) # For text processing. 

# Convert to Latin text: 

vn_province_df %>% mutate(id_latin = stri_trans_general(id, "Latin-ASCII")) -> vn_province_df

pci2018_63 %>% mutate(id_latin = stri_trans_general(id, "Latin-ASCII")) -> pci2018_63

Kiểm tra sự khác biệt về tên tỉnh bằng hàm setdiff() như sau:

setdiff(vn_province_df$id_latin, pci2018_63$id_latin)
## [1] "Ba Ria - Vung Tau" "Ho Chi Minh"       "Thua Thien Hue"

Từ phân tích này phương án xử lí là giữ nguyên tên tỉnh của vn_province_df và rename tên cho 3 tỉnh của bộ pci2018_63 bằng sử dụng hàm case_when() như sau:

pci2018_63 %>% 
  mutate(id_new = case_when(id_latin == "BRVT" ~ "Ba Ria - Vung Tau", 
                            id_latin == "TP.HCM" ~ "Ho Chi Minh", 
                            id_latin == "TT-Hue" ~ "Thua Thien Hue", 
                            TRUE ~ id_latin)) -> pci2018_63

# Remove id, id_latin column: 

pci2018_63 %>% select(-id, -id_latin) -> pci2018_63_homo

Lúc này chúng ta có thể join hai bộ dữ liệu một cách an toàn:

full_join(vn_province_df, pci2018_63_homo, by = c("id_latin" = "id_new")) -> df_for_mapping

Với dữ liệu “sạch” này chúng ta có thể tạo một Choropleth Map như sau:

df_for_mapping %>% 
  ggplot(aes(x = long, y = lat, group = group, fill = pci_score)) + 
  geom_polygon()

Làm đẹp cho PCI 2018 Map:

library(viridis) # For using viridis color palette. Ref: https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
library(showtext) # For using new font. 

# Select font for our map: 

my_font <- "Roboto Condensed"

font_add_google(name = my_font, family = my_font)

showtext_auto()

# Draw map: 

df_for_mapping %>% 
  ggplot(aes(x = long, y = lat, group = group, fill = pci_score)) + 
  geom_polygon() + 
  geom_path(color = "white", size = 0.1) + 
  coord_map(projection = "albers", lat0 = 30, lat1 = 40) + 
  labs(title = "Vietnam PCI Index 2018",
       subtitle = "Vietnam's Paracel and Spratly Islands\nare not shown in this map.",
       caption = "Data Source: http://pci2018.pcivietnam.vn") + 
  theme(axis.line = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid = element_blank(),
        plot.background = element_rect(fill = "white", color = NA),
        panel.background = element_rect(fill = "white", color = NA),
        legend.background = element_rect(fill = "white", color = NA),
        panel.border = element_blank()) +
  theme(plot.title = element_text(family = my_font, color = "grey20", size = 23)) + 
  theme(plot.subtitle = element_text(family = my_font, size = 13, color = "gray30")) + 
  theme(plot.caption = element_text(family = my_font, size = 10, colour = "grey40", face = "italic")) + 
  theme(legend.text = element_text(family = my_font, color = "grey40", size = 10)) + 
  theme(legend.title = element_text(family = my_font, color = "grey20", size = 10.5)) + 
  theme(plot.margin = unit(c(0.5, 0.3, 0.3, 0.3), "cm")) + 
  # theme(legend.position = "top") + 
  theme(legend.position = c(0.9, 0.7)) +
  scale_fill_viridis(direction = -1, 
                     option = "D", 
                     name = "CPI Index", 
                     guide = guide_colourbar(direction = "horizontal",
                                             barheight = unit(2, units = "mm"),
                                             barwidth = unit(35, units = "mm"),
                                             title.hjust = 0.5,
                                             label.hjust = 0.5, 
                                             title.position = "top"))

Mapping from sf Objects

Spatial Data mà chúng ta vừa sử dụng ở trên chưa có dữ liệu của Hoàng Sa và Trường Sa. Rất may là OpenDevelopmentMeKong cung cấp dữ liệu cho hai quần đảo này thuộc Việt Nam ở mọi cấp độ. Trước hết download dữ liệu, ví dụ, cho cấp huyện của Việt Nam từ đây sau đó sử dụng hàm sf::read_sf() để load bộ dữ liệu này:

# Load data: 

sf::read_sf("diaphanhuyen.geojson") -> geo_district_level

Dữ liệu load vào R là một sf object. Chúng ta có thể vẽ bản đồ các huyện của Việt Nam bằng sử dụng hàm ggplot2::geom_sf() như sau:

ggplot() + 
  geom_sf(data = geo_district_level) + 
  labs(title = "Mapping from sf Object")

Choropleth Map dân số theo cấp huyện:

geo_district_level %>% 
  ggplot(aes(fill = Dan_So)) + 
  geom_sf()

Hoặc vẽ một bản đồ hành chính cấp tỉnh với các huyện của nó:

geo_district_level %>% 
  ggplot(aes(fill = Ten_Tinh)) + 
  geom_sf(show.legend = FALSE)

Dữ liệu (ở dạng sf object) cũng có thể được load trực tiếp. Chẳng hạn dưới đây là vị trí (cũng như các thông tin khác) về tất cả các nhà máy thủy điện tại Việt Nam tính đến năm 2020:

my_link_data <- "https://data.vietnam.opendevelopmentmekong.net/geoserver/ODVietnam/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=ODVietnam%3Ahydropower_2020&outputFormat=application%2Fjson"

sf::read_sf(my_link_data) -> hydro_geoData

A Real Case 2: Mapping Hydroelectric Power Plants

Với những dữ liệu đã có chúng ta có thể hiển thị vị trí các nhà máy thủy điện tại Việt Nam tính đến năm 2020 bằng các điểm màu đỏ:

# Make a draft of map: 

ggplot() + 
  geom_path(data = vn_province_df, aes(x = long, y = lat, group = group)) + 
  geom_point(data = hydro_geoData, aes(x = X, y = Y), color = "red") + 
  coord_map(projection = "albers", lat0 = 30, lat1 = 40)

Chúng ta có thể lồng Bubble Chart vào bản đồ để thể hiện công suất của các nhà máy:

ggplot() + 
  geom_path(data = vn_province_df, aes(x = long, y = lat, group = group)) + 
  geom_point(data = hydro_geoData, aes(x = X, y = Y, size = Wattage_PL), color = "red", alpha = 0.3) + 
  coord_map(projection = "albers", lat0 = 30, lat1 = 40) + 
  scale_size(range = c(1, 10))

Hoàn thiện bản đồ này:

ggplot() + 
  geom_path(data = vn_province_df, aes(x = long, y = lat, group = group), color = "grey50") + 
  geom_point(data = hydro_geoData, aes(x = X, y = Y, size = Wattage_PL), color = "red", alpha = 0.3) + 
  scale_size(range = c(1, 10)) + 
  # coord_map(projection = "albers", lat0 = 30, lat1 = 40)  + 
  theme_minimal() + 
  labs(title = "Hydroelectric Stations by Power Capacity",
       subtitle = "This map shows hydroelectric power stations that are generating power using\nthe conventional dammed method. This list includes power stations that are\nlarger than 4 MW in maximum net capacity.",
       caption = "Data Source: https://data.opendevelopmentmekong.net/dataset/") + 
  theme(axis.line = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid = element_blank(),
        plot.background = element_rect(fill = "seashell", color = NA),
        panel.background = element_rect(fill = "seashell", color = NA),
        legend.background = element_rect(fill = "seashell", color = "grey50"),
        strip.background = element_rect(fill = "seashell"), 
        panel.border = element_blank()) +
  theme(plot.title = element_text(family = my_font, color = "grey20", size = 18)) + 
  theme(plot.subtitle = element_text(family = my_font, size = 11, color = "gray30")) + 
  theme(plot.caption = element_text(family = my_font, size = 10, colour = "grey40", face = "italic")) + 
  theme(legend.text = element_text(family = my_font, color = "grey40", size = 10)) + 
  theme(legend.title = element_text(family = my_font, color = "grey20", size = 11)) + 
  theme(plot.margin = unit(c(0.5, 0.3, 0.3, 0.3), "cm")) + 
  theme(legend.position = c(0.3, 0.45)) + 
  guides(size = guide_legend(title = "Net Capacity"))

Summary

Toàn bộ khóa học định hướng thực hành này được thiết kế với thời lượng học là 8h với mục tiêu khiêm tốn là vẽ những kiểu bản đồ thông dụng nhất thường gặp trong các báo cáo và nghiên cứu liên quan đến các vấn đề kinh tế - xã hội chỉ bằng sử dụng ba hàm của ggplot2 là geom_polygon(), geom_path() và geom_sf().

Không cần phải sử dụng những tools chuyên biệt, chỉ sử dụng R chúng ta cũng có thể tạo ra những bản đồ có chất lượng đủ tốt để đưa vào các báo cáo cũng như nghiên cứu. R codes dưới đây là một ví dụ:

library(osmdata) # For getting Open Street Map data. 

# Area that we want to collect spatial data: 

city_location <- "Ho Chi Minh city Vietnam"

# OSM objects that we want to get spatial: 

all_street_types <- c("motorway", "primary", "secondary", "tertiary")

secondary_streets <- c("residential", "living_street", "unclassified", "service", "footway")

# Collect spatial data for OSM objects selected: 

streets <- getbb(city_location) %>%
  opq() %>%
  add_osm_feature(key = "highway", value = all_street_types) %>%
  osmdata_sf()

small_streets <- getbb(city_location) %>%
  opq() %>%
  add_osm_feature(key = "highway", value = secondary_streets) %>%
  osmdata_sf()

river <- getbb(city_location)%>%
  opq() %>%
  add_osm_feature(key = "waterway", value = "river") %>%
  osmdata_sf()

# Long and Lat limits: 

geo_spatial_limits <- getbb(city_location)

geo_spatial_limits

# Select colours and font for mapping: 

color1 <- "cyan"
color2 <- "#ffbe7f"
back_color <- "#1d1330"


# Make a draft version of streetmap for Ho Chi Minh city: 

ggplot() +
  geom_sf(data = streets$osm_lines, inherit.aes = FALSE, color = color1, size = 0.4, alpha = 0.8) +
  geom_sf(data = small_streets$osm_lines, inherit.aes = FALSE, color = color2, size = 0.2, alpha = 0.6) +
  geom_sf(data = river$osm_lines, inherit.aes = FALSE, color = "blue", size = 1) +
  coord_sf(xlim = c(106.60, 106.75), ylim = c(10.75, 10.86), expand = FALSE) + 
  theme(plot.background = element_rect(fill = back_color, colour = NA),
        panel.background = element_rect(fill = back_color, colour = NA), 
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        panel.grid = element_blank()) + 
  theme(plot.title = element_text(family = my_font, size = 22, color = "white")) + 
  theme(plot.margin = unit(rep(0.5, 4), "cm")) + 
  labs(title = "Ho Chi Minh City from Open Street Map Database")

Để render ra bản đồ này cần các dữ liệu về đường lớn (xanh nhạt), đường nhỏ (vàng) và đường sông (xanh đậm) với tổng dung lượng là 3.6 GB nên mất một thời gian khá lâu (chừng 20 phút hoặc hơn tùy cấu hình máy tính và chất lượng mạng).

Do thời lượng có hạn nên khóa học không thể cover mọi khía cạnh của Mapping cũng như Spatial Data Visualization. Mục References dưới đây liệt kê những tài liệu tham khảo đáng tin cho những ai muốn nâng cao kĩ năng mapping và Spatial Data Visualization