R繪制中國地圖，并展示流行病學數據
本文作者：姜曉東，博士畢業于上海交通大學，目前任教于湖南師范大學醫學院，專業神經毒理學。
流行病學的數據講究“三間分布”，即人群分布、時間分布和空間分布。其中的“空間分布”最好是在地圖上展示，才比較清楚。R軟件集統計分析與高級繪圖于大成，是最適合做這項工作了。關于地圖的繪制過程，謝益輝、邱怡軒和陳麗云等人都早有文章講述，開R地圖中文教程之先河。由于目前指導畢業論文用到，因此研究了一下。本來因為網上教程很多，曾打消了寫些文字的計劃，但怡軒版主鼓勵說“教程者眾，整合者鮮”，所以才戰勝拖延癥，提起拙筆綜述整合一下，并對DIY統計GIS地圖提出了一點自己的想法。
1 地圖GIS數據的來源與R繪制軟件包
中國地圖GIS數據的官方數據可以在國家基礎地理信息中心的網站（??http://nfgis.nsdi.gov.cn??）里面可以免費下載。官方公開的數據包括：地圖數據，及居住地、交通、河流等輔助數據。今年6月開始，官方正組織開始制作新版數據。老數據暫時無法下載，讀者要自行百度搜索，本文以舊版數據為例。舊版地圖數據中部分地名和地市區劃已經過時，使用時需注意。
地圖數據有4個壓縮文件：bou1_4m.zip、bou2_4m.zip、bou3_4m.zip和bou4_4m.zip。bou代表邊界的意思，數字1~4代表國家、省、市、縣的4級行政劃分；4m代表比例是400萬分之一，這個比例的圖形是公開的。每個文件解壓縮后含有兩類文件：以字母p結尾的表示多邊形數據，用來繪制區域；以字母l結尾的文件是線形數據，用來繪制邊界。但是老版數據中，市級數據中缺少繪制區域的多邊形數據，讓市級分布圖的繪制稍麻煩一些，新版中也許會有改進。
用R繪制地圖比較簡單。比如畫一下全國范圍的區域，可以用如下代碼：
library(maptools)
mydat = readShapePoly("maps/bou1/bou1_4p.shp")
plot(mydat)
但是，可以看出這樣繪制的地圖的形狀有些扁平。這是因為，在繪圖的過程中，默認把經度和緯度作為普通數據，均勻平等對待，繪制在笛卡爾坐標系上造成的。其實，地球的球面圖形如何映射到平面圖上，在地理學上是有一系列不同的專業算法的。地圖不應該畫在普通的笛卡爾坐標系上，而是要畫在地理學專業的坐標系上。在這一點上，R的ggplot2包提供了專門的??coord_map()??函數。所以推薦R的ggplot2包來繪制地圖。
library(ggplot2)
mymap = ggplot(data = fortify(mydat)) +
geom_polygon(aes(x = long, y = lat, group = id), colour = "black",
fill = NA) +
theme_grey()
print(mymap + coord_map())
這次中國地圖的形狀與百度地圖一樣了。
ggplot2包的??coord_map??函數默認的映射類型是mercator。如果有其他需要，可以使用其他的映射類型來繪制地圖，如：
mymap + coord_map(projection = "azequidistant")
??coord_map??函數的映射類型及其含義可以通過下列代碼查詢幫助，一般我們用默認的就可以。
library(mapproj)
?mapproject
2 GIS地圖的數據結構及省市地圖的繪制
GIS地圖有很多種存儲格式，其中shapefile格式（.shp）可以通過R的maptools包打開。其他格式可以去R官網查詢相應的軟件包。
地圖數據基本可以分為點、線、面三種數據，在maptools包內分別有對應的函數來讀?。??readShapePoints???、??readShapeLines???和??readShapePoly???函數）。首先以面（poly）型數據介紹。先看代碼，通過??readShapePoly??函數讀入省級地圖：
library(maptools)
mydat = readShapePoly("maps/bou2/bou2_4p.shp")
此時，??mydat???中保存的是各個省/直轄市的多邊形面圖，數據類型是??SpatialPolygonsDataFrame??。我們可以：
length(mydat)
## [1] 925
names(mydat)
## [1] "AREA" "PERIMETER" "BOU2_4M_" "BOU2_4M_ID" "ADCODE93"
## [6] "ADCODE99" "NAME"
可以發現??mydat??中有925條記錄，每條記錄中含有面積（AREA）、周長（PERIMETER）、各種編號、中文名（NAME）等字段。其中中文名（NAME）字段是以GBK編碼的。
這個??SpatialPolygonsDataFrame???類型并不是真正的??data.frame???類型，而是一個sp包定義的類，只不過重載了 ??[]??? 和 ??$??? 運算符，使得一些行為上與??data.frame??相類似。
可以進一步統計一下，每個省/直轄市的多邊形數目。
table(iconv(mydat$NAME, from = "GBK"))
##
## 上海市 云南省 內蒙古自治區 北京市
## 12 1 1 1
## 臺灣省 吉林省 四川省 天津市
## 57 1 1 1
## 寧夏回族自治區 安徽省 山東省 山西省
## 1 1 86 1
## 廣東省 廣西壯族自治區 新疆維吾爾自治區 江蘇省
## 154 6 1 5
## 江西省 河北省 河南省 浙江省
## 1 9 1 179
## 海南省 湖北省 湖南省 甘肅省
## 79 1 1 1
## 福建省 西藏自治區 貴州省 遼寧省
## 168 1 2 94
## 重慶市 陜西省 青海省 香港特別行政區
## 1 1 1 53
## 黑龍江省
## 1
我的環境是UTF-8，所以需要??iconv??函數轉化一下才能正常顯示。
結果顯示多數省的地圖都是由一個多邊形構成，少數臨海省/直轄市由于有很多附屬島嶼，多邊形數目比較多。
利用與??data.frame???相似的 ??[]??? 和 ??$?? 運算符操作，我們可以迅速提取出一個省市的數據，比如上海及附屬崇明島：
Shanghai = mydat[mydat$ADCODE99 == 310000,]
plot(Shanghai)
?
?
其中ADCODE99是國家基礎地理信息中心定義的區域代碼，共有6位數字，由省、地市、縣各兩位代碼組成。
為了進一步在ggplot2包中繪圖，需要把??SpatialPolygonsDataFrame???數據類型轉化為真正的??data.frame???類型才可以。ggplot2包專門針對地理數據提供了特化版本的??fortify??函數來做這個工作：
head(fortify(Shanghai))
## long lat order hole piece group id
## 1 121.3 31.85 1 FALSE 1 208.1 208
## 2 121.3 31.85 2 FALSE 1 208.1 208
## 3 121.3 31.85 3 FALSE 1 208.1 208
## 4 121.3 31.85 4 FALSE 1 208.1 208
## 5 121.3 31.84 5 FALSE 1 208.1 208
## 6 121.4 31.83 6 FALSE 1 208.1 208
#p#
3 在地圖上展示流行病學數據
3.1 一地名對應一區域，長沙為例
首先把長沙所轄地區找到，這個可以根據ADCODE99編碼的前4位定位長沙，去查表就可以了。但是這個地名是99年的標準，新版正在制定過程中，隨時會變。我們權且以此為例。如果找不到表，可以通過代碼在命令行下手工查找：
mydat = readShapePoly("maps/bou4/BOUNT_poly.shp")
tmp = iconv(mydat$NAME99, from = "GBK")
grep("長沙", tmp, value = TRUE)
## [1] "長沙縣" "長沙市市轄區"
grep("長沙", tmp)
## [1] 2122 2183
mydat$ADCODE99[grep("長沙", tmp)]
## [1] 430121 430101
## 2368 Levels: 0 110100 110112 110113 110221 110224 110226 110227 ... 820000這樣我們就知道了長沙ADCODE99編碼的前4位是4301，其中43代表湖南省，01就是長沙市。接著就可以篩選出長沙的地圖數據：
Changsha = mydat[substr(as.character(mydat$ADCODE99), 1, 4) == "4301",]
mysh = fortify(Changsha, region = 'NAME99')
mysh = transform(mysh, id = iconv(id, from = 'GBK'), group = iconv(group, from = 'GBK'))
head(mysh, n = 2)
## long lat order hole piece group id
## 1 113.1 28.18 1 FALSE 1 長沙市市轄區.1 長沙市市轄區
## 2 113.1 28.18 2 FALSE 1 長沙市市轄區.1 長沙市市轄區
names(mysh)[1:2] = c("x","y") #這句是不得已而為之的黑魔法
接著我們給一串隨機數當成是流行病學數據，并用顏色填充到地圖上。
myepidat = data.frame(id = unique(sort(mysh$id)))
myepidat$rand = runif(length(myepidat$id))
myepidat
## id rand
## 1 寧鄉縣 0.98076
## 2 望城縣 0.32123
## 3 瀏陽市 0.66957
## 4 長沙縣 0.09655
## 5 長沙市市轄區 0.19437
csmap = ggplot(myepidat) +
geom_map(aes(map_id = id, fill = rand), color = "white", map = mysh) +
scale_fill_gradient(high = "darkgreen",low = "lightgreen") +
expand_limits(mysh) + coord_map()
print(csmap)
?
?
接下來的工作就是添加地名，sp包提供了??coordinates??函數，來計算地圖的中心坐標：
tmp = coordinates(Changsha)
print(tmp)
## [,1] [,2]
## 2121 113.2 28.32
## 2134 113.7 28.23
## 2136 112.8 28.29
## 2149 112.3 28.13
## 2182 113.0 28.17
tmp = as.data.frame(tmp)
tmp$names = iconv(Changsha$NAME99, from = 'GBK')
print(tmp)
## V1 V2 names
## 2121 113.2 28.32 長沙縣
## 2134 113.7 28.23 瀏陽市
## 2136 112.8 28.29 望城縣
## 2149 112.3 28.13 寧鄉縣
## 2182 113.0 28.17 長沙市市轄區
csmap + geom_text(aes(x = V1,y = V2,label = names), family = "GB1", data = tmp)
?
?
如果需要支持更多字體，可以配合使用showtext包。
3.2 內地省份的地市級圖的情況
如果國家基礎地理信息中心的GIS地圖數據的地市文件bou3_4m.zip中含有polygon文件，那么我們就可以根據上一節的內容繪制省內地市級分布圖了。官方恰恰缺少了這個文件，給繪圖造成了麻煩。解決方案有兩個：一個是另辟蹊徑，從非官方的??www.gadm.org??下載一份shp格式的中國地圖來繪制；另一個解決方案是從官方發布的縣級地圖入手，根據ADCODE99編碼適當合并，繪制省內地市分布圖，同時利用bou3_4m.zip僅存的邊界文件繪制邊界。
相信官方新版本的GIS地圖數據會包含舊版本所缺失的這份文件。目前還是建議暫時使用gadm的省級地圖。舊版官方地圖信息比較陳舊落后，比如湖南沒有標注出湘西州的規劃。
3.3 一地名對應多區域，上海為例
中國很多沿海省/直轄市有很多附屬島嶼，導致地名和區域（Polygon）存在一對多的情況。這種情況下，在??fortify??處理數據的時候一定要特別注意索引與多邊形一一對應，同時又要保持地名信息，黑魔法在代碼中：
# mydat = readShapePoly("maps/bou4/BOUNT_poly.shp")
Shanghai = mydat[substr(as.character(mydat$ADCODE99), 1, 2) == '31',]
mysh = fortify(Shanghai, region = 'NAME99')
mysh = transform(mysh, id = iconv(id, from = 'GBK'), group = iconv(group, from = 'GBK'))
head(mysh)
## long lat order hole piece group id
## 1 121.2 31.85 1 FALSE 1 崇明縣.1 崇明縣
## 2 121.3 31.85 2 FALSE 1 崇明縣.1 崇明縣
## 3 121.3 31.85 3 FALSE 1 崇明縣.1 崇明縣
## 4 121.3 31.85 4 FALSE 1 崇明縣.1 崇明縣
## 5 121.3 31.85 5 FALSE 1 崇明縣.1 崇明縣
## 6 121.3 31.84 6 FALSE 1 崇明縣.1 崇明縣
# 黑魔法在此
names(mysh)[c(1, 2, 6, 7)] = c("x", "y", "id", "code")
myepidat = data.frame(id = unique(sort(mysh$id)))
# 隨機數字替代數據
myepidat$rand = runif(length(myepidat$id))
# 官方地圖區劃比較落后過時，目前上海是16區1縣，神碼“市直轄5區”的稱呼已經過時。
myepidat
## id rand
## 1 上海市市轄區.1 0.21673
## 2 上海市市轄區.2 0.74173
## 3 上海市市轄區.3 0.02462
## 4 上海市市轄區.4 0.20619
## 5 上海市市轄區.5 0.89970
## 6 南匯縣.1 0.77084
## 7 嘉定區.1 0.21771
## 8 奉賢縣.1 0.91729
## 9 崇明縣.1 0.04879
## 10 崇明縣.2 0.02462
## 11 崇明縣.3 0.03397
## 12 崇明縣.4 0.72591
## 13 崇明縣.5 0.72059
## 14 崇明縣.6 0.43981
## 15 松江區.1 0.18296
## 16 金山區.1 0.78371
## 17 金山區.2 0.88552
## 18 閔行區.1 0.54186
## 19 青浦縣.1 0.12003
ggplot(myepidat) + geom_map(aes(map_id = id, fill = rand), map = mysh) +
expand_limits(mysh) + coord_map()?
?
3.4 其他問題
如果需要縣級以下的地圖GIS數據，比如街道、鄉村的地圖，國家地理信息中心并不提供。要么去民政部索取，要么自己繪制。
另外，提醒大家，流行病學數據并不是僅僅畫在地圖上就完事了。針對空間數據，R里面有很多空間數據的分析軟件包。推薦Roger S. Bivand的《Applied Spatial Data Analysis with R》，尤其是里面第11章“Disease Mapping”，對醫學背景同學很有益處。如果能找到一個地理資源環境學院的研究生一同討論的話就更好了。畢竟，它山之石可以攻玉，我們要承認自己的不足。
#p#
4 自己繪制簡單的GIS地圖
在制作流行病學統計地圖的過程中，對于很多區、街道、鄉村級別的地圖，無法獲得GIS數據。很多人的做法是到百度地圖上用繪圖軟件摹描出區域線圖，然后再把自己的數據計算成相應顏色，再手工填充顏色繪成統計地圖。這個過程枯燥繁瑣，而且數據映射成顏色的時候容易出錯。不如把你已經描好的線圖，制成shp格式的GIS數據地圖，分享給大家用。辛苦你一個，幸福千萬家。這個過程其實有專業的GIS軟件可以做，若你能找到專業人士，就直接“幸福千萬家”了。
如果地圖結構簡單，我們可以“土法”來做。先去NIH（美國國立衛生研究院）網站下載一個免費的圖像軟件ImageJ，用來采集地區邊界數據。然后再把采集好的數據在R軟件里面把像素坐標換算成地理坐標，在利用R軟件sp包和maptools的函數整合成??SpatialPolygonsDataFrame??，最后保存為shp格式的地圖文件。
我們以起點中文網小說《江山美人志》開篇所附地圖為例，繪制虛擬世界里面“中南郡”的GIS地圖。為了和實際問題類似，我在地圖中畫上了參考坐標線。
?
?
利用ImageJ“點”工具，同時按住Shift鍵一次批量多點采樣，再點擊分析菜的測量，最后保存結果。
ImageJ采集的點坐標是位圖像素相對坐標，為了能換算為地理經緯度坐標。我們先采集圖上參考坐標線上的經緯交點坐標，在R中建立換算關系：
mg_pos = data.frame(x = c(103,103,403,403), y = c(75,275,75,275))
real_pos = data.frame(x = c(105,105,115,115), y = c(27,20,27,20))
data_x = data.frame(img = img_pos$x, rel = real_pos$x)
data_y = data.frame(img = img_pos$y, rel = real_pos$y)
lm_x = lm(rel~img, data = data_x)
lm_y = lm(rel~img, data = data_y)
mytrans_x = function(myimg) {
predict(lm_x, newdata = data.frame(img = myimg))
}
mytrans_y = function(myimg) {
predict(lm_y, newdata = data.frame(img = myimg))
}
然后，再利用ImageJ軟件對中南郡的每個區域輪廓線單獨描邊采樣，這樣做的缺點就是兩個區域相鄰邊會有些不一致，出現小幅的咬合錯位現象，但這個對美觀影響不大。優點是大大節省時間。
把每個區域的邊界保存在單獨的文件中。然后在R中把這些數據轉化為GIS數據，保存為shp格式的標準地圖文件。關于代碼中函數的意義及范例（比我的代碼更清晰），請參考sp和maptools包的幫助文件。
library(maptools)
myfiles = c("Jiana.xls", "Kutedan.xls", "Miyaluo.xls", "Woda.xls", "Yada.xls")
mypolys = lapply(myfiles,
function(x) {
tmp = read.table(paste0("data/", x));
tmp = rbind(tmp, tmp[1, ]);
tmp$X = mytrans_x(tmp$X);
tmp$Y = mytrans_y(tmp$Y);
tmp
})
mynames = sub(".xls$", "", myfiles)
names(mypolys) = mynames
myPolygons = lapply(mynames,
function(x) {
tmp = mypolys[[x]];
Polygons(list(Polygon(cbind(tmp$X, tmp$Y))), x)
})
mySpn = SpatialPolygons(myPolygons)
myCNnames = c("嘉納", "庫特丹", "米亞洛", "沃達", "雅達")
myshpdata = SpatialPolygonsDataFrame(mySpn,
data = data.frame(
Names = mynames,
CNnames = myCNnames,
row.names = row.names(mySpn)))
# 我們要注意到：SpatialPolygonsDataFrame類的data成員的字段是可以自定義的，
# 這個是暴露給names函數以及$、[]運算符的。
writePolyShape(x = myshpdata, fn = "data/myDIYmap_poly")
這樣我們在就成功保存了shp格式的地圖文件（一共生成三個文件，一個shp文件，兩個輔助文件）。生成的地圖文件可以留給別人用，也可以正常打開繪圖了。
mydat = readShapePoly("data/myDIYmap_poly.shp")
plot(mydat)?
?
可以發現，在區域相鄰的邊界，有咬合分離現象，這是由于我們采樣的時候，每個區單獨描邊，產生了共享邊的不一致。不過，我們繪制地圖是為了展示流行病學數據，這個誤差是可以接受的。
library(ggplot2)
mysh = fortify(mydat, region = "CNnames")
names(mysh)[1:2] = c("x", "y")
myepidat = data.frame(id = unique(sort(mysh$id)))
myepidat$rand = runif(length(myepidat$id))
tmp = coordinates(mydat)
tmp = as.data.frame(tmp)
tmp$names = mydat$CNnames
ggplot(myepidat) + geom_map(aes(map_id = id, fill = rand), color = "white", map = mysh) +
geom_text(aes(x = V1,y = V2,label = names), family = "GB1", data = tmp)+
scale_fill_gradient(high = "red", low = "yellow") +
expand_limits(mysh) + coord_map()
?
?
如上，畫成統計地圖，還算美觀。
如果非要消除這種邊界交錯的不完美，就需要預先制定規劃，在位圖上分段采集邊界線，再拼接組合成區域輪廓。由于共享邊只采集一次，你能得到邊界完美的地圖。問題是，隨著地圖區域增多，你將在輪廓的拼接組合上，面臨幾何級數增長的復雜度。不過，離開現實的功利和脅迫，去追求完美，不也是推動這個世界前進的原動力么？
5 小結
盡管我在寫作中使用了這個星球上最強大的knitr軟件包來保證本文的可重復性，但是隨著官方新版數據在未來的發布，數據的字段名稱甚至組織布局將會有些變化，也會使本文代碼無法直接拷貝運行。還是希望讀者能自己掌握R，以無招勝有招。
喜歡讀統計之都主頁文章的結尾部分，因為常在此部分讀到作者“不著調”的話，發人深省。最愛楊燦兄改編的這段：
問：世間是否此山最高，或者另有高處比天高？
答：在世間自有山比此山更高，Open-mind要比天高。
參考文獻
- 謝益輝，2007，??http://yihui.name/cn/2007/09/china-map-at-province-level/??
- 邱怡軒，2009，??http://cos.name/2009/07/drawing-china-map-using-r/??
- 陳麗云，2011，??http://www.loyhome.com/用R畫（中國）地圖-2/??
- 寫長城的詩，2012，??http://www.r-bloggers.com/lang/chinese/1010??
- 楊燦，2011，??http://cos.name/2011/12/stories-about-statistical-learning??
附：本文所用地圖數據??下載??
本文轉自：??http://cos.name/2014/08/r-maps-for-china/#comment-6149??
