「都道府県別新規感染報告数の推移を描くコード」に触ってみる

中澤先生の「鐵人三国誌」【第685回】https://minato.sip21c.org/im3r/20210319.htmlでも紹介されている
「都道府県別新規感染報告数の推移を描くコード」を使ってみようと思い，ついでに色々書き足してみたりした。
描画対象とする都道府県を都道府県コードのベクトルで指定したり，片対数グラフの他に普通メモリでも描けるようにしたりが大きな変更である。そのほか，プログラムを関数にしてデータとプログラムをなるべく分離したので，保守しやすくなったと思うがどうだろうか。
呼出方は，以下のようである。描画対象やパラメータを変えて描画するときに，1行だけ入力すればよいので便利かな？

draw.covid19(c(13, 27, 28, 4), event.data)
draw.covid19(11:14, event.data, log="")

 

# Plotting new confirmed cases for Japan against date in semilog scale.
# November 24, 2020: Made
# December 22, 2020: Change color palette
# January 8, 2021: 3 prefs were added
# January 17, 2021: New action date has been added
# February 8, 2021: Emergency is off at Tochigi pref
# February 28, 2021: Emergency is off at several prefs
# March 23, 2021: rewrite easy to use. by 'r-de-r'

draw.covid19 <- function(prefs, event.data,
                         log="y", # plot()
                         bty="n", # legend()
                         lty="dotted", lwd=0.5, col="gray", # segments()
                         lwd2=NULL, lty2=NULL, # lines()
                         pos=4 # text()
                         ) {

    if (!require(COVID19)) {
        install.packages("COVID19", dep=TRUE)
        library(COVID19)
    }

    # 都道府県名（都道府県コード順）
    dict.en <-  c("Hokkaido", "Aomori", "Iwate", "Miyagi", "Akita",
                  "Yamagata", "Fukushima", "Ibaraki", "Tochigi", "Gunma",
                  "Saitama", "Chiba", "Tokyo", "Kanagawa", "Niigata",
                  "Toyama", "Ishikawa", "Fukui", "Yamanashi", "Nagano",
                  "Gifu", "Shizuoka", "Aichi", "Mie", "Shiga",
                  "Kyoto", "Osaka", "Hyogo", "Nara", "Wakayama",
                  "Tottori", "Shimane", "Okayama", "Hiroshima", "Yamaguchi",
                  "Tokushima", "Kagawa", "Ehime", "Kochi", "Fukuoka",
                  "Saga", "Nagasaki", "Kumamoto", "Oita", "Miyazaki",
                  "Kagoshima", "Okinawa")
    dict.jpn <- c("北海道", "青森", "岩手", "宮城",   "秋田",
                  "山形",   "福島", "茨城", "栃木",   "群馬",
                  "埼玉",   "千葉", "東京", "神奈川", "新潟",
                  "富山",   "石川", "福井", "山梨",   "長野",
                  "岐阜",   "静岡", "愛知", "三重",   "滋賀",
                  "京都",   "大阪", "兵庫", "奈良",   "和歌山",
                  "鳥取",   "島根", "岡山", "広島",   "山口",
                  "徳島",   "香川", "愛媛", "高知",   "福岡",
                  "佐賀",   "長崎", "熊本", "大分",   "宮崎",
                  "鹿児島", "沖縄")

    # ラベル，タイトル
    xlab <- "時間経過（年）"
    ylab <- "新規確定患者報告数"
    main <- "COVID-19の都道府県別新規確定患者報告数の推移
               [Ref.] Guidotti E, Ardia D (2020) COVID-19 Data Hub.
              Working paper https://doi.org/10.13140/RG.2.2.11649.81763"

    # 描画範囲（横軸）
    from <- sprintf("2020-%02d-01", 1:12)
    to   <- sprintf("2021-%02d-01", 1:3)

    # 発生事象描画関数
    event.draw <- function(x) {
        factor <- function(x) if (log == "") log(x)*350 - 250 else x
        date.type = as.Date(x[1,1])
        segments(date.type, 1, date.type, factor(8000), lty=lty, lwd=lwd, col=col)
        text(date.type, factor(x[1,3]), x[1,2],
            cex=0.8, offset=0, pos=pos, xpd=TRUE)
    }

    options(scipen=7)

    # RStudio/Mac OS の場合に文字化け回避のために必要
    par(family= "Hiragino Mincho Pro W3") # ほかに "HiraKakuProN-W3", "Hiragino Kaku Gothic Pro W3"

    # 折線の太さ・線種の定義
    np <- length(prefs)
    if (is.null(lwd2) || length(lwd2) != np) {
        lwd2 <- rep(0.5, np)
    }
    if (is.null(lty2) || length(lty2) != np) {
        lty2 <- rep("solid", np)
    }

    # 都道府県名取り出し
    prefsJ <- dict.jpn[prefs]
    prefs  <- dict.en [prefs]

    # cols <- rainbow(length(prefs))
    cols <- 1:length(prefs)
    palette("Okabe-Ito")
    y <- covid19(country="JPN", level=2, start="2020-01-01", end=Sys.Date())
    y <- subset(y, date < max(date))
    y$prefs <- as.factor(y$administrative_area_level_2)
    cfd <- c(y$confirmed[1], diff(y$confirmed))
    cfd <- ifelse(cfd < 1, NA, cfd)
    plot(cfd ~ as.Date(date, "%d/%m/%y"), data=y,
        xlim=c(min(y$date), max(y$date) + 20), ylim=c(1, ifelse(log=="", 2800, 10000)),
        axes=FALSE, log=log, type="n", frame=FALSE,
        xlab=xlab, ylab=ylab, main=main,
        sub=sprintf("%s まで", max(y$date)))
    axis(1, c(as.Date(from), as.Date(to)), c(from, to))
    if (log == "y") {
        axis(2, 10^(0:4))
    } else {
        axis(2)
    }
    for (i in 1:nrow(event.data)) {
        event.draw(event.data[i,])
    }
    for (i in levels(y$prefs)) {
        dd <- subset(y, prefs == i)
        dd$confirmedpp <- c(dd$confirmed[1], diff(dd$confirmed))
        dd$confirmedpp <- ifelse(dd$confirmedpp < 0, 0, dd$confirmedpp)
        if (i %in% prefs) {
            index <- grep(i, prefs)
            lines(dd$date, dd$confirmedpp, col=cols[index],
                lwd=lwd2[index], lty=lty2[index])
            text(dd$date[length(dd$date)]+1,
                max(dd$confirmedpp[length(dd$confirmedpp)], 1),
                prefsJ[index], col=cols[index], pos=4, xpd=TRUE)
        }
    }
    legend("topleft", col=cols, lwd=2, legend=prefsJ, bty=bty)
}

# 発生事象データ（随時追加）
event.data = data.frame(
    date =  c("2020-04-08", "2020-05-14", "2020-06-19", "2020-07-22", "2020-10-01",
              "2021-01-07", "2021-01-14", "2021-02-08", "2021-02-28"),
    event = c("緊急事態宣言発出", "39県解除", "東京解除", "GoToキャンペーン",
              "GoTo東京追加，\nGoToイート開始", "一都三県\n緊急事態宣言発出",
              "二府五県\n緊急事態宣言追加", "栃木県\n緊急事態宣言解除",
              "二府四県\n緊急事態宣言解除"),
    pos =   c(5000, 3000, 1700, 5000, 3000,
              5000, 3000, 1700, 1000))

# 描画都道府県は draw.covid19() の第1パラメータで c(13, 27, 4) のように都道府県コードで指定する
# 折線グラフの詳細指定（省略可）
lwd2 = rep(0.5, 4); lwd2[c(1, 4)] = 2 # 太さ
lty2 = rep("solid", 4)          # 線種

draw.covid19(c(13, 27, 28, 4), event.data, lwd2=lwd2, lty2=lty2)
draw.covid19(c(13, 27, 28, 4), event.data, log="", lwd2=lwd2, lty2=lty2)