Nhu cầu vẽ đồ thị bản đồ

Trong nghiên cứu y học nhu cầu vẽ đồ thị bản đồ thường có trong những nghiên cứu dịch tễ, biểu diễn tần suất bịnh liên quan đến địa phương.

Biểu đồ Barplot cho ta ý niệm về tần suất (frequency), phần trăm, tuy nhiên thông tin về các vùng địa lí không được thể hiện. Vì thế không gợi ý cho ta mối liên quan về tần suất bệnh và địa phương.

Trong biểu đồ Barplot bên trái, các địa phương Hà Nội, Hòa Bình, Thái Nguyên và Bắc Giang nằm xa nhau. Vì thế không gợi ý cho ta mối liên hệ có thể có giữa các địa phương này với vấn đề dịch tễ đang được nghiên cứu.

Trong bản đồ bên phải, trực quan đặt cho ta câu hỏi vì sao các địa phương gần nhau ở trên đều có vấn đề dịch tễ này.

Điều đó nói lên tầm quan trọng của vẽ bản đồ trong nghiên cứu y học.

Chuẩn bị

Để vẽ được bản đồ dịch tễ ta cần chuẩn bị:

Vẽ với tmaps

Đọc dataset vào R

library("tmap")
library("tmaptools")
library("sf")
library("leaflet")
library("readxl")
data <- read_xlsx("D:/Vnprovinces.xlsx")
head(data)
## # A tibble: 6 x 3
##      id Provinces         Events
##   <dbl> <chr>              <dbl>
## 1  1.00 Ha Noi City|Hanoi   7.00
## 2  2.00 Thai Nguyen         3.00
## 3  3.00 Hoa Binh            1.00
## 4  4.00 Nghe An             4.00
## 5  5.00 Thai Binh           2.00
## 6  6.00 Yen Bai             3.00

Đọc bản đồ vào R

vnshapefile <- "F:/R/Data/Vietnam_shapefile_maps/National Only/Provinces.shp"
vngeo <- read_shape(file = vnshapefile, as.sf = TRUE)

Kiểm tra bản đồ đã đọc được chính xác chưa. Nếu đọc đúng, sẽ vẽ được một bản đồ trống

qtm(vngeo)

head(vngeo)
## Simple feature collection with 6 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 103.8718 ymin: 10.13411 xmax: 108.0788 ymax: 22.73976
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##   ISOCODE_2                  NAME_1 ID_2      NAME_2       VARNAME_2
## 0     VN-07                Ðông B?c   34 Tuyên Quang     Tuyen Quang
## 1     VN-53                Ðông B?c   26     B?c K?n Bac Kan|Bac Can
## 2     VN-45 d?ng b?ng sông C?u Long    1   Ð?ng Tháp       Dong Thap
## 3     VN-13                Ðông B?c   32  Qu?ng Ninh      Quang Ninh
## 4     VN-22            B?c Trung B?   45     Ngh? An         Nghe An
## 5     VN-54                Ðông B?c   25   B?c Giang       Bac Giang
##   TYPE_2 ENGTYPE_2                       geometry
## 0   T?nh  Province MULTIPOLYGON (((105.1785 22...
## 1   T?nh  Province MULTIPOLYGON (((105.7406 22...
## 2   T?nh  Province MULTIPOLYGON (((105.4328 10...
## 3   T?nh  Province MULTIPOLYGON (((107.151 20....
## 4   T?nh  Province MULTIPOLYGON (((105.7822 18...
## 5   T?nh  Province MULTIPOLYGON (((106.1623 21...

Đến đây ta đã đưa được bản đồ vào R. Việc kế tiếp là kết nối dữ liệu với bản đồ.

Nhìn vào kết quả của head(vngeo) ta thấy VARNAME_2 chính là variable để kết nối với variable Provinces trong dataset. Chúng ta kết nối dữ liệu với bản đồ như sau:

vnmap <- append_data(vngeo, data, key.shp = "VARNAME_2", key.data ="Provinces")

Một message xuất hiện:

“Under coverage: 1 out of 63 shape features did not get appended data. Run under_coverage() to get the corresponding feature id numbers and key values.”

Nó cho biết có 1 địa phương không xuất hiện trong dữ liệu, chúng ta phải tìm ra địa phương này để check lại dataset bằng cách chạy under_coverage()

under_coverage()
## $result
## [1] "Under coverage: 1 out of 63 shape features did not get appended data."
## 
## $call
## [1] "append_data(shp = vngeo, data = data, key.shp = \"VARNAME_2\", "
## [2] "    key.data = \"Provinces\")"                                  
## 
## $id
## [1] 49
## 
## $value
## [1] Quang Tri
## 63 Levels: An Giang Ba Ria - VTau|Ba Ria-Vung Tau ... Yen Bai

Trong danh sách của file bản đổ (VARNAME_2) có Quang Tri, nhưng trong dataset không có Quang Tri. Phải kiểm tra lại dataset của bạn.

Khi bạn chắc về error trên, chúng ta vẽ bản đồ:

qtm(vnmap, "Events")

Chúng ta đã có một bản đồ bằng tmaps package đầy đủ các tỉnh. Trong đó màu sắc thể hiện frequency của Events từ 0 đến 7. Trong đó dữ liệu của tỉnh Quảng Trị bị missing.

Nếu muốn tỉnh Quảng Trị không bị mất trên bản đồ, bạn mở dataset nhập vào dòng Quang Tri với Events = 0. Nhưng ở đây tôi không nhập Quang tri để muốn minh họa cho các bạn thấy khi một địa phương bị mất dữ liệu thì trên bản đồ thể hiện ra sao.

vẽ với ggplots

ggplot2 cho ta nhiều tính năng vẽ bản đồ sinh động hơn về màu sắc, scale về sắc độ, text…

Bây giờ chúng ta hãy thử vẽ với ggplot2

library(rgeos)
library(maptools)
library(rgdal)
library(gpclib) 
library(ggplot2)
vnshapefile <- "F:/R/Data/Vietnam_shapefile_maps/National Only/Provinces.shp"
vngeo2 <- readOGR(vnshapefile)
## OGR data source with driver: ESRI Shapefile 
## Source: "F:/R/Data/Vietnam_shapefile_maps/National Only/Provinces.shp", layer: "Provinces"
## with 63 features
## It has 7 fields
vngeo2 <- fortify(vngeo2, region = "VARNAME_2")

Dữ liệu đã được đọc vào R với packages khác tương hợp với ggplot2

Vẽ với ggplot2

ggplot() + geom_map(data = data, aes(map_id = Provinces, fill = Events), 
                    map = vngeo2) + expand_limits(x = vngeo2$long, y = vngeo2$lat)

Có những tồn tại làm cho bạn chưa thật hài lòng với bản đồ này. Chúng ta phải cải thiện màu sắc, sắc độ của nó thôi. Tỉnh Quảng Trị bị mất trong dataset nên không được thể hiện trên bản đồ.

library(scales)
ggplot() + geom_map(data = data, aes(map_id = Provinces, fill = Events), 
  map = vngeo2) + expand_limits(x = vngeo2$long, y = vngeo2$lat) + 
  scale_fill_gradient2(low = muted("red"), 
  mid = "white", midpoint = 2, high = muted("blue"), limits = c(0, 7))

Màu xanh đậm thể hiện frequency cao nhất, màu nâu ngược lại, màu trắng cho freqnency trung gian. Chưa logic lắm.

Nếu bạn không thích logic này, thay vào đó logic màu sắc khác, bạn viết code như sau:

# color changing
ggplot() + geom_map(data = data, aes(map_id = Provinces, fill = Events), 
                             map = vngeo2) + expand_limits(x = vngeo2$long, y = vngeo2$lat) + 
  scale_fill_gradientn(colours = c("white","green" , "yellow", "brown", "darkred"))

Màu nâu đậm đại diện cho vùng có frequency cao nhất, màu trắng chỉ vùng có freqnency thấp nhất. Chấp nhận được.

Còn nữa, nếu bạn muốn một màu mà thôi thì thử làm như sau:

# one color
ggplot() + geom_map(data = data, aes(map_id = Provinces, fill = Events), 
  map = vngeo2) + expand_limits(x = vngeo2$long, y = vngeo2$lat) + 
  scale_fill_gradient(low="white", high="blue")

Nếu bạn hài lòng với cách thể hiện sắc độ này thì phần tiếp theo là đưa text vào bản đồ cho dễ hiểu hơn.

library(plyr)
provincecenters <- ddply(vngeo2, .(id), summarize, clat = mean(lat), clong = mean(long))
ggplot() + geom_map(data = data, aes(map_id = Provinces, fill = Events), 
  map = vngeo2) + expand_limits(x = vngeo2$long, y = vngeo2$lat) + 
  scale_fill_gradient(low="white", high="blue")    + 
  geom_text(data = provincecenters, aes(x = clong, y = clat, label = id), size = 2.7)+
  ggtitle("Events By Provices")

Nếu bạn muốn có 1 bản đồ không có trục x, y, không có kinh tuyến, vĩ tuyến thì bạn làm như sau:

library(ggthemes)
ggplot() + geom_map(data = data, aes(map_id = Provinces, fill = Events), 
  map = vngeo2) + expand_limits(x = vngeo2$long, y = vngeo2$lat) + 
  theme_few()+
  theme(legend.position = "right",
        axis.ticks = element_blank(), 
        axis.title = element_blank(), 
        axis.text =  element_blank()) +
  scale_fill_gradient(low="white", high="blue") +
  guides(fill = guide_colorbar(barwidth = 0.6, barheight = 10)) + 
  geom_text(data = provincecenters, aes(x = clong, y = clat, label = id), size = 2.9)+
  ggtitle("Events By Provices")

Bạn có thể thay đổi các chỉ số để có được màu sắc, font size như mong muốn.

Chúc các bạn có những trãi nghiệm thú vị.