R語言做地圖上的分析

R和ggplot可視化功能非常強大，了解了一下其中的地圖做法，發(fā)現(xiàn)R做世界地圖、美國地圖非常容易，但做中國地圖就太麻煩了，需要自己DIY。

DIY也有多種方式，但網(wǎng)絡(luò)上各種帖子教程的出圖效果都不太理想，達(dá)不到工作用要求。下面是我的摸索過程，記錄如下備忘，也請教于R老師們。

0、引子

R里有個繪制地圖的maps包，加載后即可繪制地圖，試一下以下命令：

library(maps)

map()

即可畫出一個世界地圖。再試試：

map("state")

map("county")

可分別畫出美國分州和分county的地圖，真是不費吹灰之力。

可是，要畫一幅中國地圖，就沒有這么容易了，需要先加載一個mapdata包：

library(mapdata)

map("china")

但發(fā)現(xiàn)居然還是沒有重慶的地圖，太坑爹了，沒法用，只有自己構(gòu)造中國地圖了。下面開始進(jìn)入正題，如何用R繪制中國分省熱力地圖。

1、準(zhǔn)備地圖數(shù)據(jù)

先要找到中國地圖數(shù)據(jù)文件。到國家基礎(chǔ)地理信息中心的網(wǎng)站（http://nfgis.nsdi.gov.cn）下載官方空間文件，但這個政府網(wǎng)站總是打不開！你可到微信公眾號 iamExcelPro 發(fā)送 shapefile 獲得下載地址，我是雷鋒。解壓到 c:/rstudy 目錄，有3個文件，都是需要的。

加載maptools包，讀取空間文件：

library("maptools")

china_map = readShapePoly("c:/rstudy/bou2_4p.shp")       # 讀取地圖空間數(shù)據(jù)

plot一下看看，是一幅中國地圖，有重慶，還包括南海的島嶼，政府?dāng)?shù)據(jù)就是嚴(yán)謹(jǐn)一些：

plot(china_map)

但地圖投影方式不對，看起來太扁了，完全不是我們常見的昂首雄雞狀。

加載ggplot2包，用ggplot繪制，并使用polyconic投影方式，顯示正常。

library(ggplot2)

ggplot(china_map,aes(x=long,y=lat,group=group)) +

    geom_polygon(fill="white",colour="grey") +

    coord_map("polyconic")

現(xiàn)在地圖是可用的了，但還需要加載和拼接行政信息，以便能與業(yè)務(wù)數(shù)據(jù)映射。

x <- china_map@data           #讀取行政信息

xs <- data.frame(x,id=seq(0:924)-1)          #含島嶼共925個形狀

library(ggplot2)

china_map1 <- fortify(china_map)           #轉(zhuǎn)化為數(shù)據(jù)框

library(plyr)

china_map_data <- join(china_map1, xs, type = "full")       #合并兩個數(shù)據(jù)框

提示：Joining by: id

看不懂？沒關(guān)系，過了就行。

2、準(zhǔn)備業(yè)務(wù)數(shù)據(jù)

網(wǎng)上教程居然都是在命令行里輸入數(shù)據(jù)，也很坑爹。我們還是從Excel表格轉(zhuǎn)存來得方便。

按以下格式準(zhǔn)備好指標(biāo)數(shù)據(jù)，并存為csv格式文件。不直接讀取xlsx文件是因為需要裝的包比較麻煩。

注意第1列的字段名為NAME，各省名稱也是要固定一致的，是為了和地圖數(shù)據(jù)框里的省名一致，便于合并。各省名稱是用以下命令查看并記下的。

> unique(china_map@data$NAME)

[1] 黑龍江省         內(nèi)蒙古自治區(qū)     新疆維吾爾自治區(qū) 吉林省

[5] 遼寧省           甘肅省           河北省           北京市

[9] 山西省           天津市           陜西省           寧夏回族自治區(qū)

[13] 青海省           山東省           西藏自治區(qū)       河南省

[17] 江蘇省           安徽省           四川省           湖北省

[21] 重慶市           上海市           浙江省           湖南省

[25] 江西省           云南省           貴州省           福建省

[29] 廣西壯族自治區(qū)   臺灣省           廣東省           香港特別行政區(qū)

[33] 海南省

下面讀取業(yè)務(wù)指標(biāo)數(shù)據(jù)，并與地圖數(shù)據(jù)合并：

mydata <- read.csv("c:/rstudy/geshengzhibiao.csv")          #讀取指標(biāo)數(shù)據(jù)，csv格式

china_data <- join(china_map_data, mydata, type="full")          #合并兩個數(shù)據(jù)框

提示：Joining by: NAME

3、繪制地圖

現(xiàn)在可以開始試試畫填色地圖了：

ggplot(china_data, aes(x = long, y = lat, group = group, fill = zhibiao)) +

    geom_polygon(colour="grey40") +

    scale_fill_gradient(low="white",high="steelblue") +  #指定漸變填充色，可使用RGB

    coord_map("polyconic")        #指定投影方式為polyconic，獲得常見視角中國地圖

好，看到填色地圖了，但圖中的背景色、坐標(biāo)軸、經(jīng)緯線都是不需要的，圖例也可以放到左下角，用theme命令清除：

ggplot(china_data, aes(x = long, y = lat, group = group,fill = zhibiao)) +

    geom_polygon(colour="grey40") +

    scale_fill_gradient(low="white",high="steelblue") +  #指定漸變填充色，可使用RGB

    coord_map("polyconic") +       #指定投影方式為polyconic，獲得常見視角中國地圖

    theme(               #清除不需要的元素

         panel.grid = element_blank(),

         panel.background = element_blank(),

         axis.text = element_blank(),

         axis.ticks = element_blank(),

         axis.title = element_blank(),

         legend.position = c(0.2,0.3)

         )

4、添加省名標(biāo)簽

有時候需要顯示省名標(biāo)簽，特別是給老領(lǐng)導(dǎo)看?？筛鶕?jù)每個省形狀的經(jīng)緯度平均值求近似中心位置，標(biāo)注省名。

midpos <- function(x) mean(range(x,na.rm=TRUE))#取形狀內(nèi)的平均坐標(biāo)

centres <- ddply(china_data,.(province),colwise(midpos,.(long,lat)))

ggplot(china_data,aes(long,lat))+       #此處語法與前面不同，參考ggplot2一書P85

    geom_polygon(aes(group=group,fill=zhibiao),colour="black")+

    scale_fill_gradient(low="white",high="steelblue") +

    coord_map("polyconic") +

    geom_text(aes(label=province),data=centres) +

    theme(

         panel.grid = element_blank(),

         panel.background = element_blank(),

         axis.text = element_blank(),

         axis.ticks = element_blank(),

         axis.title = element_blank()

         )

但發(fā)現(xiàn)海南兩個字跑到南海去了，不行。下面改使用省會城市的經(jīng)緯度數(shù)據(jù)標(biāo)注省名。

province_city <- read.csv("c:/rstudy/chinaprovincecity.csv")  #讀取省會城市坐標(biāo)

ggplot(china_data,aes(long,lat))+

    geom_polygon(aes(group=group,fill=zhibiao),colour="grey60")+

    scale_fill_gradient(low="white",high="steelblue") +

    coord_map("polyconic") +

geom_text(aes(x = jd,y = wd,label = province), data =province_city)+

    theme(

         panel.grid = element_blank(),

         panel.background = element_blank(),

         axis.text = element_blank(),

         axis.ticks = element_blank(),

         axis.title = element_blank()

         )

至此，終于DIY繪制出了一份中國分省的熱力地圖，真是夠折騰夠找虐的，好在圖形很精準(zhǔn)，以后也可以復(fù)用代碼。

所以結(jié)論就是，一般非專業(yè)用戶還是使用BingMap、PowerMap、Tableau或者《用地圖說話》中的Excel模板，直接填數(shù)據(jù)出地圖吧，其中Excel模板方式是最簡單、便攜，office協(xié)同性最好的。

如果你想用這些代碼出圖，可準(zhǔn)備好用到的3個文件到相應(yīng)目錄，其中一個是你的業(yè)務(wù)指標(biāo)csv文件，一步步運行以上代碼，應(yīng)該就可得到一幅中國地圖。