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