COVID-19学术科普专题：绘制疫情分布地图_

提及空间数据可视化，人们对ArcGIS，Qgis耳熟能详，但R语言对于空间数据的处理同样也十分强大，常用的语言包有sp，maptools，rgdal，ggplot2等。本文以广东省新冠肺炎疫情数据为例（截至4月14日），用R语言绘制疫情地图，识别高中低风险地区，协助开展分区分级管理工作。

R绘制疫情地图的四个步骤

Step 1

准备地图

下载和读取地图

广东省的地图可以从国家基础地理信息中心（http://www.ngcc.cn/ngcc）或者资源环境数据云平台（ http://www.resdc.cn）免费下载。下载的地图数据包括存储矢量数据的空间信息文件.shp，存储矢量数据的属性信息文件.dbf，要素几何特征的索引文件.shx，以及存储坐标系统的文件.prj。下载的地理数据文件如下图：

图1 下载的地图文件

地图数据的读取代码和结果如下所示：

library(rgdal)

library(maptools)

library(sp)

#有时因编码问题，需要加上代码use_iconv=TRUE, encoding = "GBK”进行转码

图2 地图数据的读取结果

转为数据框地图

下载的地图是一个sp包定义的SpatialPolygons DataFrame类，并不是在R常用的data.frame类型，因此需要用fortify函数将地图转为data.frame类型，才能在用ggplot2画图。代码如下：

library(ggplot2)

map

ggplot(map, aes(x = long, y = lat, group = group))+ geom_path()+ coord_fixed()

绘制的广东省地图展示如下：

图3 fortify后绘制的广东省地图

再将转换后的地图与原来数据合并，代码如下：

#转为数据框地图

GD@data$id

#数据框与原来数据合并

GDnew

#地图转为了data.frame类型

class(GDnew)

Step 2

整理疫情数据

疫情数据包括城市中心的名称、经纬度（可百度或数据框地图中获取）以及确诊病例数（若需要，还可以增加死亡数、出院数等）。数据格式如下所示，CityCN为城市名，其中Lng为经度，Lat为纬度，Case为该城市的累计病例数。

CityCN
Lng
Lat
Case
广州市
113.5345876
23.35273052
488
深圳市
114.1310512
22.65540855
458
珠海市
113.2515266
22.16070788
103
东莞市
113.8758797
23.01688407
100
......
云浮市
111.7934509
22.91578903

Step 3

合并地理数据和疫情数据

我们的目的是将疫情数据展示在地图上，因此需要将疫情数据整合到地理数据。这里，我们通过两者共有的城市名称（也可以是城市代码或其他变量）用left_join进行合并。代码如下所示：

library(dplyr)

GDnew1

由此，我们通过两个数据共有的城市名称进行链接，合并后的数据不但包括地图数据，也包括了疫情数据（红色方框）。合并后数据如下图所示。

图4 合并后数据展示

Step 4

使用ggplot2画图

数值型变量

在合并后的数据里，确诊病例数“Case”是数值型变量（numeric），最大值是488，最小值是0。若直接画图，代码（暂不美化）及输出结果如下：

library(RColorBrewer)

ggplot() +

geom_polygon(data=GDnew1, aes(x = long, y =lat, group = group, fill= GDnew1$Case),colour = "white")+

scale_fill_gradient(low = "rosybrown1",high = "coral2")

图5 数值型变量绘制的地图

由图可见，除了广州和深圳外的其他城市都是鲜红色，为了能更好展示差异，可以考虑重新设置色带（未来可视化专题将介绍该部分内容），或是将病例数转化为因子型变量（factor)再绘图。

因子型变量

首先，使用cut函数将数值型的Case变成因子型变量，代码具体如下：

GDnew1$Case_fac

labels=c("0","1~9","10~14","15~24","25~49","50~99",">100"),

include.lowest=TRUE, right=FALSE)

然后，用ggplot2绘图并美化，变成如下结果：

图6 因子型变量绘制的地图

主要代码如下：

ggplot() +

geom_polygon(data=GDnew1, aes(x = long, y =lat, group = group, fill= GDnew1$Case_fac),colour = "white")+

scale_fill_brewer(palette="Reds")+#离散型调色

guides(fill = guide_legend(title='Confirmed Case',reverse=TRUE))+ #图例顺序倒置

coord_fixed()+

labs(x=NULL, y=NULL) +

geom_text(aes(x = Lng,y = Lat), label = paste(case$CityEN), size = 4.5,data = case)+#显示城市名

geom_text(aes(x = Lng,y = Lat-0.04), label = paste("n",case$Case), size = 4.5,data = case)+#显示病例数

theme(panel.border = element_blank(),.....)#优化主题

随着新冠肺炎疫情态势向好，开展地区疫情风险评估，对落实分区分级管理，精准推进复工复产具有重要意义。下图是基于广东省卫生健康委发布的广东省新冠肺炎疫情风险等级分区分级名单（2020年4月5日），制作全省分区县的风险等级地图。风险分为高中低三个级别，名单中无高风险地区，中风险地区为越秀区、白云区、宝安区和惠来县，其余为低风险地区。绘制步骤与上述因子型变量一致，不再赘述。