example of spatial data which may be desirable to # overlay on a map (for Twitter thread)

example_geo_data.R

# example_geo_data.R: example of spatial data which may be desirable to
# overlay on a map.
# Cameron Patrick <cameron.patrick@unimelb.edu.au>, 22 Jan 2022.

library(tidyverse)
library(rnaturalearth)
library(ragg)
# You will also need the natural earth data:
# remotes::install_github("ropensci/rnaturalearthhires")

n_rows <- 6
n_cols <- 8
n_times <- 5

### example geographic region
aus_polys <- fortify(ne_countries(country = "Australia", scale = 10))
grid_centres <-
  expand_grid(grid_row = seq_len(n_rows),
              grid_col = seq_len(n_cols)) %>%
  mutate(grid = paste0(LETTERS[grid_col], grid_row),
         lat = -37.5 - 0.25 * grid_row + 0.125,
         long = 144 + 0.25 * grid_col - 0.125)
ggplot(aus_polys, aes(x = long, y = lat, group = group)) +
  geom_polygon(colour = "black", fill = "cornsilk") +
  geom_text(data = grid_centres,
            aes(x = long, y = lat, label = grid),
            size = 9 / .pt,
            inherit.aes = FALSE) +
  geom_hline(data = tibble(lat = seq(-37.5, -39, by = -0.25)),
             aes(yintercept = lat),
             colour = "grey60",
             linetype = "dashed",
             inherit.aes = FALSE) +
  geom_vline(data = tibble(long = seq(144, 146, by = 0.25)),
             aes(xintercept = long),
             colour = "grey60",
             linetype = "dashed",
             inherit.aes = FALSE) +
  annotate("point", x = 144.96, y = -37.81) +
  annotate("text", x = 144.96, y = -37.79, label = "Melbourne", 
           hjust = 0.5, vjust = 0, size = 11 / .pt) +
  coord_quickmap(xlim = c(144, 146),
                 ylim = c(-39, -37.5)) +
  labs(x = "Longitude", y = "Latitude") +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "lightskyblue", colour = NA),
        panel.grid = element_blank())
ggsave("fig1_map.png", width = 8, height = 6, dpi = 300, device = agg_png)

ggplot(aus_polys, aes(x = long, y = lat, group = group)) +
  geom_polygon(colour = "black", fill = "cornsilk") +
  coord_cartesian(xlim = c(144, 146),
                  ylim = c(-39, -37.5)) +
  labs(x = "Longitude", y = "Latitude") +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "lightskyblue", colour = NA),
        panel.grid = element_blank())
ggsave("fig1b_map_stretched.png",
       width = 8, height = 6, dpi = 300, device = agg_png)

### simulate some example data
set.seed(1234)
grid_data <-
  expand_grid(grid_row = seq_len(n_rows),
              grid_col = LETTERS[seq_len(n_cols)]) %>%
  mutate(grid = paste0(grid_col, grid_row),
         grid_intercept = rnorm(n(), 3, 0.5),
         grid_slope_1 = rnorm(n(), 0, 0.2),
         grid_slope_2 = rnorm(n(), -0.5, 0.2)) %>%
  filter(!(grid %in% c("A1", "H3",
                       "D3", "C4", "A5", "B5", "C5", "D5", "E5", "F5",
                       "A6", "B6", "C6", "D6", "E6", "F6", "G6")))
grid_time_data <- grid_data %>%
  expand_grid(time = letters[seq_len(n_times)]) %>%
  mutate(time_num = as.numeric(factor(time)) - 1,
         group_1 = rnorm(n(), grid_intercept + grid_slope_1 * time_num, 0.2),
         group_2 = rnorm(n(), grid_intercept + grid_slope_2 * time_num, 0.2))
full_data <- grid_time_data %>%
  pivot_longer(group_1:group_2,
               names_to = "intervention",
               values_to = "log_estimate") %>%
  mutate(estimate = exp(log_estimate),
         conf_low = exp(log_estimate - 0.4),
         conf_high = exp(log_estimate + 0.4))

### plot simulated data
full_data %>%
  ggplot(aes(x = time, y = estimate, ymin = conf_low, ymax = conf_high,
             colour = intervention, group = intervention)) +
  geom_line() +
  geom_pointrange(fatten = 1.25) +
  facet_grid(rows = vars(grid_row),
             cols = vars(grid_col)) +
  scale_colour_manual(values = c("dodgerblue3", "firebrick3")) +
  labs(x = "Time",
       y = "Mean outcome (95% CI)",
       colour = "Intervention") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        legend.position = c(0.01, 0.01),
        legend.justification = c(0, 0))
ggsave("fig2_data.png", width = 8, height = 6, dpi = 300, device = agg_png)

### current best approach
empty_grids <- tribble(
  ~grid, ~fill,
  "A1", "cornsilk",
  "H3", "cornsilk",
  "D3", "lightskyblue",
  "C4", "lightskyblue",
  "A5", "lightskyblue",
  "B5", "lightskyblue",
  "C5", "lightskyblue",
  "D5", "lightskyblue",
  "E5", "lightskyblue",
  "F5", "lightskyblue",
  "A6", "lightskyblue",
  "B6", "lightskyblue",
  "C6", "lightskyblue",
  "D6", "lightskyblue",
  "E6", "lightskyblue",
  "F6", "lightskyblue",
  "G6", "lightskyblue",
) %>%
  mutate(grid_col = str_sub(grid, 1, 1),
         grid_row = str_sub(grid, 2, 2))
full_data %>%
  ggplot(aes(x = time, y = estimate, ymin = conf_low, ymax = conf_high,
             colour = intervention, group = intervention)) +
  geom_line() +
  geom_pointrange(fatten = 1.25) +
  geom_rect(data = empty_grids,
            xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf,
            aes(fill = fill),
            inherit.aes = FALSE) +
  facet_grid(rows = vars(grid_row),
             cols = vars(grid_col)) +
  scale_fill_identity() +
  scale_colour_manual(values = c("dodgerblue3", "firebrick3")) +
  labs(x = "Time",
       y = "Mean outcome (95% CI)",
       colour = "Intervention") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        legend.position = c(0.01, 0.01),
        legend.justification = c(0, 0))
ggsave("fig2b_data_panel_bg.png", width = 8, height = 6, dpi = 300, device = agg_png)
ggsave("fig2c_data_trans_bg.png",
       plot = last_plot() + 
         theme(plot.background = element_rect(fill = "transparent", colour = NA)),
       width = 8, height = 6, dpi = 300, device = agg_png)

Here you go!

library(tidyverse)
library(rnaturalearth)
library(ragg)
library(png)
library(grid)
# You will also need the natural earth data:
# remotes::install_github("ropensci/rnaturalearthhires")

n_rows <- 6
n_cols <- 8
n_times <- 5

### example geographic region
aus_polys <- fortify(ne_countries(country = "Australia", scale = 10))
grid_centres <-
  expand_grid(grid_row = seq_len(n_rows),
              grid_col = seq_len(n_cols)) %>%
  mutate(grid = paste0(LETTERS[grid_col], grid_row),
         lat = -37.5 - 0.25 * grid_row + 0.125,
         long = 144 + 0.25 * grid_col - 0.125)


ggplot(aus_polys, aes(x = long, y = lat, group = group)) +
  geom_polygon(colour = "black", fill = "cornsilk") +
  coord_cartesian(xlim = c(144, 146),
                  ylim = c(-39, -37.5)) +
  labs(x = "", y = "") +
  scale_x_continuous(breaks = c(144.0, 146.0), 
                     labels = c("", "")) + 
  scale_y_continuous(breaks = c(-37.5, -39.0), 
                     labels = c("", "")) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "lightskyblue", colour = NA),
        panel.grid = element_blank())
ggsave("fig1b_map_stretched.png",
       width = 8, height = 6, dpi = 300, device = agg_png)



### simulate some example data
set.seed(1234)
grid_data <-
  expand_grid(grid_row = seq_len(n_rows),
              grid_col = LETTERS[seq_len(n_cols)]) %>%
  mutate(grid = paste0(grid_col, grid_row),
         grid_intercept = rnorm(n(), 3, 0.5),
         grid_slope_1 = rnorm(n(), 0, 0.2),
         grid_slope_2 = rnorm(n(), -0.5, 0.2)) %>%
  filter(!(grid %in% c("A1", "H3",
                       "D3", "C4", "A5", "B5", "C5", "D5", "E5", "F5",
                       "A6", "B6", "C6", "D6", "E6", "F6", "G6")))
grid_time_data <- grid_data %>%
  expand_grid(time = letters[seq_len(n_times)]) %>%
  mutate(time_num = as.numeric(factor(time)) - 1,
         group_1 = rnorm(n(), grid_intercept + grid_slope_1 * time_num, 0.2),
         group_2 = rnorm(n(), grid_intercept + grid_slope_2 * time_num, 0.2))
full_data <- grid_time_data %>%
  pivot_longer(group_1:group_2,
               names_to = "intervention",
               values_to = "log_estimate") %>%
  mutate(estimate = exp(log_estimate),
         conf_low = exp(log_estimate - 0.4),
         conf_high = exp(log_estimate + 0.4))



# load background image from above
img <- png::readPNG("fig1b_map_stretched.png")

# create grid of interventions
p <- full_data %>%
  ggplot(aes(x = time, y = estimate, ymin = conf_low, ymax = conf_high,
             colour = intervention, group = intervention)) +
  geom_line() +
  geom_pointrange(fatten = 1.25) +
  facet_grid(rows = vars(grid_row),
             cols = vars(grid_col)) +
  scale_colour_manual(values = c("dodgerblue3", "firebrick3")) +
  labs(x = "Time",
       y = "Mean outcome (95% CI)",
       colour = "Intervention") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        plot.background = element_rect(fill=NA),
        legend.position = c(0.01, 0.01),
        legend.justification = c(0, 0)) 


# drop empty facets
grob <- ggplotGrob(p);
idx <- which(grob$layout$name %in% c("panel-1-1", 
                                     "panel-3-4", "panel-3-8",
                                     "panel-4-3",
                                     "panel-5-1", "panel-5-2", "panel-5-3", "panel-5-4", "panel-5-5", "panel-5-6",
                                     "panel-6-1", "panel-6-2", "panel-6-3", "panel-6-4", "panel-6-5", "panel-6-6", "panel-6-7"));
for (i in idx) grob$grobs[[i]] <- nullGrob()

grid.draw(gList(rasterGrob(img, width = unit(1,"npc"), height = unit(1,"npc")), 
                grob))

# update: Remembered how to save
png("fig2b_panel_overlayed_map.png",
    width = 8, height = 6, units = "in", res = 300)
grid.draw(gList(rasterGrob(img, width = unit(1,"npc"), height = unit(1,"npc")), 
                grob))
dev.off()

Edit: Remembered how to save, so added it to the code.