walkingsdbgの日記

データ分析技術(機械学習・統計等)についてのメモを書きます。

新型コロナウイルス陽性者数の時間変化を日本地図に描く


はじめに

 データ可視化手法の勉強のため、新型コロナウイルス陽性者数の時間変化を日本地図にヒートマップで描きます。本記事は備忘録的な意味合いが強いので、読みにくい点が多いかと思います。あらかじめご了承ください。

作図方法

 以下の手順に従ってグラフを作成します。

  1. 日本地図のデータを取得する
  2. 新型コロナウイルスの陽性者数のデータを取得する
  3. 日本地図と陽性者数のデータを結合する
  4. 日本地図に陽性者数のヒートマップを載せる
  5. ヒートマップの時間変化をアニメーションで表現する

なお、グラフ作成に使用したパッケージを先に示しておきます。

library(jsonlite)
library(dplyr)
library(sf)
library(animation)
library(ggplot2)
library(lubridate)

日本地図のデータを取得する

 Rのパッケージ{NipponMap}のjpn.shpに日本地図のデータがあります。パッケージのインストール後、インストール先フォルダから以下のように読込みます。*1

d_map <- read_sf(system.file("shapes/jpn.shp", package = "NipponMap")[1],
                 crs = "+proj=longlat +datum=WGS84")

新型コロナウイルスの陽性者数のデータを取得する

 2020年4月末からの陽性者数のデータが以下のサイトにjson形式で公開されています。
https://opendata.corona.go.jp/api/Covid19JapanAll
ダウンロード後、以下のように読込みます。

d_covid_json <- fromJSON("data/Covid19JapanAll.json")

このデータには日次陽性者数が含まれていない(累計陽性者数(npatient)として保存されている)ため、以下のように累計陽性者数の差分を取ることで日次陽性者数を計算します*2

d_covid <- d_covid_json$itemList %>%
  rename(pref = name_jp,
         npatients_acc = npatients) %>%
  mutate(date = as.Date(date),
         npatients_acc = as.numeric(npatients_acc)) %>%
  group_by(pref) %>%
  arrange(date) %>%
  mutate(npatients = c(min(npatients_acc),diff(npatients_acc))) %>%
  ungroup()

さらに、作図用のデータとして、新規陽性者数の週次平均を以下のように計算します。

d_covid <- d_covid %>%
  mutate(date_week = floor_date(date,"week")) %>%
  group_by(date_week,SP_ID,pref) %>%
  summarise(npatients = mean(npatients) + 10^(0)) %>% #後で対数表記を行うので、log(0)を避けるために微小な値を足しておく
  ungroup()

日本地図と陽性者数のデータを結合する

日本地図のデータと陽性者数のデータには都道府県の項目が共通しているので、これをキーにして結合します。ただし、日本地図のデータ内の都道府県はID表記である一方、陽性者数のデータ内の都道府県は漢字表記であるため、陽性者数のデータの都道府県をIDに変換してから結合します。

d_map$SP_ID <- as.numeric(d_map$SP_ID) #陽性者数のデータと結合するため、numeric型に揃えておく
d_covid_map <- d_covid %>%
    mutate(
         SP_ID = case_when(
           pref == "北海道" ~ 1,
           pref == "青森県" ~ 2,
           pref == "岩手県" ~ 3,
           pref == "宮城県" ~ 4,
           pref == "秋田県" ~ 5,
           pref == "山形県" ~ 6,
           pref == "福島県" ~ 7,
           pref == "茨城県" ~ 8,
           pref == "栃木県" ~ 9,
           pref == "群馬県" ~ 10,
           pref == "埼玉県" ~ 11,
           pref == "千葉県" ~ 12,
           pref == "東京都" ~ 13,
           pref == "神奈川県" ~ 14,
           pref == "新潟県" ~ 15,
           pref == "富山県" ~ 16,
           pref == "石川県" ~ 17,
           pref == "福井県" ~ 18,
           pref == "山梨県" ~ 19,
           pref == "長野県" ~ 20,
           pref == "岐阜県" ~ 21,
           pref == "静岡県" ~ 22,
           pref == "愛知県" ~ 23,
           pref == "三重県" ~ 24,
           pref == "滋賀県" ~ 25,
           pref == "京都府" ~ 26,
           pref == "大阪府" ~ 27,
           pref == "兵庫県" ~ 28,
           pref == "奈良県" ~ 29,
           pref == "和歌山県" ~ 30,
           pref == "鳥取県" ~ 31,
           pref == "島根県" ~ 32,
           pref == "岡山県" ~ 33,
           pref == "広島県" ~ 34,
           pref == "山口県" ~ 35,
           pref == "徳島県" ~ 36,
           pref == "香川県" ~ 37,
           pref == "愛媛県" ~ 38,
           pref == "高知県" ~ 39,
           pref == "福岡県" ~ 40,
           pref == "佐賀県" ~ 41,
           pref == "長崎県" ~ 42,
           pref == "熊本県" ~ 43,
           pref == "大分県" ~ 44,
           pref == "宮崎県" ~ 45,
           pref == "鹿児島県" ~ 46,
           pref == "沖縄県" ~ 47)) %>%
  left_join(.,d_map %>% select(SP_ID,geometry),key = "SP_ID")

日本地図に陽性者数のヒートマップを載せる

Rのパッケージ{ggplot2}のgeom_sf()を使えば、日本地図をプロットすることができます。ただし、入力データは予めdata.frameからsfオブジェクト(都道府県と属性情報がセットになった47行のオブジェクト)に変換する必要があります。そこで、パッケージ{sf}をインストールし*3、このパッケージに含まれる関数st_as_sf()を用いてsfオブジェクトに変換後、以下のように日本地図上に陽性者数のヒートマップを描きます。

#地図上の沖縄県の位置を右上に移動する
    d_covid_map$geometry[grep("沖縄県", d_covid_map$pref)] <-
        d_covid_map$geometry[grep("沖縄県", d_covid_map$pref)] + c(5, 17)
#沖縄県の移動先に補助線を引く
    layer_autoline_okinawa <- function(
                    x = c(129, 132.5, 138),
                    xend = c(132.5, 138, 138),
                    y = c(40, 40, 42),
                    yend = c(40, 42, 46),
                    size = ggplot2::.pt / 15
                  ){
                     ggplot2::annotate("segment",
                         x = x,
                         xend = xend,
                         y = y,
                         yend = yend,
                         size = .pt / 15
                  )}

#日本地図上に陽性者数のヒートマップを作成
    plotdate_s <- as.Date('2020-04-26') #適当に日付を設定
    plotdate_e <- plotdate_s + 6
    sf_plot <- st_as_sf(d_covid_map %>% filter(date_week == plotdate_s))
    p <- ggplot(sf_plot)
    p <- p + geom_sf(aes(fill = npatients))
    p <- p + layer_autoline_okinawa()
    p <- p + scale_fill_gradient(low="white", high="#E60012", trans="log10", limits=c(10^(0),10^4), name="陽性者数")
    p <- p + labs(title = paste0("都道府県別新規陽性者数(週次平均)\n",as.character(plotdate_s),"~",as.character(plotdate_e))) 
    p <- p + theme(
      panel.background = element_rect(fill = "midnightblue"),
      panel.grid  = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
    )
    plot(p)

ここで、プログラム冒頭で、日本地図上の沖縄県の位置を南西から北西に変更しています。これは、沖縄県の位置が南西のままだと日本地図のサイズが大きくなり、本州が見づらくなるためです。なお、沖縄県の位置変更については、以下の記事・つぶやきを参考にさせていただきました。
rpubs.com

ヒートマップの時間変化をアニメーションで表現する

前節で作成した日本地図上のヒートマップは、ある特定の時期における陽性者数を表すものです。このヒートマップを時期に応じて変化させるため、以下のような手順でアニメーションを作成します。

  1. 前節のヒートマップを、時期を変化させてN枚作成する。
  2. N枚のヒートマップを、パッケージ{animation}を利用して1つのgifファイルにする

これは以下のプログラムで実現できます。

#GIFアニメーションのオプション設定
ani.options(loop = TRUE, outdir=getwd(), convert = 'program/ImageMagick-7.1.0-4-portable-Q16-x64/convert.exe')

#GIFアニメーションの作成
plotdate_list <- d_covid_map %>% arrange(date_week) %>% distinct(date_week) %>% pull(date_week)
saveGIF({
  for (n in 1:length(plotdate_list)){
    #各週の陽性者数のグラフを描く
    plotdate_s <- plotdate_list[n]
    plotdate_e <- plotdate_s + 6
    sf_plot <- st_as_sf(d_covid_map %>% filter(date_week == plotdate_s))
    p <- ggplot(sf_plot)
    p <- p + geom_sf(aes(fill = npatients))
    p <- p + layer_autoline_okinawa()
    p <- p + scale_fill_gradient(low="white", high="#E60012", trans="log10", limits=c(10^(0),10^4), name="陽性者数")
    p <- p + labs(title = paste0("都道府県別新規陽性者数(週次平均)\n",as.character(plotdate_s),"~",as.character(plotdate_e))) 
    p <- p + theme(
      panel.background = element_rect(fill = "midnightblue"),
      panel.grid  = element_blank(),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
    )
    plot(p)
  }
}, movie.name="covid19_JapanPref.gif", interval=0.25)

このプログラムでは、saveGIF()というgifファイルを作成・保存する関数を使用し、引数として、各週の陽性者数(週次平均)のヒートマップをforループで作成しています。冒頭のani.option()ではgifの詳細設定(gifのループ有無・出力先フォルダ指定・gif作成に必要なexeファイルの読込)を行っています。ここで注意しないといけないのは、ani.option()で指定しているように、gif作成にはImageMagickというソフトが別途必要であることです。
imagemagick.org
このソフトには色々なバージョンがあってややこしいのですが、以下のサイトのアドバイスに従い、portable版をダウンロード・解凍しました。ImageMagickを利用するには、上記プログラムのように、ani.option()にてImageMagickのconvert.exeへのパスを指定すればOKです。
higuma.github.io

実行結果

 これまでの手順を実行すると、以下のようなgifアニメーションが作成できます。期間は2020年5月~2021年7月です。東京等から全国に感染が広がる様子と、緊急事態宣言による感染者数の増減がよく見えます。

f:id:walkingsdbg:20210812000422g:plain
新型コロナウイルス陽性者数の時間変化を表すアニメーション

まとめ

 Rのパッケージ{NipponMap}・{sf}・{animation}とgif作成ソフトImageMagicを用いて、新型コロナウイルスの新規陽性者数(週次平均)の変化を2020年5月~2021年7月の期間で地図上に描きました。ある指標の地域別での時間変化を概観するには便利な方法だと感じたので、仕事でも機会があれば使ってみようと思います。

補足

本記事のプログラムは以下にまとめております。詳細に興味がある方はこちらをご覧ください。
github.com
 また、本記事の作成にあたって、以下の記事を参考にさせていただきました。
dinov.sakura.ne.jp
totech.hateblo.jp

*1:インストール後、library(NipponMap)のようにパッケージを読込む必要は無いことに注意してください。

*2:データの中には、なぜか累計陽性者数が時間とともに減少する(=日次陽性者数を計算すると負になる)部分があるので、その日の陽性者数は0で置き換えています。

*3:パッケージ{sf}を読込む前に、最新版の{Rcpp}をインストールすることを推奨します。私の環境では1.0.7より古いバージョンのRcppを使うと、sfの読込時にエラーが出ました。