利用丁香园数据生成疫情分布地图（R语言）| 博文精选

来源 | CSDN 博客

作者 | 万里写入胸怀间

责编 | Carol

出品 | CSDN云计算（ID：CSDNcloud）

疫情牵动大家，除了做好分内工作，管好自己不给社会添乱，也就是只能持续关注疫情了。现在各大门户平台都上线了疫情实时地图，但基本上数据都只展示到省这个级别。

当然现在各地存在感染病例的小区也都及时披露了，相信疫后应该会有高精度的传染路径、传染扩散分布的分析。就目前而言，丁香园上的数据应该是目前聚合数据里面精度最高的了，可以获取到设区市级别的数据。但丁香园的数据存在地名不规范的问题，以及省、市、县三级行政区划并存的问题，不进行清洗不便做图。

底图制备

一张合适的底图是制作地图的基础。与行政区划相关的底图，最权威的数据源当然是从自然资源部的数据中心下载相关矢量数据再进行加工，最近微信公众号数读城市就在做这个浩大的工程。

我的方法是用高德API获取行政区划边界数据然后自行制图。高德数据的优点在于可以获取从全国、省、市、区（直辖市可以到镇）各级别行政区划边界的数据，对于不同级别行政区数据并存的情况，可以按需灵活组合，快速生成需要的底图。

丁香园公布的数据以设区市为主，台湾、香港、澳门公布总体数字，县（市）的数据较少。根据数据特点，我决定省级行政区划就保留港、澳、台三地，其余31个省、市、自治区保留设区市、自治州、省代管县级市以及直辖市的区、县。也就是高德API中，各省、自治区、直辖市子一级的行政单元。我是用sf包进行处理的，底图对象叫china_city数据框结构如下。

adcode即城市编码，下一步里从丁香园上下载的数据里有个locationId,和它是同一个意思。高德的adcode基于民政部的行政区划代码，并在民政部代码的基础上加了市辖区一级的代码，会干扰后续的处理，需要把它清理掉。

Amap_adcode=readxl::read_xlsx('AMap_adcode_citycode.xlsx') %>% 
            select(name=中文名,adcode=adcode) %>% mutate_at('adcode',as.numeric) %>% 
            filter(!grepl('辖区',name))

此外，坐标我也转换成wgs84。

疫情数据读取

方法来自于《R语言获取丁香园疫情数据动态网站》这篇帖子。这里就简单展示一下代码。

library(pacman)
p_load(sf,tidyverse,rvest,jsonlite)
setwd() #这里要自行设定工作路径
#read web
system('~/2019_nConv/phantomjs-2.1.1-windows/phantomjs-2.1.1-windows/bin/phantomjs.exe my_js_dxy.js')
url='dxy_web.html'
withjs=read_html(url) %>% html_nodes('script#getAreaStat') %>% html_text()
mydata=str_remove_all(withjs,'try \\{ window.getAreaStat = ') %>% str_remove_all('\\}catch\\(e\\)\\{\\}')%>% fromJSON()

《R语言获取丁香园疫情数据动态网站》：

https://blog.csdn.net/yuanzhoulvpi/article/details/104099382

以上和链接帖子里完全一样，但我要获取的是城市的数据，所以还要再加一行：

city_data=mydata$cities %>% rlist::list.rbind()

清洗数据

清洗数据目的只有一个，那就是把citydata的locationId转换成底图数据adcode中包含的。这样就可以链接两张表并进行作图了。主要有两个方面。

1.手动纠正locationId值为（0，-1）地区的locationId值。

这部分数据的产生有几个类别。

• 新疆生产建设兵团的数据。新疆生产建设兵团数据是以师团为单位上报的。由于师团往往分布在好几个省下一级行政区中，我就以师部所在地所属的行政区为准。比如第八师师部所在地石河子市，属于阿勒泰地区，则将这一行的locationId改为阿勒泰地区的行政区划代码654300。

• 无行政区划代码的开发区、高新区等。如重庆的两江新区、南昌的赣江新区等。部分高新区涉及到不同的行政区，这里就不详加考证，划入组成排序第一的城区。

• 因地名不规范而没有提供区划代码的地区。如云南大理州、楚雄州等，也需要进行手动纠正。

要注意的是，这部分数据随着每天新增感染病例在空间上的分散而不断变化，以下代码仅针对2月13日中午丁香园公布的数据。此外，还剩下一些数据因为地区有待明确，或被列为“外地来京”“外地来沪”等类型的病例，无法明确空间上的归属，出于严谨的考虑，这部分数据在制作分布图的时候就不采用了。用city_data=city_data %>% filter(locationId>100)剔除掉。

#
#change adcode manually
{ 
  city_data$locationId[grepl('第四师',city_data$cityName)]=654300
  city_data$locationId[grepl('第六师',city_data$cityName)]=654300
  city_data$locationId[grepl('第七师',city_data$cityName)]=654000
  city_data$locationId[grepl('第八师',city_data$cityName)]=654300
  city_data$locationId[grepl('第九师',city_data$cityName)]=654200
  city_data$locationId[grepl('第十二师',city_data$cityName)]=650100  
  city_data$locationId[grepl('宁东',city_data$cityName)]=640100
  city_data$locationId[grepl('杨凌',city_data$cityName)]=610400
  city_data$locationId[grepl('两江新区',city_data$cityName)]=500105
  city_data$locationId[grepl('高新区',city_data$cityName)]=500106
  city_data$locationId[grepl('保亭',city_data$cityName)]=469029
  city_data$locationId[grepl('昌吉',city_data$cityName)]=652300
  city_data$locationId[grepl('酉阳',city_data$cityName)]=500242
  city_data$locationId[grepl('楚雄',city_data$cityName)]=532300
  city_data$locationId[grepl('文山',city_data$cityName)]=532600
  city_data$locationId[grepl('赣江',city_data$cityName)]=360100
  city_data$locationId[grepl('万盛',city_data$cityName)]=500110
  city_data$locationId[grepl('德宏',city_data$cityName)]=533100
  city_data$locationId[grepl('大理',city_data$cityName)]=532900
  city_data$locationId[grepl('红河',city_data$cityName)]=532500

  city_data=city_data %>% filter(locationId>100)
}

2.把县（市）locationId转换为上级行政区的locationId

这里的县（市）指的是受设区市或自治区管辖的县（市），省直管县级市不在此列。经过转换之后，就可以在下一步的分类汇总中，将县（市）数据汇总到上级行政区中。他们的转换需要两个步骤。

• 怎么识别县（市）？经过上一步骤，所有的地区都有了合法的行政区划代码。那么，代码不在底图中的地区，就是我们需要转换的县（市）了。于是，代码就是row.no=which(!city_data$locationId %in% china_city$adcode)

• 其次，怎么转换。比如我们得到的数据是石家庄市桥西区，其行政区划代码是130104，我们希望把它换成上级行政区也就是石家庄市的代码130100。那我们就要在Amap_adcode这张表里，以桥西区所在行为基准，往上查找距其最近的可以被100整除的（如果用字符型的话，那就是距其最近的代码最后两位是’00’）的代码。查找和替换的核心语句如下:

posi=grep(city_df$cityName[i],Amap_adcode$name)
while(Amap_adcode$adcode[posi]%%100!=0)  posi=posi-1 
city_df$locationId[i]=Amap_adcode$adcode[posi]

该部分完整的代码如下：

#R语言中行和列的待遇大不一样，好像还没有像select和mutate对列操作那样，对行在选择的同时进行操作的函数。所以这里就专门写了一个函数对行进行操作。
row_change=function(dataframe,row.no,fun){
  df_unfilter=dataframe[-row.no,]
  df_filter=dataframe[row.no,] %>%
    fun() %>%
    rbind(df_unfilter)
  return(df_filter)
}
#定义一个对县locationId进行向上替换的函数
upper_class_adcode=function(city_df)
{ 
  for (i in 1:nrow(city_df)){
    if(city_df$locationId[i]>100){
        posi=grep(city_df$cityName[i],Amap_adcode$name)
        while(Amap_adcode$adcode[posi]%%100!=0)  posi=posi-1 
        city_df$locationId[i]=Amap_adcode$adcode[posi]
      }
  }
  return(city_df)
}
#得到县（市）数据所在行
row.no=which(!city_data$locationId %in% china_city$adcode)
#利用以上两个自定义函数进行替换
city_data_new=city_data %>% row_change(row.no,upper_class_adcode) 
#从原始数据中单独提取港澳台数据
GAT=mydata[lapply(mydata$cities,is_empty) %>% unlist,] %>% select(cityName=provinceName,currentConfirmedCount,confirmedCount,suspectedCount,curedCount,deadCount,locationId)
#合并形成完整的城市疫情数据
city_data_new =city_data_new  %>% rbind(GAT)

分类汇总、生成包含数据的sf对象

得到完整的城市疫情数据之后，以locationId为分类变量进行汇总，就可以把同一设区市的县级数据相加，形成设区市的汇总数据。鉴于丁香园上的地名并不规范，我们通过与Amap_code表格进行联结，用规范的地名替换原地名。

city_data_sum=city_data_new %>%
                         select(-cityName)%>%  
                         group_by(locationId) %>% 
                         summarise_all(sum)  %>% 
                        left_join(Amap_adcode,by=c('locationId'='adcode'))

再下一步，将底图的sf对象与汇总后的城市疫情数据链接，就可以生成包含数据的sf对象。

city_data_sf=china_city %>%
             left_join(,city_data_sum,by=c('adcode'='locationId'))%>%
              st_as_sf()
city_data_sf$confirmedCount[is.na(city_data_sf$confirmedCount)]=0              #将未报告患病人数的城市确诊数设为0

作图

这个方面，大神太多，我也是用的其他帖子《2019-nCoV疫情地图动态可视化》。基本上原封不动的抄过来。在我的电脑上，显示中文字体时R语言会崩溃，需要使用showtext包中的showtext_auto这个函数。

mybreaks <- c(0, 1, 10, 50, 100, 500, 1000, 5000, 100000)
mylabels <- c("0", "1-9", "10-49", "50-99", "100-499",
              "500-999", "1000-4999", ">=5000")
city_data_sf=mutate(city_data_sf,conf_level=cut(confirmedCount,breaks=mybreaks,labels=mylabels,include.lowest = T,right=F,ordered_result = T))
ggplot(city_data_sf) + 
  geom_sf(aes(fill = conf_level)) + 
  coord_sf() + 
  scale_fill_brewer(palette = "YlOrRd",) + 
  guides(fill = guide_legend(title = "确诊病例", reverse = T)) + 
  labs(title = "2019-ncov确诊患者空间分布",
       subtitle = Sys.Date(),
       caption = "数据来自丁香园") + 
  theme(
    # 标题
    plot.title = element_text(face = "bold", hjust = 0.5,
                              color = "black"),
    plot.subtitle = element_text(face = "bold", hjust = 0.5, size = 20,
                                 color = "red"),
    plot.caption = element_text(face = "bold", hjust = 1,
                                color = "blue"),
    # 图例
    legend.title = element_text(face = "bold",
                                color = "black"),
    legend.text = element_text(face = "bold",
                               color = "black"),
    legend.background = element_rect(colour = "black"),
    legend.key = element_rect(fill = NA), # 图???????薇???
    legend.position = c(0.85, 0.2),
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    # 绘图面板
    panel.background = element_blank(),
    panel.border = element_rect(color = "black", linetype = "solid", size = 1, fill = NA)
  )

推荐阅读：

真香，朕在看了！猛戳“阅读原文”，填写中国远程办公-调查问卷