Data and methods

Libraries and my function

> library(readr)
> library(readxl)
> library(lubridate)
> library(dplyr)
> library(tidyr)
> library(magrittr)
> library(sp)
> library(scanstatistics)
> library(classInt)
> library(RColorBrewer)
> library(sf)
> library(jpndistrict)
> library(animation)
> 
> findc2 <- function(dat,bk,col){
+   rtn <- rep(col[1],length(dat))
+   for(i in 1:length(dat)){
+     for(j in 2:(length(bk)-1)){
+       if(bk[j]<=dat[i]){
+         rtn[i] <- col[j]
+       } 
+     }
+   }
+   rtn
+ }

Preparation

library(dplyr)
library(tidyr)
library(magrittr)
library(estatapi)
appId <- "**********************"
StatsList <- estat_getStatsList(appId = appId, searchWord = "", statsCode = "00200521")
grep("昼夜間", StatsList$TITLE)
StatsList$`@id`[grep("昼夜間", StatsList$TITLE)[3]]
statsDataId <- "0003179188"
meta_info <- estat_getMetaInfo(appId = appId, statsDataId = statsDataId)
df <- estat_getStatsData(appId = appId, statsDataId = statsDataId)
dn.tokyo <- left_join(
  df %>% filter(`常住地又は従業地・通学地_2015` == "常住地による人口（夜間人口）") %>% filter(area_code >13000, area_code <14000) %>% select(area_code, pop = value),
  df %>% filter(`常住地又は従業地・通学地_2015` == "（別掲）昼夜間人口比率") %>% filter(area_code >13000, area_code <14000)%>% select(area_code, dnratio = value),
  by="area_code"
)
dn.tokyo <- dn.tokyo %>% mutate(ccode = area_code, dnratio = dnratio/100) %>% select(ccode, pop, dnratio)
write.csv(dn.tokyo, "data/day_night_tokyo.csv", row.names = FALSE)

data <- read_excel("data/suzuki/tokyo_municipality_infected_0502.xlsx")
ld <- ncol(data)-1
colnames(data) <- c("cname", as.character(ymd("2020-3-31")+0:(ld-1)))
pop <- read_csv("data/suzuki/popmaster.csv", locale = locale(encoding = "SHIFT-JIS"))
dn.tokyo <- read_csv("data/day_night_tokyo.csv")
pop <- left_join(pop, dn.tokyo, by="ccode") %>% select(-pop)
data <- left_join(pop, data, by="cname")
jugyo <- read_excel("data/suzuki/ssdse.xlsx") %>% select(cname = 市区町村, jugyo = '従業者数/1000')
jugyo <- jugyo[-63,]
jugyo  %<>% mutate(cname=substr(cname,1, str_length(cname)-1))
data %<>% left_join(jugyo, by="cname")
data %<>% mutate(pop1 = popk, pop2 = (1+dnratio)/2*popk, pop3 = (popk + jugyo)/2) %>% select(ccode,pop1, pop2, pop3,`2020-03-31`:`2020-05-02`)
write.csv(data, "tokyo20200502.csv", row.names = FALSE)

Load data and data cleaning

> data <- read_csv("tokyo20200502.csv")
> ld <- ncol(data)-4
> colnames(data) <- c("ccode","pop1", "pop2", "pop3", as.character(ymd("2020-3-31")+0:(ld-1)))
> 
> data.counts <- data[,c(1,5:(ncol(data)))]
> data.counts.long <- gather(data.counts, key = "date", value = "count", -ccode)
> 
> data.counts.long <- left_join(data.counts.long, data %>% select(ccode, pop2), by="ccode")
> data.counts.long <- data.counts.long %>% mutate(rate = count/pop2)
> data.counts.long$rate[is.nan(data.counts.long$rate)] <- 0
> 
> nclass <- 6
> rdpu6 <- brewer.pal(nclass,"RdPu")
> cls <- classIntervals(as.numeric(data.counts.long$count), nclass, style = "jenks")
> cls.rate <- classIntervals(as.numeric(data.counts.long$rate), nclass, style = "jenks")
> 
> ip <- 13
> sf_pref <- jpn_pref(ip, district = TRUE) 
> sf_pref <- filter(sf_pref, city_code<13360)
> 
> leg.plc <- "bottomleft"
> map.data <- data.frame(ccode = as.integer(sf_pref$city_code))
> dd <- as.character(ymd("2020-3-31")+0:(ld-1))
> 
> ############ scanstatistics
> 
> #data.counts.long <- data.counts.long %>% select(-rate)
> counts <- data.counts.long %>% 
+   df_to_matrix(time_col = "date", location_col = "ccode", value_col = "count")
> 
> population <- data.counts.long %>%
+   df_to_matrix(time_col = "date", location_col = "ccode", value_col = "pop2")
> 
> points <- read_csv("data/points.csv",locale = locale(encoding = "SHIFT-JIS"))
> points <- points %>% mutate(ccode = as.integer(jiscode))
> points <- left_join(left_join(map.data, data %>% filter(pop2 != 0),by="ccode") %>% select(ccode), points, by="ccode") 
> points <- points %>% select(ccode, lon, lat)
> zones <- points %>% select(lon, lat) %>% 
+   as.matrix %>% 
+   spDists(x = ., y = ., longlat = TRUE) %>% 
+   dist_to_knn(k = 10) %>% 
+   knn_zones
saveGIF({
    for(id in 1:ld){
    par(mfrow=c(3,1))
    data.color <- left_join(map.data, data.counts.long %>% filter(date == dd[id])%>% select(ccode,count), by="ccode")
    idxz <- c(which(data.color$count == 0))
    col.this <- findc2(data.color$count, cls$brks, rdpu6)
    col.this[idxz] <- "#FFFFFF"
    text.lg <- paste("(",cls$brks[-(nclass+1)],",",cls$brks[-1],"]",sep="")
    text.lg <- c("[0 or NA]", text.lg)
    main <- paste("Number of infected for covid19 (daily) (", dd[id], ")", sep="")
    sf_pref %>% st_geometry() %>% plot(col = col.this, main=main)
    legend(leg.plc, legend=text.lg, fill = c("#FFFFFF",rdpu6), cex=1)
    
    
    data.color.rate <- left_join(map.data, data.counts.long %>% filter(date == dd[id])%>% select(ccode,rate), by="ccode")
    #idxz <- c(which(data.color$count == 0))
    col.this.rate <- findc2(data.color.rate$rate, cls.rate$brks, rdpu6)
    col.this.rate[idxz] <- "#FFFFFF"
    text.lg.rate <- paste("(",round(cls.rate$brks[-(nclass+1)],2),", ",round(cls.rate$brks[-1],2),"]",sep="")
    text.lg.rate <- c("[0 or NA]", text.lg.rate)
    main <- paste("Rate of infected (/1k pop) for covid19 (daily) (", dd[id], ")", sep="")
    sf_pref %>% st_geometry() %>% plot(col = col.this.rate, main=main)
    legend(leg.plc, legend=text.lg.rate, fill = c("#FFFFFF",rdpu6), cex=1)
    
    set.seed(1)
    # remove pop = 0
    idxpnz <- which(population[1,]!=0)
    pb_poisson_result <- scan_pb_poisson(counts = counts[id,idxpnz], 
                                         zones = zones, 
                                         population = population[id,idxpnz],
                                         n_mcsim = 999,
                                         max_only = FALSE)
    
    #points[pb_poisson_result$MLC$locations,]$ccode
    idx.hs <- map.data %>% mutate(idx = 1:nrow(map.data)) %>% filter(ccode %in% points[pb_poisson_result$MLC$locations,]$ccode) %>% select(idx)
    idx.hs <- as.integer(idx.hs$idx)
    
    col.this.hs <- rep("#FFFFFF", nrow(map.data))
    col.this.hs[idx.hs] <- "#FF0000"
    text.lg <- c("not MLC","MLC")
    main <- paste("Hot spot of infected for covid19 (daily) (", dd[id], ")", sep="")
    sf_pref %>% st_geometry() %>% plot(col = col.this.hs, main=main)
    legend(leg.plc, legend=text.lg, fill = c("#FFFFFF","#FF0000"), cex=1)
  }
},movie.name = "20200502.gif", img.name = "Tokyo0331to0502",ani.width = 550,ani.height=1200)

References