ggspatial_with_inset.R

# Inset method based on https://gis.stackexchange.com/questions/222799/create-an-inset-map-in-r/222877#222877

# libraries
library(sf)
library(raster)
library(tidyverse)
library(rosm)
library(ggspatial)
library(ggsn)
library(rmapshaper)

options(stringsAsFactors = FALSE)

### MAIN MAP
# sites dataset isn't online, but its just an sfc_POINT layer with an ID field and a T/F field for whether the site
# has lab data or not
sites <- st_read(file.path(getwd(), 'sites', 'all_sites.gpkg'))

# order lab data signifier factor by false then true and sort data so lab points always plot last (on top)
all_sites <- all_sites[order(all_sites$LAB_YN), ]
all_sites$LAB_YN <- factor(all_sites$LAB_YN, levels = c(FALSE,TRUE))
# short names make the mess below slightly more readable
ss <- as(all_sites, 'Spatial')

# plot bounds. A rounded bbox can make ticks look better, and it makes positioning elements simpler
pext <- round(raster::extent(bbox(ss)),1)

sites_plot <- ggplot() +
  geom_osm(type = 'cartolight', x = pext) +
  geom_spatial(data = as(all_sites, 'Spatial'), aes(colour = LAB_YN, size = LAB_YN), alpha = 0.7) + 
  # match labels on scales or legend duplicates itself
  scale_colour_manual(values = c('#ff9600', '#214a7e'), labels = c('Description only', 'Lab data')) +
  scale_size_manual(values = c(0.3, 1.6), labels = c('Description only', 'Lab data')) +
  # default y ticks were a bit derpy
  scale_y_continuous(breaks = seq(-26.5, -24.5, 0.25)) +
  # give ggsn fns plot extent manually when using geom_spatial
  ggsn::north(x.min = pext[1], x.max = pext[2], y.min = pext[3], y.max = pext[4],
              symbol = 10, scale = 0.1,
              anchor = c(x = pext[1] + 0.5, 
                         y = pext[3] + 0.15)) +
  ggsn::scalebar(x.min = pext[1], x.max = pext[2], y.min = pext[3], y.max = pext[4],
                 anchor = c(x = pext[1] + 0.3, 
                            y = pext[3] + 0.05), 
                 height = 0.015, dist = 20, st.size = 3, dd2km = TRUE, model = 'WGS84') +
  labs(x = 'Longitude', y = 'Latitude') +
  ggtitle('Soil site locations', subtitle = 'Burnett-Mary region, South East Queensland') +
  theme(legend.position = c(0.25, 0.2),
        # use fill = alpha('white', 0) for transparent bg
        legend.background = element_rect(fill = '#f5f5f3'),
        legend.key = element_rect(size = 0),
        # use this if you don't want a legend title and don't have a transparent legend.background
        legend.title = element_blank()) +   
  coord_fixed()
  
### INSET MAP
# I used a 1:100,000 scale outline of QLD available at http://qldspatial.information.qld.gov.au/catalogue/custom/search.page?q=%22geographic+features+-+queensland+series%22
# because it was already saved to my HD, rnaturalearth would be less work
qld_outline <- st_read(file.path(dirname(getwd()), 'area_boundaries', 'state_outline.gpkg')) %>%
  mutate(STATE = 'QLD') %>%
  group_by(STATE) %>%
  summarise() %>%
  st_buffer(., dist = 0L) %>%
  as(., 'Spatial') %>%
  rmapshaper::ms_simplify(., keep = 0.0005) # trial and error for the keep factor
  
# get extent of main map and make it into a polygon
main_ext <- as(pext, 'SpatialPolygons')
  
# save inset as a grob. Matched the colour scheme to OSM Cartolight using good ol' MS Paint - screencap, paste, dropper tool
inset_map <- ggplotGrob(
  ggplot() +
  theme_void() +
  theme(plot.background = element_rect(fill = '#d1d5d7', colour = '#f5f5f3')) +
  geom_spatial(data = qld_outline, fill = '#f5f5f3', colour = NA) +
  geom_spatial(data = main_ext,    fill = '#E4535E', colour = NA, alpha = 0.7) +
  coord_fixed()
)

### FINAL COMBINATION
sites_plot_ins <- sites_plot +
  # location on plot by trial and error; aspect ratio of inset plot is preserved
  # so don't worry about squaring up the coords below
  annotation_custom(grob = inset_map, 
                    xmin = pext[1] + 0.8, ymin = pext[4] + 0.15, 
                    xmax = pext[2], ymax = pext[3] + 1.4)

#  \o/ https://github.com/obrl-soil/bits-n-pieces/blob/master/images/BM_sites_plot_titled.png \o/