In annecori/mRIIDSprocessData: Epidemiological forcasting using real-time data

output: pdf_document: citation_package: natbib keep_tex: true fig_caption: true latex_engine: pdflatex title: Mapping the risk profile author: - name: Sangeeta Bhatia affiliation: Imperial College London abstract: keywords: date: "r format(Sys.time(), '%B %d, %Y')" geometry: margin=1in fontfamily: mathpazo fontsize: 11pt spacing: double bibliography: biblio-style: apsr endnote: no ---

knitr::opts_chunk$set(fig.width = 12,
                      fig.height = 6,
                      echo = FALSE,
                      warning = FALSE,
                      message = FALSE,
                      fig.path = "figures/")

library(dplyr)
library(ggplot2)
library(ggthemes)
library(sf)
library(ggrepel)
library(patchwork)

Read in the weighted risk

Read in the weighted risk.

wtd_risk_pops <- here::here("output",
                            "weighted_rel_risk_pops_alt.csv") %>%
    readr::read_csv() 

wtd_risk_pops$adm0_iso3 <- countrycode::countrycode(wtd_risk_pops$adm0,
                                               "country.name",
                                               "iso3c")


wtd_risk_area <- here::here("output",
                            "weighted_rel_risk_area_alt.csv") %>%
    readr::read_csv() 


wtd_risk_area$adm0_iso3 <- countrycode::countrycode(wtd_risk_area$adm0,
                                               "country.name",
                                               "iso3c")

Add the relative risk across the subdivisions?

cum_risk_pops <- group_by(wtd_risk_pops, adm0) %>%
    summarise(cum_risk = sum(wtd_rel_risk, na.rm = TRUE)) 
cum_risk_pops$adm0_iso3 <- countrycode::countrycode(cum_risk_pops$adm0,
                                                    "country.name",
                                                    "iso3c")



cum_risk_area <- group_by(wtd_risk_area, adm0) %>%
    summarise(cum_risk = sum(wtd_rel_risk, na.rm = TRUE)) 
cum_risk_area$adm0_iso3 <- countrycode::countrycode(cum_risk_area$adm0,
                                                    "country.name",
                                                    "iso3c")

DRC

DRC has to read in separately because of the mismatch between GADM and meta-data files.

drc <- here::here("data/Geography/drc_moritz/shp",
                  "congo_angola.shp") %>%
    st_read() 

drc <- filter(drc, ISO == "COD")

drc_2 <- split(drc, drc$NAME_2)
dims <- purrr::map(drc_2, nrow)
idx <- which(dims > 0)
drc_2 <- drc_2[idx] %>%
    lapply(st_union) %>%
    lapply(st_sf) 

adm2 <- names(drc_2)
drc_popsrisk <- purrr::map(adm2, function(x) {
    shp <- drc_2[[x]]
    shp <- rename(shp, "geometry" = "X..i..")
    shp$wtd_rel_risk <- filter(wtd_risk_pops, orig_name == x) %>%
        pull(wtd_rel_risk)
    shp
})

drc_arearisk <- purrr::map(adm2, function(x) {
    shp <- drc_2[[x]]
    shp <- rename(shp, "geometry" = "X..i..")
    shp$wtd_rel_risk <- filter(wtd_risk_area, orig_name == x) %>%
        pull(wtd_rel_risk)
    shp
})

Prepare the names of the other files to be read at ADM2 and ADM1 resolutions respectively.

Matching by population

pops_adm1_neighbors <- c("Angola", "Burundi",
                         "Uganda", "Republic of Congo",
                         "Central African Republic",
                         "South Sudan",
                         "Rwanda",
                         "Zambia") %>%
    countrycode::countrycode("country.name", "iso3c")

pops_shp_dirs <- paste0("data/Geography/gadm36_", pops_adm1_neighbors, "_shp")
names(pops_shp_dirs) <- pops_adm1_neighbors
pops_shp_adm1 <- paste0("gadm36_", pops_adm1_neighbors, "_1.shp")

Read in files.

pops_drc_nghbrs <- purrr::map2(pops_shp_dirs,
                               pops_shp_adm1, 
                               ~ st_read(here::here(.x, .y)))

Now join the data sets.

pops_drc <- filter(wtd_risk_pops, adm0_iso3 == "COD") %>%
    left_join(drc, ., by = c("NAME_2" = "orig_name"))


pops_nghbrs_adm1_risk <- purrr::map(pops_adm1_neighbors, function(x) {
    shp <- pops_drc_nghbrs[[x]]
    risk <- filter(wtd_risk_pops, adm0_iso3 == x)
    shp <- left_join(shp, risk, by = c("NAME_1" = "orig_name"))
    shp
})

Matching by area

First fix DRC

area_drc <- filter(wtd_risk_area, adm0_iso3 == "COD") %>%
    left_join(drc, ., by = c("NAME_2" = "orig_name"))

area_adm1_neighbors <- c("Angola",
                         "Central African Republic",
                         "Republic of Congo",
                         "South Sudan",
                         "Uganda",                         
                         "Zambia") %>%
    countrycode::countrycode("country.name", "iso3c")
##                    "Tanzania")

area_adm0_neighbors <- c("Burundi", 
                         "Rwanda") %>%
    countrycode::countrycode("country.name", "iso3c")
drc_nghbrs_iso <-  c(area_adm1_neighbors, area_adm0_neighbors)

area_shp_dirs <- paste0("data/Geography/gadm36_", drc_nghbrs_iso, "_shp")
names(area_shp_dirs) <- drc_nghbrs_iso
area_shp_adm1 <- paste0("gadm36_", area_adm1_neighbors, "_1.shp")
area_shp_adm0 <- paste0("gadm36_", area_adm0_neighbors, "_0.shp")
area_shp_files <- c(area_shp_adm1, area_shp_adm0)

Read in files.

area_drc_nghbrs <- purrr::map2(area_shp_dirs,
                               area_shp_files, 
                               ~ st_read(here::here(.x, .y)))

And join the data sets.

area_nghbrs_adm1_risk <- purrr::map(area_adm1_neighbors, function(x) {
    shp <- area_drc_nghbrs[[x]]
    risk <- filter(wtd_risk_area, adm0_iso3 == x)
    shp <- left_join(shp, risk, by = c("NAME_1" = "orig_name"))
    shp
})


area_nghbrs_adm0_risk <- purrr::map(area_adm0_neighbors, function(x) {
    shp <- area_drc_nghbrs[[x]]
    risk <- filter(wtd_risk_area, adm0_iso3 == x)
    shp <- left_join(shp, risk, by = c("NAME_0" = "orig_name"))
    shp
})

area_all_nghbrs_risk <- c(area_nghbrs_adm1_risk, area_nghbrs_adm0_risk)

Read in boundaries

shp_adm0 <- paste0("gadm36_", drc_nghbrs_iso, "_0.shp")
drc_nghbrs_boundary <- purrr::map2(area_shp_dirs,
                                   shp_adm0, 
                                   ~ st_read(here::here(.x, .y)))

Cumulative risk

nghbrs_cumrisk_area <- purrr::map(drc_nghbrs_iso, function(x) {
    shp <- drc_nghbrs_boundary[[x]]
    risk <- filter(cum_risk_area, adm0_iso3 == x)
    shp <- left_join(shp, risk, by = c("NAME_0" = "adm0"))
    shp
})

nghbrs_cumrisk_pops <- purrr::map(drc_nghbrs_iso, function(x) {
    shp <- drc_nghbrs_boundary[[x]]
    risk <- filter(cum_risk_pops, adm0_iso3 == x)
    shp <- left_join(shp, risk, by = c("NAME_0" = "adm0"))
    shp
})

Read in centroids to place labels

centroids_lon <- purrr::map_dfr(drc_nghbrs_boundary,
                                 ~st_centroid(.x)$geometry[[1]][1])

centroids_lon <- tidyr::gather(centroids_lon, adm0_iso3, centroid_lon)
centroids_lat <- purrr::map_dfr(drc_nghbrs_boundary,
                                 ~st_centroid(.x)$geometry[[1]][2])
centroids_lat <- tidyr::gather(centroids_lat, adm0_iso3, centroid_lat)

centroids_adm0 <- left_join(centroids_lat, centroids_lon)
drc_labels <- data.frame(adm0_iso3 = "COD",
                         centroid_lat = -4,
                         centroid_lon = 21)

centroids_adm0 <- rbind(centroids_adm0, drc_labels)
centroids_adm0$country_name <- countrycode::countrycode(centroids_adm0$adm0_iso3,
                                                        origin = "iso3c",
                                                        destination = "country.name")

idx <- which(centroids_adm0$country_name == "Congo - Brazzaville")
centroids_adm0$country_name[idx] <- "Republic of Congo"

idx <- which(centroids_adm0$country_name == "Congo - Kinshasa")
centroids_adm0$country_name[idx] <- "Democratic Republic of the Congo"

## manually adjust some labels
idx <- which(centroids_adm0$adm0_iso3 == "RWA")
centroids_adm0$centroid_lon[idx] <- 32
centroids_adm0$centroid_lat[idx] <- -2.5

idx <- which(centroids_adm0$adm0_iso3 == "BDI")
centroids_adm0$centroid_lon[idx] <- 31.5

idx <- which(centroids_adm0$adm0_iso3 == "ZMB")
centroids_adm0$centroid_lat[idx] <- -15

idx <- which(centroids_adm0$adm0_iso3 == "COD")
centroids_adm0$centroid_lon[idx] <- 22

A plain map with only country outline and label.

drc_outline <- st_union(drc) %>% st_sf
p <- ggplot() + geom_sf(data = drc_outline, fill = "white")
p <- p + theme(
             panel.ontop = FALSE,
             panel.grid = element_blank(), 
             line = element_blank(), 
             rect = element_blank())
for (df in drc_nghbrs_boundary) {
    p <- p + geom_sf(data = df, fill = "white")
}
p <- p + coord_sf(datum = NA)
p <- p + geom_text_repel(data = centroids_adm0,
                         arrow = arrow(length = unit(0.03, "npc"),
                                        type = "closed",
                                        ends = "first"),
                          aes(label = country_name,
                              x = centroid_lon,
                              y = centroid_lat),
                          size = 3,
                          fontface = "bold",
                          inherit.aes = FALSE)
p <- p + xlab("") + ylab("")
ggpubr::ggexport(p, filename = "central_africa.png", res = 100)

Plot

Common stuff

m1 <- filter(cum_risk_pops,
             adm0 != "Democratic Republic of the Congo") %>%
    pull(cum_risk)

m2 <- filter(cum_risk_area,
             adm0 != "Democratic Republic of the Congo") %>%
    pull(cum_risk)

max_risk <- max(m1, m2)
p <- ggplot()
p <- p + geom_sf(data = drc_nghbrs_boundary[[1]], 
                 col = "blue", lwd = 1.05)
p <- p + geom_sf(data = drc_nghbrs_boundary[[2]],
                 col = "blue", lwd = 1.05)
p <- p + geom_sf(data = drc_nghbrs_boundary[[3]],
                 col = "blue", lwd = 1.05)
p <- p + geom_sf(data = drc_nghbrs_boundary[[4]],
                 col = "blue", lwd = 1.05)
p <- p + geom_sf(data = drc_nghbrs_boundary[[5]],
                 col = "blue", lwd = 1.05)
p <- p + geom_sf(data = drc_nghbrs_boundary[[6]],
                 col = "blue", lwd = 1.05)
p <- p + geom_sf(data = drc_nghbrs_boundary[[7]],
                 col = "blue", lwd = 1.05)
p <- p + geom_sf(data = drc_nghbrs_boundary[[8]],
                 col = "blue", lwd = 1.05)
p <- p + scale_fill_distiller(palette = "Spectral",
                              name = "Relative Risk",
                              limits = c(0, max_risk))
##                              labels = waiver())
p <- p + xlab("") + ylab("")

Matching by population.

p1 <- p
p1 <- p1 + theme(
               panel.ontop = FALSE,
               panel.grid = element_blank(), 
               line = element_blank(), 
               rect = element_blank())
for (df in drc_popsrisk) {
    p1 <- p1 + geom_sf(data = df,
                       aes(fill = wtd_rel_risk),
                       lwd = 0,
                       alpha = 0.8) 
}
for (df in pops_nghbrs_adm1_risk) {
    p1 <- p1 + geom_sf(data = df, aes(fill = wtd_rel_risk),
                     lwd = 0.1,
                     alpha = 0.8)
}
p1 <- p1 + coord_sf(datum = NA)
p1 <- p1 + geom_label_repel(data = centroids_adm0,
                          arrow = arrow(length = unit(0.03, "npc"),
                                        type = "closed",
                                        ends = "first"),
                          aes(label = adm0_iso3,
                              x = centroid_lon,
                              y = centroid_lat),
                          size = 3,
                          fontface = "bold",
                          inherit.aes = FALSE)

Same as before, just plot aggregated risk.

p2 <- p
p2 <- p2 + theme(
               panel.ontop = FALSE,
               panel.grid = element_blank(), 
               line = element_blank(), 
               rect = element_blank(), 
               text = element_blank())
for (df in drc_popsrisk) {
    p2 <- p2 + geom_sf(data = df,
                     aes(fill = wtd_rel_risk),
                     lwd = 0,
                     alpha = 0.8) 
}
for (df in nghbrs_cumrisk_pops) {
    p2 <- p2 + geom_sf(data = df, aes(fill = cum_risk),
                     lwd = 0.1,
                     alpha = 0.8)
}
p2 <- p2 + coord_sf(datum = NA)
p2 <- p2 + geom_label_repel(data = centroids_adm0,
                          arrow = arrow(length = unit(0.03, "npc"),
                                        type = "closed",
                                        ends = "first"),
                          aes(label = adm0_iso3,
                              x = centroid_lon,
                              y = centroid_lat),
                          size = 3,
                          fontface = "bold",
                          inherit.aes = FALSE)

Now for area

p3 <- p
p3 <- p3 + theme(
               panel.ontop = FALSE,
               panel.grid = element_blank(), 
               line = element_blank(), 
               rect = element_blank(), 
               text = element_blank())
for (df in drc_arearisk) {
    p3 <- p3 + geom_sf(data = df,
                     aes(fill = wtd_rel_risk),
                     lwd = 0,
                     alpha = 0.8) 
}

for (df in area_all_nghbrs_risk) {
    p3 <- p3 + geom_sf(data = df, aes(fill = wtd_rel_risk),
                     lwd = 0.1,
                     alpha = 0.8)
}
p3 <- p3 + coord_sf(datum = NA)
p3 <- p3 + geom_label_repel(data = centroids_adm0,
                          arrow = arrow(length = unit(0.03, "npc"),
                                        type = "closed",
                                        ends = "first"),
                          aes(label = adm0_iso3,
                              x = centroid_lon,
                              y = centroid_lat),
                          size = 3,
                          fontface = "bold",
                          inherit.aes = FALSE)

Same as before, just plot aggregated risk.

p4 <- p
p4 <- p4 + theme(
               panel.ontop = FALSE,
               panel.grid = element_blank(), 
               line = element_blank(), 
               rect = element_blank(), 
               text = element_blank())
for (df in drc_arearisk) {
    p4 <- p4 + geom_sf(data = df,
                     aes(fill = wtd_rel_risk),
                     lwd = 0,
                     alpha = 0.8) 
}

for (df in nghbrs_cumrisk_area) {
    p4 <- p4 + geom_sf(data = df, aes(fill = cum_risk),
                     lwd = 0.1,
                     alpha = 0.8)
}
p4 <- p4 + coord_sf(datum = NA)
p4 <- p4 + geom_label_repel(data = centroids_adm0,
                          arrow = arrow(length = unit(0.03, "npc"),
                                        type = "closed",
                                        ends = "first"),
                          aes(label = adm0_iso3,
                              x = centroid_lon,
                              y = centroid_lat),
                          size = 3,
                          fontface = "bold",
                          inherit.aes = FALSE)

Finally

ggpubr::ggarrange(p1, p2, p3, p4,
                  ncol=2, nrow=2, common.legend = TRUE, legend="bottom") %>%
    ggpubr::ggexport(filename = "relative_risks2.pdf")

ggpubr::ggarrange(p1, p2,
                  ncol = 2, nrow = 1, common.legend = TRUE, legend="bottom") %>%
    ggpubr::ggexport(filename = "bypops_relative_risks.pdf")

ggpubr::ggarrange(p3, p4,
                  ncol = 2, nrow = 1, common.legend = TRUE, legend="bottom") %>%
    ggpubr::ggexport(filename = "byarea_relative_risks.pdf")

library(cowplot)
prow <- plot_grid(p1 + theme(legend.position = "none"),
                  p2 + theme(legend.position = "none"),
                  p3 + theme(legend.position = "none"),
                  p4 + theme(legend.position = "none"),
                  align = 'vh',
                  labels = c("A", "B", "C", "D"),
                  hjust = -1,
                  nrow = 2
                  )

legend <- get_legend(p1 + theme(legend.position = "bottom"))
p5 <- plot_grid(prow, legend, ncol = 1, rel_heights = c(1, .2))
save_plot("relative_risks_cowplot.pdf", p5)