2022年以来的中国本土新型冠状病毒肺炎新增病例情况

作者：

日期： 2022-05-10 时间： 16:56

主题背景

2022年以来，中国31个省（自治区、直辖市）和新疆生产建设兵团报告的新型冠状病毒COVID-19肺炎疫情不容乐观。每天早餐看到的疫情新闻、每天工作群里的各种时空伴随摸排、时紧时松的核酸检测，让每一个老百姓切身感受到了不容乐观的疫情。尤其是3、4月的吉林和上海，更是让全中国忧心忡忡。时间来到4月上旬，我突然希望看到今年以来的全国各省疫情的变化情况，如果能像最近学习的山脊图那样呈现出来就再好不过了。可去哪里看呢？上网找不到（也没仔细找），自己想想办法。

数据获取

大多数头部网站都是站在访问者的角度进行的数据展现，最常见的就是获取访问者所在的省份信息作为默认条件提供当日的疫情信息，同时肯定会提供全国的总体情况。每个省的情况也可通过额外的查询条件获取，但同样是当日的信息，很少有历史数据的展示。结果便是没有办法通过爬虫工具抓取经过整理的规则数据，有历史数据的公开资源在哪里？找来找去发现能够满足要求的只有国内疫情数据的源头——中华人民共和国国家卫生健康委员会。但问题是，这里的数据是新闻通报格式，也就是纯文字描述的形式，它长这个样子：

1月1日0—24时，31个省（自治区、直辖市）和新疆生产建设兵团报告新增确诊病例191例。其中境外输入病例60例（上海18例，广东10例，福建9例，广西8例，天津5例，四川4例，浙江2例，山东2例，辽宁1例，云南1例），含6例由无症状感染者转为确诊病例（广东2例，四川2例，浙江1例，山东1例）；本土病例131例（陕西123例，其中西安市122例、延安市1例；浙江7例，均在宁波市；河南1例，在洛阳市），含1例由无症状感染者转为确诊病例（在河南）。无新增死亡病例。新增疑似病例1例，为境外输入病例（在上海）。

于是…从4月8日开始，一直到5月8日，我花了整整一个月的时间，每天晚上20分钟，从卫健委官方网站上一个页面一个页面的把2022.1.1~2022.5.8的疫情通报里的一个数据一个数据给抄录下来，记录到一张Excel表格中。这个采集数据的过程真痛苦。采集数据的方法，我应该是不得要领，“愚公移山”。

数据处理

抄录的数据有以下3个特点：

  1. 每天上榜的城市不尽相同，因为每日有新增确诊病例才会上榜
  2. 每天的数据包括境外输入、本土新增两种类型
  3. 少了2个省份数据：西藏、宁夏。（这两个省的疫情防控做得真是漂亮）

由于山脊图一般是多组在时间/空间维度上连续的观测值连线在纵向上的叠加排布。比如2017年Henrik Lindberg绘制的那张流行的一天中运动和休闲的高峰时间。

运动与休闲

因此，原始抄录数据还需要进行处理，包括：

  1. 去除境外输入数据。出于这部分人群入境后有比较好的管控闭环，且数量相对少的原因。在抄录数据的过程中，我感觉其与本土新增病例数据变化也不正相关，对各省疫情变化不构成主要影响因素，其规律后续再单独分析。这个操作使用Excel自带的筛选功能即可完成。
  2. 复制筛选后的数据。将剩下的本土新增病例数据复制到一个新的Sheet，作为后续加工的原始数据集。这个操作直接使用Excel复制-新工作表-粘贴即可完成。

  3. 按城市补齐所有日期的数据。从2022.1.1~2022.5.8，一个城市一个城市的来，当天没有本土新增病例补零即可。这个处理也在Excel中完成，因会改变原有的数据结构，操作稍微复杂一点。

    • 在Excel中新建一张Sheet

    • 将第一列作为日期，使用填充功能按日期从2022.1.1开始，填充日期序列，到2022.5.8截止，共128行。

    • 将第二列作为省份，从原Sheet中的省份拷贝后分别填入，一个时间周期（128天）内的省份是同一个。

    • 按上面的两个步骤，将原Sheet中的30个省份和兵团同样的方式扩展数据，得到128*30=3840行数据。

    • 将第三列作为本土新增病例数据，编写公式将第一列和第二列作为条件到原Sheet中匹配查找对应的新增病例数，未匹配到代表当天无该城市新增病例，取0。输入以下公式后记得Ctrl+Shift+Enter
      =IFERROR(VLOOKUP([@Date]&[@Province],IF({1,0},表4[Date]&表4[Province],表4[New Case]),2,0),0)

    • 新增第四列计算本土7天移动平均日增确诊病例数，希望可以让可视化后的曲线更平滑一些。在每个城市的2022.1.7对应的7 Days Smooth单元格输入以下公式，并向下填充至2022.5.8
      =AVERAGE(C2:C8)

    • 将第四列每个城市2022.1.1~2022.1.6对应的7 Days Smooth单元格填充0。

最后处理完成的数据表，我把它共享在这里COVID-19 Data。

数据可视化

这个部分使用R，过程中有一些细节的处理，我在代码说进行注释。

  1. 导入XLSX数据
library("xlsx")
covid <- read.xlsx("COVID-19_2022.xlsx", 3, stringsAsFactors = FALSE, encoding = "UTF-8")
  1. 绘图
# 由于7天移动平均数存在大量小于1的数值，而上海和吉林两个地区的数据峰值接近3000，
# 在实际绘图过程中由于这两者数值差异巨大，同一坐标参考系下无法让其他省份被明示，
# 通过取log对数来缩小数值间的差异是可行的办法，但小于1取对数会出现负值，明显不合适。
# 综上原因，舍弃对7天移动平均数进行绘图，仍选择本土新增绝对值。
# 
# 求数据集中的本土新增最大值
max_all <- max(covid[, 3])
# 定义绘图重叠值，是指与没有重叠相比，一个图将覆盖另一个图的相对量。
max_overlap <- 10
# 获取数据集中的省份数量，此例子为30
provnum <- length(unique(covid$Province))
# 定义绘图的初始高度，比如需要绘制30个省份，那么可以考虑将整个绘图区域分成
# 40等分，一个省份的绘图高度为10等分
initheight <- 1 / (provnum + max_overlap)
# 定义30个省份的绘制颜色
ncolors <- 30
pal <- colorRampPalette(c("#ff35de", "#2be826", "#d6e826", "#ff0a0b"))
pcolor <- sample(pal(ncolors))
# 设定画布参数
par(mar = c(0, 0, 0, 0), bg = "black")
# 绘制空白区域
plot(0:1, 0:1, type = "n", xlab = "", ylab = "", axes = FALSE, asp = 1)
# 按本土日新增病例累积从小到大排序，上海排最后
for (i in (1:provnum)) {
  j <- (i - 1) * 128 + 1
  k <- j + 127
  sumscases[i] <- sum(covid[c(j:k),3])
}
provorder <- order(sumscases, decreasing = FALSE)
# 绘制山脊图
for (i in (provnum:1)) {
  # i从大到小，先绘制上方，视觉上是最后面的山脊
  par(new = TRUE, plt = c(0, 1, initheight * (i - 1), initheight * (i - 1) + initheight * max_overlap))
  # 根据本土日新增病例累积数排序从大到小，根据排序对应的索引值，找到该省份在原始数据中的开始行和结束行
  j <- (provorder[i] - 1) * 128 + 1
  k <- j + 127
  # 设置当前的绘图空间，尤其是xlim和ylim
  plot(log(ifelse(covid[c(j:k), 3]==0,1,covid[c(j:k), 3])), xlim = c(-1, 130), ylim = c(-2, log(max_all)+1), xlab = "", ylab = "", type = "n", axes = FALSE)
  if (i == provnum) {
    # 绘制右上角坐标轴，纵坐标刻度等分，对应实际值是以指数级近似增加。
    axis(4, at=round(seq(0, 8.61, length.out = 10), 2), labels = paste(c(0, 0.003, 0.01, 0.02, 0.05, 0.1, 0.4, 0.8, 2, 5.49), "K", sep = ""), 
         line = -3, lwd = 0, lwd.ticks = .5, cex.axis = .7, col = "#f0f0f0", col.axis = "#f0f0f0", las = 2)
  }
  if (i == 1) {
    # 绘制横坐标轴
    axis(1, at=c(1, 32, 60, 91, 121, 128), labels = c("2022.1.1", "2022.2.1", "2022.3.1", "2022.4.1", "2022.5.1", "2022.5.8"), line = -1.8, lwd = .2, lwd.ticks = .2, cex.axis = .7, col = "#f0f0f0", col.axis = "#f0f0f0", las = 0)
  }
  # 绘制多边形，实际新增为0时，考虑取对数将计算异常，因此判断日新增为0时，输出1后在取log对数，结果仍得到0
  # 实际新增为1时，取对数后等于0，因此绘图最终会有偏差。
  polygon(c(1:128,128:1),c(log(ifelse(covid[c(j:k), 3]==0,1,covid[c(j:k), 3])), rep(0, 128)), col = pcolor[i], border = NA)
  # 在多边形的上侧绘制灰色线条
  lines(log(ifelse(covid[c(j: k), 3]==0,1,covid[c(j:k), 3])), lwd = .8, col = "#f0f0f0")
  # 在图形左侧对应输出省份名
  text(1.7, 0, covid[j, 2], pos = 2, cex = .7, col = "#f0f0f0")
}

图像加工

按照惯例将图形作为PDF文件输出，然后使用Adobe Illustrator进行画板尺寸、文字说明等相关编辑后，输出JPEG文件。

日新增确诊情况最终可视化效果

小结

从接触山脊图，到最终自己成图，前后花了1个月的时间。收集和整理数据，尤其是收集数据是数据可视化最花时间和精力的环节。数据是客观世界的真实反映，山脊图是反映时间轴上事物发生的变化，其背后的数据不但需要从客观世界的海量信息中进行提取、核实，还需要经过较长时间的积累。但这是一种很直观的反映趋势变化，以及不同类别变化趋势之间的区别，我个人认为是很棒的图表类型。

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注