library(here)
library(raster)
library(sf)
library(rnaturalearth)
library(plotly)
library(grid)
library(cowplot)
library(wesanderson)
library(tidyverse)

Read in summary stats files

spread_path <- here("output", "spreadsheets")

sumstats_100 <- read_csv(file.path(spread_path, "sumstats_100.csv"))
sumstats_200 <- read_csv(file.path(spread_path, "sumstats_200.csv"))
sumstats_300 <- read_csv(file.path(spread_path, "sumstats_300.csv"))

Read in rasters at the three resolutions to provide a template for the genetic diversity maps.

template_100 <- stack(here("data", "climate", "rasters_100km", "chelsa_100km.tif"))
template_200 <- stack(here("data", "climate", "rasters_200km", "chelsa_200km.tif"))
template_300 <- stack(here("data", "climate", "rasters_300km", "chelsa_300km.tif"))

world_basemap <- ne_countries(scale = "medium", returnclass = "sf") %>% 
  st_transform(crs = crs(template_100))

Resolution

100 km

Hill number 1 (corrected) map.

hill_1_100 <- raster(template_100)

hill_1_100[sumstats_100$cells] <- sumstats_100$hill_1_c

hill_1_100_poly <- rasterToPolygons(hill_1_100) %>% 
  st_as_sf() 

pal <- wes_palette("Zissou1", 100, type = "continuous")
hill_1_100_map <- ggplot() +
  geom_sf(data = world_basemap) +
  geom_sf(data = hill_1_100_poly, aes(fill = layer), size = 0.1) +
  scale_fill_gradientn(colors = pal,  guide = NULL) + 
  ggthemes::theme_map()


ggplotly(hill_1_100_map) 
hill_1_100_leg <- ggplot() +
  geom_sf(data = world_basemap) +
  geom_sf(data = hill_1_100_poly, aes(fill = layer), size = 0.1) +
  scale_fill_gradientn(colors = pal,  name = "Hill 1") 

hill_1_100_leg <- get_legend(hill_1_100_leg) 
ggdraw(hill_1_100_leg)

# # 
# # grid.newpage()
# # grid.draw(hill_1_100_leg)

Average pi map.

avg_pi_100 <- raster(template_100)

avg_pi_100[sumstats_100$cells] <- sumstats_100$avg_pi

avg_pi_100_poly <- rasterToPolygons(avg_pi_100) %>% 
  st_as_sf() %>% 
  filter(layer < 0.02)

pal <- wes_palette("Zissou1", 100, type = "continuous")
avg_pi_100_map <- ggplot() +
  geom_sf(data = world_basemap) +
  geom_sf(data = avg_pi_100_poly, aes(fill = layer), size = 0.1) +
  scale_fill_gradientn(colors = pal,  name = "Avg Pi") + 
  ggthemes::theme_map() +
  NULL

avg_pi_100_map

200 km

Hill number 1 (corrected) map.

hill_1_200 <- raster(template_200)

hill_1_200[sumstats_200$cells] <- sumstats_200$hill_1_c

hill_1_200_poly <- rasterToPolygons(hill_1_200) %>% 
  st_as_sf()
pal <- wes_palette("Zissou1", 200, type = "continuous")
ggplot() +
  geom_sf(data = world_basemap) +
  geom_sf(data = hill_1_200_poly, aes(fill = layer), size = 0.2) + 
  scale_fill_gradientn(colors = pal, name = "Hill 1") + 
  ggthemes::theme_map()

300

Hill number 1 (corrected) map.

hill_1_300 <- raster(template_300)

hill_1_300[sumstats_300$cells] <- sumstats_300$hill_1_c

hill_1_300_poly <- rasterToPolygons(hill_1_300) %>% 
  st_as_sf()

pal <- wes_palette("Zissou1", 300, type = "continuous")
ggplot() +
  geom_sf(data = world_basemap) +
  geom_sf(data = hill_1_300_poly,
          aes(fill = layer),
          size = 0.3) +
  scale_fill_gradientn(colors = pal, name = "Hill 1 Corrected") +
  ggthemes::theme_map()

Hill number 1 (corrected) map.

avg_pi_300 <- raster(template_300)

avg_pi_300[sumstats_300$cells] <- sumstats_300$avg_pi

avg_pi_300_poly <- rasterToPolygons(avg_pi_300) %>% 
  st_as_sf() %>% 
  filter(layer < 0.02)

pal <- wes_palette("Zissou1", 300, type = "continuous")
avg_pi_300_map <- ggplot() +
  geom_sf(data = world_basemap) +
  geom_sf(data = avg_pi_300_poly, aes(fill = layer), size = 0.1) +
  scale_fill_gradientn(colors = pal,  name = "Avg Pi") + 
  ggthemes::theme_map() +
  NULL

avg_pi_300_map

# ggplotly(avg_pi_300_map) 
# # 
# # grid.newpage()
# # grid.draw(avg_pi_300_leg)