流行病學(xué)中發(fā)病情況的地區(qū)分布圖是很常見的一類做圖，R中的ggplot2是一個(gè)非常優(yōu)雅的繪圖工具，通過原始的地理測(cè)繪數(shù)據(jù)配合ggplot2繪制地圖是比較硬核的一種方法，國家測(cè)繪中心早年免費(fèi)提供地理測(cè)繪數(shù)據(jù)，最近的免費(fèi)數(shù)據(jù)是2004年的測(cè)繪數(shù)據(jù)( https://uploads.cosx.org/2009/07/chinaprovinceborderdata_tar_gz.zip )[1]，2014年新版的測(cè)繪數(shù)據(jù)則是收費(fèi)的。兩者最直觀的區(qū)別在于重慶。如果僅是練手可以使用2004年的測(cè)繪數(shù)據(jù)。
流行病學(xué)中比較常用的是氣泡圖，這個(gè)在網(wǎng)上已經(jīng)有很多教程。但用這種方法繪制地圖時(shí)直接按省份進(jìn)行染色的教程竟然一篇都沒有，摸索了一下，繪制出了相對(duì)完美的分布圖（發(fā)病人數(shù)使用的是1月25日的數(shù)據(jù)）：
plot_zoom.png
這個(gè)地圖的制圖難點(diǎn)在于原始的測(cè)繪數(shù)據(jù)雖然是34個(gè)省市自治區(qū)的（原數(shù)據(jù)實(shí)際沒有標(biāo)出澳門），但因?yàn)橛械氖》菝娣e太大，所以一共花了924個(gè)數(shù)據(jù)用于描述這34個(gè)省市自治區(qū)，這被稱為描述層。為了畫出描述層，一共有9萬多條幾何層的數(shù)據(jù)，這些幾何層的數(shù)據(jù)是一些坐標(biāo)點(diǎn)，相互連接起來之后就畫出了各省份的邊界，而同一組描述層內(nèi)的幾何層數(shù)據(jù)會(huì)拼接成一個(gè)省份。
如果是普通的34個(gè)數(shù)據(jù)，那是最容易繪圖的了，直接生成34個(gè)顏色就可以，但原始的數(shù)據(jù)是個(gè)S4類數(shù)據(jù)，里面的data子數(shù)據(jù)包含了924個(gè)描述層數(shù)據(jù)，通過fortify函數(shù)解析這個(gè)S4數(shù)據(jù)可以獲得9萬多條幾何層的繪圖數(shù)據(jù)。R自帶的plot函數(shù)是直接通過描述層調(diào)用幾何層，因此需要924個(gè)顏色去填充，正是通過一篇plot函數(shù)繪圖的文章讓我想到了ggplot2的做圖方法[2]。ggplot2中是直接按9萬多條原始繪圖數(shù)據(jù)去制圖的，所以需要9萬多個(gè)顏色，我直接套用了原文中生成顏色的函數(shù)，做出來的圖已經(jīng)可以用于展示。代碼中已添加了詳細(xì)注釋：
# 加載包
library(ggplot2) # 繪圖
library(rgdal) # 讀取地圖數(shù)據(jù)
library(plyr) # 數(shù)據(jù)處理
# 全局參數(shù)
options(stringsAsFactors = FALSE) # 這個(gè)改了省心，r喜歡按因子讀取字符
windowsFonts("Times New Roman" = windowsFont("Times New Roman")) # 設(shè)置字體
# 讀取文件，設(shè)置工作路徑
setwd("E:/Store/code/03.R/03.ggplot2/06.map/test") # 根據(jù)實(shí)際情況改自己的工作路徑
vData <- readOGR("E:/Store/code/03.R/code/06.maps/map/2014/maps/bou2/bou2_4p.shp") # 按地圖文件的實(shí)際位置修改
# 數(shù)據(jù)處理
vTemp <- vData@data; # S4類數(shù)據(jù),用@取子集,得到的為描述層
vDescribeData <- data.frame(vTemp, id = seq(0:924) - 1) # 為描述層添加id，方便后續(xù)內(nèi)聯(lián)數(shù)據(jù)
vFormatData <- fortify(vData) # 得倒的為幾何層數(shù)據(jù)
vMapData<-join(vFormatData, vDescribeData, type = "full") # 合并描述層和幾何層數(shù)據(jù)（普通做圖這一步并不是必須的，通常只需要幾何層就可以做出地圖，合并起來更直觀一些）
# 讀取各省份的發(fā)病人數(shù)
vNum <- read.csv("num.csv", stringsAsFactors = FALSE)
vProvName <- vNum$vProvName # 省份名稱
vNum <- vNum$vNum # 發(fā)病人數(shù)
vProvCol = vNum # 將發(fā)病人數(shù)傳給顏色，用于表示顏色的深淺
# 顏色函數(shù)，從他人博客套用的函數(shù)，用處是給每一個(gè)描述層生成一個(gè)顏色（描述層有924條數(shù)據(jù)）
fGetColor <- function(mapdata, vProvName, vProvCol, othercol){
f = function(x, y) ifelse(x %in% y, which(y == x), 0)
colIndex = sapply(mapdata@data$NAME, f, vProvName)
fg = c(othercol, vProvCol)[colIndex + 1]
return(fg)
}
# 使用自定義的顏色函數(shù)，為描述層生成顏色
vCol <- fGetColor(vData, vProvName, vProvCol, 0)
# 為描述層顏色添加id
vRealCol <- data.frame(vCol, id = seq(0:924) - 1)
# 通過id將描述層的顏色數(shù)據(jù)和之前的幾何層地圖數(shù)據(jù)內(nèi)聯(lián)，簡單的說就是將924條顏色數(shù)據(jù)根據(jù)id映射到9萬多條幾何層數(shù)據(jù)
vMapData <- join(vMapData, vRealCol, type = "full")
# ggplot2常規(guī)做圖，通過scale_fill_continuous函數(shù)產(chǎn)生漸變色
ggplot(vMapData, aes(x = long, y = lat, group = group, fill = as.numeric(vMapData$vCol)))+
geom_polygon() +
scale_fill_continuous(low = "white", high = "red") +
geom_path(colour = "black") +
ggtitle("新型肺炎確診病例分布示意圖") +
theme(plot.title = element_text(lineheight = 0.8, face = "bold", size = 20)) +
theme(axis.text = element_text(size = 15, colour = "black", family = "Times New Roman"),
axis.title = element_blank()) +
# specify tick marks
coord_map(xlim = c(73, 136), ylim = c(5, 54)) +
scale_y_continuous(breaks = seq(0, 60, 10)) +
scale_x_continuous(breaks = seq(70, 135, 10)) +
# remove guide
guides(fill = FALSE)
其中num.csv的格式如下：
num.png
下載地址為：http://www.1x1y.top/books/lib/exe/fetch.php?media=num.zip
vNum2列是真實(shí)發(fā)病人數(shù)，我用了漸變色會(huì)導(dǎo)致人數(shù)較低時(shí)偏白色，所以我對(duì)發(fā)病人數(shù)做了一點(diǎn)轉(zhuǎn)換，實(shí)際繪圖用的是轉(zhuǎn)換后的vNum列數(shù)據(jù)= =。將發(fā)病人數(shù)生成新的一列評(píng)級(jí)數(shù)據(jù)作為染色標(biāo)準(zhǔn)時(shí)做出的圖會(huì)更好。由于數(shù)據(jù)轉(zhuǎn)換以及添加地名比較簡單，我就沒有再接著做了。只要再配合一個(gè)爬蟲爬取發(fā)病人數(shù)，就可以實(shí)時(shí)繪制最新的疫情分布圖了。地圖文件參考文中給出的鏈接下載。
參考文獻(xiàn)：
[1] 邱怡軒, https://cosx.org/2009/07/drawing-china-map-using-r/.2009.
[2] [阿蠻的杜鵑], https://www.cnblogs.com/lizhilei-123/p/6734378.html.2017.
