Creating a stacked map in R using ggplot2 and gganimate

stacked_map.R

library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(ggnewscale) # Multiple Fill and Colour Scales in 'ggplot2'
library(gganimate) # A Grammar of Animated Graphics
library(hrbrthemes) # Additional Themes, Theme Components and Utilities for 'ggplot2'
library(rgdal) # Bindings for the 'Geospatial' Data Abstraction Library
library(pals) # Color Palettes, Colormaps, and Tools to Evaluate Them
library(transformr) # Polygon and Path Transformations


# Download Colombia Map ---------------------------------------------------

col.map <- readOGR("/home/mauricio/Documents/Website/covid_3layermap/colombia_dep/colombia_dep.shp")

# Order polygons
col.map@data$id <- rownames(col.map@data)

departamentos <-
  col.map@data %>% 
  dplyr::select(id, Cod.Dep = DPTO, Departamento = NOMBRE_DPT) %>% 
  mutate(Cod.Dep = as.numeric(Cod.Dep),
         Departamento = str_to_title(Departamento))

# Reduce resolution shapefile
col.map <- rgeos::gSimplify(col.map, tol=0.01, topologyPreserve=FALSE)

# Datos Covid Colombia
covid.col <- read_csv("https://www.datos.gov.co/api/views/gt2j-8ykr/rows.csv?accessType=DOWNLOAD")

date_update <- max(covid.col$`fecha reporte web`) %>% as.Date()
 
covid.col <-
covid.col %>%
  dplyr::select(Cod.Dep = `Codigo departamento`,
                Tipo = atención,
                Sexo) %>% 
  mutate(Sexo = toupper(Sexo),
         Tipo = ifelse(Tipo == "Recuperado" | Tipo == "Fallecido", Tipo, "Activo")) %>% 
  group_by(Cod.Dep, Tipo, Sexo) %>% 
  summarise(Casos = n()) %>% 
  pivot_wider(names_from = Sexo, values_from = Casos) %>% 
  ungroup() #%>% with(sum(F, na.rm = T) + sum(M, na.rm = T))

# Expand grid create dummy dataframe
covid.col <- 
  expand_grid(Cod.Dep = unique(covid.col$Cod.Dep), 
              Tipo = unique(covid.col$Tipo)) %>% 
  full_join(covid.col, by = c("Cod.Dep", "Tipo"))

# Replace NA for zeros
covid.col[is.na(covid.col)] <- 0

# Calculate Total cases
covid.col <- 
  covid.col %>% 
  mutate(Casos = F + M) 

# Juntar datos
covid.col <- 
left_join(departamentos, covid.col, by = "Cod.Dep")

# Summary cases
resumen.casos <- 
covid.col %>% group_by(Tipo) %>% 
  summarise(min = min(Casos),
            max = max(Casos),
            total = sum(Casos))


# Create Ranges by case ---------------------------------------------------

Activo <- 
covid.col %>% 
  filter(Tipo == "Activo") %>%
  mutate(rango = cut(Casos, breaks = pretty(Casos, 10), 
                        label = paste(na.exclude(lag(pretty(Casos, 10))), "-", na.exclude(lead(pretty(Casos, 10)))),
                     include.lowest = T
                        )) 
Fallecido <- 
  covid.col %>% 
  filter(Tipo == "Fallecido") %>%
  mutate(rango = cut(Casos, breaks = pretty(Casos, 10), 
                     label = paste(na.exclude(lag(pretty(Casos, 10))), "-", na.exclude(lead(pretty(Casos, 10)))),
                     include.lowest = T
  )) 

Recuperado <- 
  covid.col %>% 
  filter(Tipo == "Recuperado") %>%
  mutate(rango = cut(Casos, breaks = pretty(Casos, 10), 
                     label = paste(na.exclude(lag(pretty(Casos, 10))), "-", na.exclude(lead(pretty(Casos, 10)))),
                     include.lowest = T
  )) 

covid.col <- 
rbind(Activo, Recuperado, Fallecido)

mapa.covid <- 
fortify(col.map, region = "id") %>% dplyr::select(-order, -hole, -piece)

mapa.covid <- 
mapa.covid %>% left_join(covid.col, by = "id")


# Translation to center of coordinate system ------------------------------

mapa.covid <- 
  mapa.covid %>% 
  mutate(long = long - mean(long),
         lat = lat - mean(lat))

# Shear effect and Rotation Matrix ------------------------------

# Shear Matrix
shear_matrix <- 
  matrix(c(1,0,0.7,0.5),2,2) %>% t

# Rotation matrix
rotation_matrix <- 
  function(theta){
    matrix(c(cos(theta*pi/180), -sin(theta*pi/180),
             sin(theta*pi/180), cos(theta*pi/180)),
           nrow = 2, ncol = 2, byrow = T) %>% 
      round(digits = 2)
  }


final.data <- list()

for (theta in seq(5,360,5)) {
  # Application of the linear transformations
  xy <- 
    mapa.covid %>% 
    dplyr::select(long, lat) %>% 
    as.matrix() %*% rotation_matrix(theta) %*% shear_matrix 
  
  final.data[[theta]] <- mapa.covid %>% 
    mutate(long = xy[,1],
           lat = xy[,2],
           theta = theta)
}

mapa.covid <- final.data %>% bind_rows()


# Plot --------------------------------------------------------------------
mapa.covid <- 
mapa.covid %>% mutate(Tipo = factor(Tipo, levels = c("Fallecido", "Activo", "Recuperado")))

# mapa.covid <- 
# mapa.covid %>%  filter(theta == 45)

plotCOVID <-
ggplot(mapping = aes(x = long)) +
  labs(title = "COVID-19 en Colombia", 
       subtitle = paste("Casos de coronavirus por departamento. En total se han reportado",
                        scales::comma(sum(resumen.casos$total)), "en todo el país.",
                        "\nÚltima actualización:", date_update), 
       caption = "Creado por: @Mauricio_Cely") +
  geom_polygon(data = mapa.covid %>% filter(Tipo == "Fallecido"),
               aes(y = lat + 8, group = group, fill = Casos), color = "black", size = .1) +
  scale_fill_gradientn(colours = brewer.reds(10)[5:9], name = "Fallecidos") +
  new_scale_fill() +
  geom_polygon(data = mapa.covid %>% filter(Tipo == "Activo"),
               aes(y = lat, group = group, fill = Casos), color = "black", size = .1) +
  scale_fill_gradientn(colours = kovesi.linear_kry_5_95_c72(10)[c(9:7, 4)], name = "Activo") +
  new_scale_fill() +
  geom_polygon(data = mapa.covid %>% filter(Tipo == "Recuperado"),
               aes(y = lat - 8,group = group, fill = Casos), color = "black", size = .1) +
  scale_fill_gradientn(colours = ocean.algae(10)[c(9:5, 4:1)], name = "Recuperados") +
  geom_text(aes(x = -15, y = 0), size=6, color="gray65", 
            label = paste0("Activos:\n", scales::comma(resumen.casos$total[1]))) +         # layer 2
  geom_text(aes(x = -15, y = 8), size=6, color="gray65", 
            label = paste0("Fallecidos:\n", scales::comma(resumen.casos$total[2]))) + # layer 1
  geom_text(aes(x = -15, y = -8), size=6, color="gray65", 
            label = paste0("Recuperados:\n", scales::comma(resumen.casos$total[3]))) +     # layer 3
  coord_quickmap(xlim = c(-18,13)) +
  theme_ft_rc(grid = F, plot_title_size = 20) +
  theme(legend.position = "right",
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.title = element_blank()) + 
  transition_states(theta)

  # ggsave(last_plot(), 
  #        filename = "featured.jpg", 
  #        device = "jpeg",
  #        width = 16, 
  #        height = 14.55, scale = .5)
  
options(gganimate.dev_args = list(res = 115))

plotCOVID %>% animate(fps = 30,
                      duration = 10,
                      # rewind = T,
                      width = 900,
                      height = 830,  
                      renderer = gifski_renderer("COVID_COL.gif"),
                      type = "cairo")