R语言与新型冠状病毒(02)地图

来自 Peng Zhao | February 1, 2020

失眠,索性起来写第二篇。

书接上回

在上一节里,为了避免 Windows 用户遭遇中文乱码,我们将系统语言设置为中文:

Sys.setlocale('LC_CTYPE', 'Chinese')
## [1] "Chinese (Simplified)_China.936"

我们安装和加载了 “ncovr” 包(请重新安装一次,因为我刚刚把它升级到 0.0.2 版本了),用它来获取了疫情数据:

require(remotes)
install_github("pzhaonet/ncovr")
require(ncovr)
ncov <- get_ncov()

我们还安装了 “tidyverse” 包,但是并未调用它。今天,我们就要启用其中的两个成员:“dplyr”和 “magrittr”。

require(dplyr)
require(magrittr)

相信大家最近这些天,对下面这张地图都很熟悉了吧:

knitr::include_graphics('https://openr.pzhao.org/img/map-ncov-netease.png')

本节将以 ncov$area 数据集为基础,绘制这张地图。

准备工作

绘制地图有很多种方法(详见《学 R》第十一章:“地图”)。这里,我们选择最酷的一种,那就是交互地图,支持缩放和弹出信息的那一种。绘图所用的包是 “leaflet” 和 “leafletCN”,在上面重装 “ncovr” 包时已经自动安装了,现在只需加载一下:

require(leaflet)
require(leafletCN)

准备工作完成。

试画地图

先浏览一下数据:

View(ncov$area)

在 `provinceName 一列里,可以发现有个“待明确地区”,这显然不是省名称,我们先把这条记录删掉:

ncov$area <- ncov$area[ncov$area$provinceName != "待明确地区",]

可能会出乎你的意料,这最酷的偏偏画起来是最简单的。只需一条指令,就可以绘制全国新型冠状病毒的确诊病例地图:

# 试画地图
geojsonMap(dat = ncov$area, 
           mapName = "china", 
           palette="Reds",
           namevar = ~provinceName,
           valuevar = ~confirmedCount,
           legendTitle = "确诊人数")

点击左上角的缩放键,就可以 zoom in/out 了。如果是触屏,可以用两个手指放大缩小。点击地图上任一位置,就会弹出对应信息(当然,如果你是在微信公众号读到了本文,是无法进行交互的。请点击文末左下角的“阅读原文”,移步网页版)。

酷吧?!

改进地图

酷是酷,这幅地图有两个问题:

  1. 代表确诊病例颜色的深浅是线性的。湖北省颜色最深,而其他省之间的差异不太明显。这是由于其他省跟湖北省差别很大而导致的。
  2. 点击各省,弹出的只有个省名称。要是能显示确诊病例数量就更好了。

对于问题 1,我们用 cut() 函数对确诊病例数进行分组;对于问题 2,只需用 geojsonMap() 函数的 popup 参数来设定:

# 修改后的地图
count_cut <- c(0, 9, 99, 999, Inf)
ncov$area$confirmed_level <- cut(
  ncov$area$confirmedCount, 
  count_cut, 
  labels = c('< 10', '10-99', '100-999', '>999'))
geojsonMap(dat = ncov$area, 
           mapName = "china", 
           colorMethod = "factor",
           palette="Reds",
           namevar = ~provinceName,
           valuevar = ~confirmed_level,
           popup =  paste(
             ncov$area$provinceName, 
             ncov$area$confirmedCount),
           legendTitle = "确诊人数")

目标地图绘制完成!

作业:仿照上面的方法,绘制全国治愈病例地图和死亡病例地图。

超额完成

目标光完成是不够的;让我们挑战一下自我,除了完成,还要超额。

上面的地图里,虽然按数量级将各省分了组,但问题是分在同一组里的省份看不出区别。这该怎么办?

大学老师教过我们,当跨越两三个数量级并且要区分低值时,可以考虑对数坐标。

然而,查阅 “leafletCN” 包的帮助文档,却发现 geojsonMap() 函数并不支持用对数坐标来显示颜色。怎么办呢?

我抓破头皮,想了个办法,就是把 geojsonMap() 函数给黑了,重新写了个函数,叫做 geojsonMap_legendless(),收在了 “ncovr” 包里(0.0.2 版)。对,这是我刚写的,所以在本文开头请你重装了“ncovr” 包。现在,可以使用这个新函数,配合 “leaflet” 包的 addLegend() 函数,就可以绘制对数坐标图:

# 超额绘制地图
ncov$area$confirmed_log <- log10(ncov$area$confirmedCount)
ncov$area$confirmed_log[ncov$area$confirmedCount == 0] <- NA
geojsonMap_legendless(
  dat = ncov$area, 
  mapName = "china",
  palette = "Reds", 
  namevar = ~provinceName,
  valuevar = ~confirmed_log,
  popup =  paste(
    ncov$area$provinceName, 
    ncov$area$confirmedCount)) %>%
  leaflet::addLegend(
    "bottomright",
    bins = 5,
    pal = leaflet::colorNumeric(
      palette = "Reds",
      domain = ncov$area$confirmed_log
    ),
    values = ncov$area$confirmed_log,
    title = "确诊人数",
    labFormat = leaflet::labelFormat(
      transform = function(x) 10 ^ x),
    opacity = 1)

对数坐标的缺点在于 0 值取对数后成为负无穷大,难以展示。这里把 0 用灰色表示,跟缺少数据(NA)的区域颜色相同。

城市地图

最后,我们来做个更难的实验:

有个朋友前些天从国外回到了浙江省温州市,发现温州有了个新名字,叫“湖北省温州市”:疫情按严重程度排名的话,温州是湖北之外疫情最严重的一个城市。然而在全国地图上看不出来这个信息。

我也想看看我的家乡河南省,各市之间有什么差别。这就需要绘制以市为单元的地图了。

ncov$area 数据集为基础绘制城市地图,有 2 个难点:

  1. 城市疫情信息藏在 cities 一列,形式是 list ,并不是绘图需要的 data.frame (详见《学R》第八章)。
  2. 然而,直辖市却存在 provinceName 一列。这就需要提取出来,跟其他城市合并。

解决这两个难点的方法是这样的:

# 城市地图
ncov_city <- dplyr::bind_rows(ncov$area$cities) %>% dplyr::select(1:5)
ncov_area <- ncov$area[, 2:6]
names(ncov_area) <- names(ncov_city)
ncov_cities <- rbind(ncov_city, ncov_area)
cities <- regionNames(mapName = 'city')
ncov_cities <- ncov_cities[ncov_cities$cityName %in% cities, ]
ncov_cities$confirmed_log <- log10(ncov_cities$confirmedCount)
ncov_cities$confirmed_log[ncov_cities$confirmedCount == 0] <- NA
geojsonMap_legendless(
  dat = as.data.frame(ncov_cities), 
  mapName = "city",
  palette = "Reds", 
  namevar = ~cityName,
  valuevar = ~confirmed_log,
  popup =  paste(ncov_cities$cityName, 
                 ncov_cities$confirmedCount)) %>%
  leaflet::addLegend(
    "bottomright", 
    bins = 5,
    pal = leaflet::colorNumeric(
      palette = "Reds",
      domain = ncov_cities$confirmed_log),
    values = ncov_cities$confirmed_log,
    title = "确诊人数",
    labFormat = leaflet::labelFormat(
      transform = function(x) 10 ^ x),
    opacity = 1)