流行病学中发病情况的地区分布图是很常见的一类做图，R中的ggplot2是一个非常优雅的绘图工具，通过原始的地理测绘数据配合ggplot2绘制地图是比较硬核的一种方法，国家测绘中心早年免费提供地理测绘数据，最近的免费数据是2004年的测绘数据( https://uploads.cosx.org/2009/07/chinaprovinceborderdata_tar_gz.zip )[1]，2014年新版的测绘数据则是收费的。两者最直观的区别在于重庆。如果仅是练手可以使用2004年的测绘数据。
流行病学中比较常用的是气泡图，这个在网上已经有很多教程。但用这种方法绘制地图时直接按省份进行染色的教程竟然一篇都没有，摸索了一下，绘制出了相对完美的分布图（发病人数使用的是1月25日的数据）：
plot_zoom.png
这个地图的制图难点在于原始的测绘数据虽然是34个省市自治区的（原数据实际没有标出澳门），但因为有的省份面积太大，所以一共花了924个数据用于描述这34个省市自治区，这被称为描述层。为了画出描述层，一共有9万多条几何层的数据，这些几何层的数据是一些坐标点，相互连接起来之后就画出了各省份的边界，而同一组描述层内的几何层数据会拼接成一个省份。
如果是普通的34个数据，那是最容易绘图的了，直接生成34个颜色就可以，但原始的数据是个S4类数据，里面的data子数据包含了924个描述层数据，通过fortify函数解析这个S4数据可以获得9万多条几何层的绘图数据。R自带的plot函数是直接通过描述层调用几何层，因此需要924个颜色去填充，正是通过一篇plot函数绘图的文章让我想到了ggplot2的做图方法[2]。ggplot2中是直接按9万多条原始绘图数据去制图的，所以需要9万多个颜色，我直接套用了原文中生成颜色的函数，做出来的图已经可以用于展示。代码中已添加了详细注释：
# 加载包
library(ggplot2) # 绘图
library(rgdal) # 读取地图数据
library(plyr) # 数据处理
# 全局参数
options(stringsAsFactors = FALSE) # 这个改了省心，r喜欢按因子读取字符
windowsFonts("Times New Roman" = windowsFont("Times New Roman")) # 设置字体
# 读取文件，设置工作路径
setwd("E:/Store/code/03.R/03.ggplot2/06.map/test") # 根据实际情况改自己的工作路径
vData <- readOGR("E:/Store/code/03.R/code/06.maps/map/2014/maps/bou2/bou2_4p.shp") # 按地图文件的实际位置修改
# 数据处理
vTemp <- vData@data; # S4类数据,用@取子集,得到的为描述层
vDescribeData <- data.frame(vTemp, id = seq(0:924) - 1) # 为描述层添加id，方便后续内联数据
vFormatData <- fortify(vData) # 得倒的为几何层数据
vMapData<-join(vFormatData, vDescribeData, type = "full") # 合并描述层和几何层数据（普通做图这一步并不是必须的，通常只需要几何层就可以做出地图，合并起来更直观一些）
# 读取各省份的发病人数
vNum <- read.csv("num.csv", stringsAsFactors = FALSE)
vProvName <- vNum$vProvName # 省份名称
vNum <- vNum$vNum # 发病人数
vProvCol = vNum # 将发病人数传给颜色，用于表示颜色的深浅
# 颜色函数，从他人博客套用的函数，用处是给每一个描述层生成一个颜色（描述层有924条数据）
fGetColor <- function(mapdata, vProvName, vProvCol, othercol){
f = function(x, y) ifelse(x %in% y, which(y == x), 0)
colIndex = sapply(mapdata@data$NAME, f, vProvName)
fg = c(othercol, vProvCol)[colIndex + 1]
return(fg)
}
# 使用自定义的颜色函数，为描述层生成颜色
vCol <- fGetColor(vData, vProvName, vProvCol, 0)
# 为描述层颜色添加id
vRealCol <- data.frame(vCol, id = seq(0:924) - 1)
# 通过id将描述层的颜色数据和之前的几何层地图数据内联，简单的说就是将924条颜色数据根据id映射到9万多条几何层数据
vMapData <- join(vMapData, vRealCol, type = "full")
# ggplot2常规做图，通过scale_fill_continuous函数产生渐变色
ggplot(vMapData, aes(x = long, y = lat, group = group, fill = as.numeric(vMapData$vCol)))+
geom_polygon() +
scale_fill_continuous(low = "white", high = "red") +
geom_path(colour = "black") +
ggtitle("新型肺炎确诊病例分布示意图") +
theme(plot.title = element_text(lineheight = 0.8, face = "bold", size = 20)) +
theme(axis.text = element_text(size = 15, colour = "black", family = "Times New Roman"),
axis.title = element_blank()) +
# specify tick marks
coord_map(xlim = c(73, 136), ylim = c(5, 54)) +
scale_y_continuous(breaks = seq(0, 60, 10)) +
scale_x_continuous(breaks = seq(70, 135, 10)) +
# remove guide
guides(fill = FALSE)
其中num.csv的格式如下：
num.png
下载地址为：http://www.1x1y.top/books/lib/exe/fetch.php?media=num.zip
vNum2列是真实发病人数，我用了渐变色会导致人数较低时偏白色，所以我对发病人数做了一点转换，实际绘图用的是转换后的vNum列数据= =。将发病人数生成新的一列评级数据作为染色标准时做出的图会更好。由于数据转换以及添加地名比较简单，我就没有再接着做了。只要再配合一个爬虫爬取发病人数，就可以实时绘制最新的疫情分布图了。地图文件参考文中给出的链接下载。
参考文献：
[1] 邱怡轩, https://cosx.org/2009/07/drawing-china-map-using-r/.2009.
[2] [阿蛮的杜鹃], https://www.cnblogs.com/lizhilei-123/p/6734378.html.2017.
网友评论