State New Cases COVID-19 3D Visualization in R

covidcases3D.R

library(rayrender)
library(sf)
library(dplyr)
library(ggplot2)
library(lubridate)
library(magick)

#Extract state data
us_states = spData::us_states

#Transform to a different coordinate system
shp_utm = sf::st_transform(us_states, 2163)

#Load NYTimes COVID19 data
coviddata = read.csv("us-states.csv")

#Calculate daily new case count
coviddata %>%
  filter(ymd(date) > ymd("2020-01-31")) %>% 
  group_by(state) %>%
  mutate(newcases = cases-lag(cases))  %>% 
  mutate(newcases = ifelse(is.na(newcases), 0, newcases)) %>%
  ungroup() ->
covid_case_data

#Generate 3D plots (takes a few hours)
tempday = lubridate::ymd("2020-02-01")
finalday = lubridate::today()
counter = 1
while(tempday != finalday) {
  covid_case_data %>%
    filter(date == as.character(tempday)) %>%
    mutate(fips = as.character(fips)) %>%
    filter(!is.na(state)) %>%
    filter(!is.na(cases)) ->
  yesterday
  
  shp_utm %>%
    left_join(yesterday, by = c("NAME"="state")) %>%
    mutate(newcases = ifelse(is.na(newcases),-121000,newcases)) %>%
    mutate(bot = -122000) ->
  newdata
  
  xz_rect(y=-0.1,zwidth=500,xwidth=500,
          material = glossy(color="grey5",gloss=0.7,reflectance=0.1)) %>%
    add_object(extruded_polygon(newdata, center=TRUE,z=-5,scale=c(1/100000,1,1/100000), data_column_top = "newcases",
                                scale_data = 1/max(covid_case_data$newcases,na.rm=TRUE)*15, data_column_bottom = "bot",
                                material=glossy(color="#754315",gradient_color="darkred",
                                                gradient_point_start = c(0,0,0),gradient_point_end = c(0,8,0),
                                                reflectance=0,gloss = 0.3))) %>%
    add_object(sphere(y=200,radius=80,material=light())) %>%
    render_scene(parallel=TRUE,lookfrom = c(0,50,-40),lookat=c(0,-5,0),samples=400,
                 filename=glue::glue("covid{counter}"),
                 fov=0,clamp_value=10,ortho_dimensions = c(50,50),width=1000,height=1000)
  counter = counter + 1
  tempday = tempday + lubridate::duration("1 day")
}


#Create inset line graph of total cases
tempday = lubridate::ymd("2020-02-01")
finalday = lubridate::today()
counter = 1
while(tempday != finalday) {
  covid_case_data %>%
    group_by(date) %>%
    summarise(totalnewcases = sum(newcases)) %>% 
    ggplot() + 
    geom_line(aes(x=ymd(date),y=totalnewcases),color="white",size=0.25) +
    geom_vline(xintercept=ymd(tempday),color="red",size=0.25) +
    scale_x_date("Date", date_breaks = "1 month",date_labels = "%b-%d",expand=c(0,0))  +
    scale_y_continuous("US Total New Cases\n(daily)",
                       expand=c(0.02,0), breaks = seq(0,50000,10000), 
                       labels = paste0(seq(0,50,10),"k")) +
    theme(text = element_text(size=6,color="white"),
          line = element_line(size=0.25,color="white"),
          axis.ticks = element_line(size=0.25,color="white"),
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          axis.ticks.length = unit(1,"pt"),
          panel.background = element_rect(fill = "transparent",colour = NA),
          plot.background = element_rect(fill = "transparent",colour = NA),
          axis.text.x  = element_text(angle=0, hjust=0.5,color="white"), 
          axis.text.y  = element_text(angle=0, hjust=0.5,color="white"), 
          plot.margin = margin(1, 2, 1, 2, "pt")) ->
    dailyplot 
  
  ggsave(dailyplot,filename=sprintf("linetotal%d.png",counter),width=1.5,height=1, bg = "transparent")
  counter = counter + 1
  tempday = tempday + lubridate::duration("1 day")
}

#Generate inset colormap ggplot of new state cases
tempday = lubridate::ymd("2020-02-01")
finalday = lubridate::today()
counter = 1
while(tempday != finalday) {
  covid_case_data %>%
    filter(date == as.character(tempday)) %>%
    mutate(fips = as.character(fips)) %>%
    filter(!is.na(state)) %>%
    filter(!is.na(cases)) ->
  yesterday
  newdata = dplyr::left_join(shp_utm,yesterday, by = c("NAME"="state")) 
  newdata %>%
    ggplot() + 
    geom_sf(aes(fill=newcases), size=0.1, color="grey50") +
    scale_fill_gradient("New\nCases",low="#ffe2bf",high="darkred",limits=c(0,12500), breaks=seq(0,12500,2500),
                        labels=sprintf("%0.1fk",seq(0,12.5,by=2.5)),na.value="grey90") +
    scale_x_continuous(expand=c(0,0)) +
    scale_y_continuous(expand=c(0,0)) +
    theme(text = element_text(size=4,color="white"),
          line = element_line(size=0.25,color="white"),
          axis.ticks = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.ticks.length = unit(1,"pt"),
          panel.background = element_rect(fill = "transparent",colour = NA),
          plot.background = element_rect(fill = "transparent",colour = NA),
          legend.key.height =  unit(5,"pt"),
          legend.key.width =  unit(5,"pt"),
          legend.text = element_text(size=4,color="white"),
          legend.margin = margin(0, 0, 0, 0, "pt"),
          legend.box.margin = margin(0, 5, 0, -10,  "pt"),
          legend.background = element_rect(fill = "transparent",colour = NA),
          axis.text.x  = element_blank(),
          axis.text.y  = element_blank(),
          plot.margin = margin(1, 2, 1, 2, "pt")) ->
    dailyplot 
  
  ggsave(dailyplot,filename=sprintf("mapdailyl%d.png",counter),width=1.5,height=1, bg = "transparent")
  counter = counter + 1
  tempday = tempday + lubridate::duration("1 day")
}

#Append two ggplot images together + add image overlay with information
for(i in 1:148) {
  lll = image_read(sprintf("mapdailyl%d.png",i))
  image_append(c(image_read(sprintf("linetotal%d.png",i)),lll)) %>%
    image_resize("x250") %>%
    image_composite(image_read(sprintf("covid%d.png",i)),., gravity = "southwest")  %>%
    image_composite(image_read("overlayinfo.png"))  %>%
    image_write(sprintf("covidfull%d.png",i))
}

#Add title and caption
tempday = lubridate::ymd("2020-02-01")
finalday = lubridate::today()
counter = 1
for(i in 1:148) {
  rayimage::add_title(sprintf("covidfull%d.png",i), title_size = 30,
                      title_bar_color="black", title_bar_alpha = 0.5, title_color="white",
                      title_text = sprintf("Feb-June Daily New COVID-19 Cases by State: %s", 
                                           as.character(tempday))) %>% 
    rayimage::add_title(title_size = 24, title_offset = c(20,85),
                        title_color="white",
                        title_text = "Daily recorded new case count mapped to the height of the extruded state polygon in 3D.", 
                        filename=glue::glue("covidfinal{i}.png"))
  tempday = tempday + lubridate::duration("1 day")
}