基于ArcGIS和R语言的世界疫情分析

“随着新冠肺炎疫情发展，目前国内疫情得到了初步控制，世界范围的疫情越来越严重，在2月28日，世卫组织将新冠肺炎疫情风险从“高”，调整为“非常高”，新冠肺炎已经开始在世界范围内扩散，防控新冠肺炎，成为了全人类的共同任务”

本文使用ArcGIS和R语言，对世界新冠肺炎疫情进行分析。

01

ArcGIS制作世界疫情分布图

—

先看一下成图，以新增确诊病例为例，制作世界各国新增确诊病例分布图。

然后看一下整理的数据，数据分成两部分，疫情属性数据和世界矢量数据，疫情表格中，有“国家”字段，和矢量数据中“CH_NAME”字段对应，可以通过挂接的方式，将两组数据连接起来。

数据挂接完成后，就可以对数据进行符号化了，在这里采用了分级设色的方法，对每个国家的新增病例进行符号化。

然后对国家名称和新增病例数进行标注，具体方法不再赘述，请参阅：

这样，在ArcGIS中的制图工作就完成了。

02

R语言进行世界疫情分析

—

重点在这一部分。

首先加载所需的一些程序包

读取EXCEL文件用的readxl包

时间序列用的lubridate包

绘图神器ggplot2包

格网grid包

数据处理dplyr包

library(readxl)library(lubridate)library(ggplot2)library(grid)library(dplyr)读取数据，由于数据有csv格式的，也有xls格式的，分成两段读取，一部分读csv，一部分读xls，在这里我将csv文件和xls文件分别放在了worldDay1和worldDay2文件夹下。
为了生成时间序列，我最后使用ymd函数生成了一个R语言可以识别的时间序列。

a = list.files("worldDay1") #list.files命令将input文件夹下所有文件名输入adir = paste("./worldDay1/",a,sep="") #用paste命令构建路径变量dirn = length(dir) #读取dir长度，也就是文件夹下的文件个数mergedata1 = read.csv(file = dir[1],header=T,sep=",") #读入第一个文件内容（可以不用先读一个，但是为了简单，省去定义data.frame的时间，我选择先读入一个文件。for (i in 2:n){ new.data = read.csv(file = dir[i], header=T, sep=",") #读取CSV文件 mergedata1 = rbind(mergedata1,new.data)}
a = list.files("worldDay2") #list.files命令将input文件夹下所有文件名输入adir = paste("./worldDay2/",a,sep="") #用paste命令构建路径变量dirn = length(dir) #读取dir长度，也就是文件夹下的文件个数mergedata2 = read_excel(path = dir[1],sheet = 1, col_names = T) #读入第一个文件内容（可以不用先读一个，但是为了简单，省去定义data.frame的时间，我选择先读入一个文件。for (i in 2:n){ new.data = read_excel(path = dir[i],sheet = 1, col_names = T) #读取EXCEL格式文件 mergedata2 = rbind(mergedata2,new.data)}MergeData <- rbind(mergedata1, mergedata2)#由于读取完的数据存在 空行，我使用na.omit函数去除空行MergeData <- na.omit(MergeData)Daydate <- ymd(MergeData$时间)WorldData <- cbind(MergeData, Daydate)

接下来让R语言自动输出一些结论：

使用paste函数可以将文字和计算结果同时打印输出出来，如下所示：

WorldData5 <- WorldData[WorldData$Daydate==max(WorldData$Daydate), ]worldData5Yestoday <- WorldData[WorldData$Daydate==max(WorldData$Daydate)-1, ]
paste("截至", max(WorldData$Daydate),"全球新增新冠肺炎病例",sum(WorldData5$新增确诊),"较昨日新增病例",sum(WorldData5$新增确诊)-sum(worldData5Yestoday$新增确诊))paste("其中中国", WorldData5$新增确诊[WorldData5$国家=="中国"])paste("外国合计", sum(WorldData5$新增确诊[WorldData5$国家!="中国"]))paste("报告新增肺炎病例的国家有",nrow(WorldData5[WorldData5$新增确诊!=0, ]) ,"个")paste("较",max(WorldData$Daydate)-1, "减少/增加", nrow(WorldData5[WorldData5$新增确诊!=0, ])-nrow(worldData5Yestoday[worldData5Yestoday$新增确诊!=0, ]) )paste("目前已有",nrow(WorldData5) ,"个国家报告新冠肺炎病例", "有", nrow(WorldData5)- nrow(worldData5Yestoday), "个国家新发现新冠肺炎确诊病例")

绘制报告新增确诊病例的国家数变化折线图，并输出一段描述：

p1 <- ggplot(WorldData2, aes(Daydate, `n()`))+
  labs(x="时间", y="新增确诊国家数")+
  theme(axis.title = element_text(size = 20),    #调整标题大小
        axis.text.x = element_text(size = 18),    #x轴标签大小
        axis.text.y = element_text(size = 18))+     #y轴标签大小
  geom_line()+
  geom_point()
p1
paste(max(WorldData2$Daydate), "共有", WorldData2$`n()`[WorldData2$Daydate==max(WorldData2$Daydate)],
      "个国家报告发现新增确诊新冠肺炎病例，和", max(WorldData2$Daydate)-1, "相比增加/减少", 
      (WorldData2$`n()`[WorldData2$Daydate==max(WorldData2$Daydate)])-(WorldData2$`n()`[WorldData2$Daydate==(max(WorldData2$Daydate)-1)]))

绘制中国新增确诊病例变化图：

WorldData3 <- WorldData[WorldData$国家!="中国", ]
WorldData4 <- WorldData[WorldData$国家=="中国", ]
p2 <- ggplot(WorldData4, aes(Daydate, 新增确诊, color=国家))+
  labs(x="时间", y="新增确诊病例数", title="中国新增确诊病例变化图")+
  theme(title = element_text(size = 22),
        axis.title = element_text(size = 20),    #调整标题大小
        axis.text.x = element_text(size = 18),    #x轴标签大小
        axis.text.y = element_text(size = 18))+     #y轴标签大小
  geom_line()+
  geom_point()
p2
paste(max(WorldData4$Daydate), "中国新增确诊病例", 
      WorldData4$新增确诊[WorldData4$Daydate== max(WorldData4$Daydate)], "例,较",max(WorldData4$Daydate)-1,
      "新增确诊病例",WorldData4$新增确诊[WorldData4$Daydate== (max(WorldData4$Daydate)-1)],
      "例,增加/减少", 
      WorldData4$新增确诊[WorldData4$Daydate== max(WorldData4$Daydate)]-WorldData4$新增确诊[WorldData4$Daydate== (max(WorldData4$Daydate)-1)],
      "例")

绘制其它国家新增确诊病例变化图：

p3 <- ggplot(WorldData3, aes(Daydate, 新增确诊, color=国家))+
  labs(x="时间", y="新增确诊病例数", title="其它国家新增确诊病例变化图")+
  theme(title = element_text(size = 22),
        axis.title = element_text(size = 20),    #调整标题大小
        axis.text.x = element_text(size = 18),    #x轴标签大小
        axis.text.y = element_text(size = 18), legend.position = "bottom")+     #y轴标签大小
  geom_line()+
  geom_point()
p3

上面绘制了好几个折线图，接下来这个要换换花样了，绘制一个新增确诊病例国家所在洲的饼图：

饼图稍微复杂一些，需要的数据也和前面不太一样。

前文提到了ArcGIS制作的疫情分布图，那么，各洲分布从哪来呢？没错，就从上文提到的挂接后的数据中来。

全球矢量数据存在一个CONTINENT字段，根据这个字段，即可统计各个洲情况。

在ArcMap中将属性表导出

在这里提供一个窍门，为了导出的数据处理读取比较方便，可以直接输出文本文件，并通过自行指定文件后缀名的方式，让ArcMap导出为csv文件

导出后的文件用R读取，可能会存在乱码的情况，这种情况下只需要用NOTEPAD++将CSV文件编码改为ANSI码即可。

导出CSV文件后，就可以绘制饼图了：

下面数据处理采用了dplyr包中的函数，具体使用方法，推荐一个公众号：

可以参阅下面公众号的文章

#新增确诊国家各洲分布
GlobalData <- read.csv(file = "全球0306.csv", header = T)
GlobalData <- na.omit(GlobalData)
GlobalData2 =
  GlobalData%>%
  filter(新增确诊!=0)%>%
  group_by(CONTINENT)%>%
  summarise(n())
ggplot(GlobalData2,aes(x = factor(1),`n()`,fill=factor(CONTINENT)))+
  geom_bar(stat="identity",position="fill")+
  coord_polar(theta="y")+
  labs(x = '', y = '', title = '') + 
  theme(axis.text = element_blank())+
  guides(fill= guide_legend(title = "新增病例所在洲"))+
  geom_text(stat="identity",aes(label = `n()`), size=4, position=position_fill(vjust = 0.5), color= "black")

绘制饼图的时候，关于饼图标注，要特别注意一下position参数，如下图所示，标注的参数需要和饼图的前身（柱状图bar）对应。

到这里文章就结束了，如果觉得不错还请转发或者点击“在看”哦

欢迎打赏和留言

本来定的定时群发，结果一看，明天3月8日了，最后留个小彩蛋

最近一直都在忙着分析数据写报告，有时候忙的饭都来不及做

这几天媳妇给我做了好几顿饭

一定要好好夸夸她，下面是她做的可乐鸡翅，很好吃。

再次祝福广大的劳动妇女们节日快乐！