我在 R 里画地图（一）

前言

大约一个月前，接到一个任务：绘制全国的患病率地图，就草草做了一个初版，当时还做了一个 ppt 和组里的小伙伴分享。直到后来，才知道这种图有一个专业的称呼：choropleth maps，中文名字是分层设色图。可以看到在 google 中检索 choropleth maps 得到的结果：

单在画图这一方面讲，其实这是一个老生常谈的话题，也有数不胜数的工具和包。但问题的关键在于：绝大多数国外提供的中国地图并不规范。关于中国地图的规范问题，在姜大伟的知乎专栏使用中国地图的正确姿势中有比较翔实的介绍。我也根据此做成了其中一张幻灯片：

所以问题摆在了我们的面前：

找到一张可供使用的规范的中国矢量地图；
根据这张地图绘制我们想要的分层设色图。

提前提醒：本文代码直接复制粘贴不能够运行，请到文末下载分享的全部程序和地图文件！

寻找标准地图

在网络上，使用最多的是一份名为 bou2_4p.shp 的地图矢量文件，这份文件来自哪里已不可考，似乎是 2012 年国家提供的 1：400 万地理信息地图，但是后来又关闭了开放。这份地图应该是目前问题最少的地图，藏南、台湾、南海诸岛也都存在。但是时间已经来到了 2019 年，这份 12 年的地图是否那么无懈可击呢？

答案是否定的。根据 b 站 up 主“地理人_zxl”的视频 ArcGis更正老式中国基础地理信息数据错误，通过天地图和 bou2_4p 的比对，还是可以发现新疆的边界存在一些问题（红线为 bou2_4p，底图为标准地图）：

同时，bou2_4p 的地图也没有澳门特别行政区。通过上面视频中提供的方法，在 ArcGis 中对 bou2_4p 进行修改，就可以得到修正版本的 bou2_4p 了。这就是我们接下来要使用的标准地图。

使用 R 绘制地图

首先要说明，在我个人使用过的工具和包里面，还是 ArcGis 这种专业 GIS 软件绘制地图最为顺滑快捷。那为什么要使用 R 语言呢？其实还是出于以下这几个原因：

不只我要出图，别人也要进行地图制作；
免费开源，安装便捷，上手简单；
只需改动患病率数据，运行代码就可出图

再次提醒，以下我不会贴出完整代码，只会贴出一些核心的重点代码。完整的代码会分享在文章末尾。

基础绘图

首先，我们通过 rgdal::readOGR 读取修正过后的中国地图，然后对数据框进行转换，再和患病率数据 da 进行合并。

患病率数据 da 长什么样子呢？你一看便懂：

NAME	rate
安徽省	15.89
北京市	12.85
福建省	15.14
黑龙江省	NA

这里的患病率数据中用 rate 写明了每个省的患病率大小。然后通过 ggplot 绘制出我们需要的基础制图。

china_map <- rgdal::readOGR("map//bou2_4p.shp", use_iconv = TRUE, encoding = "UTF-8")
xdata <- china_map@data 
xs <- data.frame(xdata, id=as.character(seq(0, 924))) 
china_map1 <- fortify(china_map) 
china_map_data <- left_join(china_map1, xs, type = "full")
china_data <- left_join(china_map_data, da, by='NAME', type="full") 
  
ggplot(china_data, aes(long, lat))+
    geom_polygon(aes(group=group, fill = rate), color = 'white', size = 0.25)

我们可以在这一步对 x 轴和 y 轴的范围进行裁剪，修改地理坐标投影系，并修改绘图主题，隐藏 x 轴和 y 轴等等。下一步再进行进一步的美化。

美化底图

接下来就是进行地图的进一步美化，比如我们想要更改色阶的颜色。这里的 p1 就是我们上一步保存的底图。第二步的核心在于 scale_fill_gradientn，可以指定无限数目的渐变色了。而 na.value 则是指定缺失值省份的颜色。

p2 <- p1 +
  scale_fill_gradientn(
    colours=c("#4C79B0", "#9FD2AE", "#F4C4B2", "#BE1E56"),
    na.value = ’#D3D3D3'
  )
p2

顺带一提，想要选取好看的颜色，可以试试在 google 图片里面搜索 choropleth map。如果看到顺眼的颜色，可以直接得到色彩的 hex 值复制到自己的图里。

绘制九段线

九段线的地图文件来自 l9.shp，这里我们想要把九段线绘制在整个地图的右下角，并且小图中出现的广东等省份也是染好色的。思路就是把整张图和九段线全部画出来，再通过经纬度截取出我们想要的那部分：

l9 <- rgdal::readOGR('map\\l9.shp')
l91 <- fortify(l9)
  
china_map <- rgdal::readOGR("map//bou2_4p.shp", use_iconv = TRUE, encoding = "UTF-8")
xdata <- china_map@data 
xs <- data.frame(xdata, id=as.character(seq(0, 924))) 
china_map1 <- fortify(china_map) 
china_map_data <- left_join(china_map1, xs, type = "full")
china_data <- left_join(china_map_data, da, by='NAME', type="full") 
  
p9 <- ggplot(china_data, aes(long, lat))+
  geom_polygon(aes(group=group, fill = rate), color = "white", size = 0.25) +
  coord_map("albers", lat0=27, lat1=45) +
  scale_fill_gradientn(
    colours=plot_color,
    guide = "none"
  ) +
  geom_line(data=l91, aes(x=long, y=lat, group=group)) +
  coord_cartesian(xlim=c(105,125),ylim=c(3,30))

得到九段线小图后，再通过 grid 包把这个小图放置在原图的右下角。

library(grid)

vie <- viewport(width=0.225, height=0.15, x=0.78, y=0.15)
print(p9) 
print(l9, vp=vie)

这样，就达到了如下图这样的效果：

添加省份名字标签

在第一个版本做好之后，老板说让我可以加上省份名称的英文标签。我想，不如就再增加一个添加省份名称标签的功能，中英文都整上。这里，如果熟悉 ggplot 的朋友，应该会觉得很简单。只要准备好了两份数据文件，分别是 cn_text 和 en_text，里面记录每个省的标签名字和标签的位置。然后使用 geom_text 即可完成：

text_da <- read.csv("map\\cn_text.csv")
p3 <- p2 +
  geom_text(data=text_da, aes(label=text))

在 geom_text 里面，修改文字的字体，大小，颜色自然也不在话下。除此之外，只要把 geom_text 替换成 geom_shadowtext，就可以使用带白边的标签，效果如下：

在我的破电脑上，如果单次开始运行整个绘图程序，到出图保存，大概需要 10s 左右。这可能是一个劣势，当然也不是不能接受的。

后记

想来这篇也拖了很久的时间，着实不太应该。本来想把所有内容写在一篇内容中，但是想想我也不易写，读者也不易读，实在是不太好。于是我只把最最核心的内容记录在博客中，主要是陈述思路，而删减了很多内容。

这里，我把所有地图文件，程序和说明文档打包在了这里。我已经尽量做到简单易用，这个系列的第二篇就写写怎么画插值热力图。之后我还会继续研究绘制地图的相关方法，有新的技巧，感想和技能还会更新。

XZY's BLOG