我在 R 里画地图(一)

前言

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

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

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

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

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

寻找标准地图

在网络上,使用最多的是一份名为 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

006y8mN6gy1g79r767fk6j31by0ts40r.jpg

顺带一提,想要选取好看的颜色,可以试试在 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) 

这样,就达到了如下图这样的效果:

006y8mN6gy1g79rfarviwj31bw0tgjtr.jpg

添加省份名字标签

在第一个版本做好之后,老板说让我可以加上省份名称的英文标签。 我想,不如就再增加一个添加省份名称标签的功能,中英文都整上。 这里,如果熟悉 ggplot 的朋友,应该会觉得很简单。 只要准备好了两份数据文件,分别是 cn_texten_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,就可以使用带白边的标签,效果如下:

006y8mN6gy1g79rg0c73zj31b80tgdir.jpg

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

后记

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

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