R语言做地图上的分析

R和ggplot可视化功能非常强大，了解了一下其中的地图做法，发现R做世界地图、美国地图非常容易，但做中国地图就太麻烦了，需要自己DIY。

DIY也有多种方式，但网络上各种帖子教程的出图效果都不太理想，达不到工作用要求。下面是我的摸索过程，记录如下备忘，也请教于R老师们。

0、引子

R里有个绘制地图的maps包，加载后即可绘制地图，试一下以下命令：

library(maps)

map()

即可画出一个世界地图。再试试：

map("state")

map("county")

可分别画出美国分州和分county的地图，真是不费吹灰之力。

可是，要画一幅中国地图，就没有这么容易了，需要先加载一个mapdata包：

library(mapdata)

map("china")

但发现居然还是没有重庆的地图，太坑爹了，没法用，只有自己构造中国地图了。下面开始进入正题，如何用R绘制中国分省热力地图。

1、准备地图数据

先要找到中国地图数据文件。到国家基础地理信息中心的网站（http://nfgis.nsdi.gov.cn）下载官方空间文件，但这个政府网站总是打不开！你可到微信公众号 iamExcelPro 发送 shapefile 获得下载地址，我是雷锋。解压到 c:/rstudy 目录，有3个文件，都是需要的。

加载maptools包，读取空间文件：

library("maptools")

china_map = readShapePoly("c:/rstudy/bou2_4p.shp")       # 读取地图空间数据

plot一下看看，是一幅中国地图，有重庆，还包括南海的岛屿，政府数据就是严谨一些：

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")

现在地图是可用的了，但还需要加载和拼接行政信息，以便能与业务数据映射。

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

xs <- data.frame(x,id=seq(0:924)-1)          #含岛屿共925个形状

library(ggplot2)

china_map1 <- fortify(china_map)           #转化为数据框

library(plyr)

china_map_data <- join(china_map1, xs, type = "full")       #合并两个数据框

提示：Joining by: id

看不懂？没关系，过了就行。

2、准备业务数据

网上教程居然都是在命令行里输入数据，也很坑爹。我们还是从Excel表格转存来得方便。

按以下格式准备好指标数据，并存为csv格式文件。不直接读取xlsx文件是因为需要装的包比较麻烦。

注意第1列的字段名为NAME，各省名称也是要固定一致的，是为了和地图数据框里的省名一致，便于合并。各省名称是用以下命令查看并记下的。

> unique(china_map@data$NAME)

[1] 黑龙江省         内蒙古自治区     新疆维吾尔自治区 吉林省

[5] 辽宁省           甘肃省           河北省           北京市

[9] 山西省           天津市           陕西省           宁夏回族自治区

[13] 青海省           山东省           西藏自治区       河南省

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

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

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

[29] 广西壮族自治区   台湾省           广东省           香港特别行政区

[33] 海南省

下面读取业务指标数据，并与地图数据合并：

mydata <- read.csv("c:/rstudy/geshengzhibiao.csv")          #读取指标数据，csv格式

china_data <- join(china_map_data, mydata, type="full")          #合并两个数据框

提示：Joining by: NAME

3、绘制地图

现在可以开始试试画填色地图了：

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命令清除：

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、添加省名标签

有时候需要显示省名标签，特别是给老领导看。可根据每个省形状的经纬度平均值求近似中心位置，标注省名。

midpos <- function(x) mean(range(x,na.rm=TRUE))#取形状内的平均坐标

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()

         )

但发现海南两个字跑到南海去了，不行。下面改使用省会城市的经纬度数据标注省名。

province_city <- read.csv("c:/rstudy/chinaprovincecity.csv")  #读取省会城市坐标

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绘制出了一份中国分省的热力地图，真是够折腾够找虐的，好在图形很精准，以后也可以复用代码。

所以结论就是，一般非专业用户还是使用BingMap、PowerMap、Tableau或者《用地图说话》中的Excel模板，直接填数据出地图吧，其中Excel模板方式是最简单、便携，office协同性最好的。

如果你想用这些代码出图，可准备好用到的3个文件到相应目录，其中一个是你的业务指标csv文件，一步步运行以上代码，应该就可得到一幅中国地图。