Make a map with an inset using Boston Harbor and the eastern seaboard as an example - also show sites.

make_map_with_inset.R

#' --------------------------
#' I recently had to make a map of a small area
#' with an inset showing context. So, here
#' is the code (this is for the Stone Living Lab
#' Living Seawalls), but you should be able to modify 
#' it for whatever you want, or even use rnaturalearth
#' and replace the blow-up of the small area with a vector
#' map instead of a raster. Note, you will need
#' an API key from maptiler for the vectors to work.
#' --------------------------

library(sf)
library(maptiles)
library(ggplot2)
library(terra)
library(tidyterra)
library(patchwork)
library(rnaturalearth)
library(tigris)
options(tigris_use_cache = TRUE)

####
## All the bounding boxes you could ever want
####

# east coast bounding box
east_coast <- st_bbox(c(xmin = -76.390579,
                        xmax = -63.766312,
                        ymax = 48.060084,
                        ymin = 38.538602),
                      crs = 4326)


# Boston Harbor Bounding Box

harbor <-  st_bbox(c(xmin = -71.085156,
                       xmax = -70.969441,
                       ymax = 42.403399,
                       ymin = 42.322624),
                     crs = 4326)

# sites
sites <- tibble::tribble(
  ~lat, ~lon, ~name,
  42.35490593013214, -71.0444413181874, "Fan Pier",
  42.38328982587925, -71.02949227110159, "Condor Street\n(2 LSWs)"
) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

####
## Get background maps
####

# Boston Harbor
source("scripts/stadia_key.R")

external_tiles <- create_provider(
  name = "maptiler_satellite-v2",
  url = paste0("https://api.maptiler.com/tiles/satellite-v2/{z}/{x}/{y}.jpg?KEY=",maptiler_key),
  citation = "© MapTiler © OpenStreetMap contributors"
)


boston <- get_tiles(harbor, provider = external_tiles,
                    crop = TRUE, cachedir = "tiles",
                    apikey = maptiler_key, zoom = 13,
                    project = TRUE) 

plotRGB(boston, stretch = "lin")

boston_vector <- ne_download( scale = 10L, type = "states", category = "cultural") |>
  st_make_valid() |>
  st_crop(harbor)

# Land area boundary (excludes water) for Suffolk County/Boston
boston_vector <- county_subdivisions("MA", c("Suffolk", "Middlesex", "Norfolk", "Plymouth")) |>
  erase_water() |>
  st_make_valid() |>
  st_transform(crs = 4326) |>
  st_crop(harbor) |>
  st_collection_extract("POLYGON") 


# East Coast
seaboard <- ne_states(country = "united states of america") |>
  st_crop(east_coast)


seaboard_raster <- get_tiles(east_coast, provider = external_tiles,
                    crop = TRUE, cachedir = "tiles/seaboard",
                    apikey = maptiler_key, zoom = 5,
                    project = TRUE, verbose = TRUE) 
  

####
## Make the maps
####

# eastern seaboard with box for boston
seaboard_plot <- 
  ggplot() +
  geom_sf(data = seaboard) +
  geom_sf(data = harbor |> 
            st_as_sfc(crs = 4326) |> 
            st_centroid(),
          color = "#e68613",
          size = 2)  +
  ggthemes::theme_map() +
  theme(panel.background = element_rect(fill = NULL))

seaboard_raster_plot <- ggplot() +
  geom_spatraster_rgb(data = seaboard_raster, stretch = "lin") +
  geom_sf(data = harbor |> 
            st_as_sfc(crs = 4326) |> 
            st_centroid(),
          color = "red",
          size = 2)  +
  coord_sf(expand = FALSE) +
  ggthemes::theme_map() +
  theme(panel.background = element_rect(fill = NULL),
        panel.border = element_rect(colour = "black", fill=NA, linewidth=1))

# boston harbor with boxes for condor and fan pier
harbor_plot <- 
  ggplot() +
  geom_spatraster_rgb(data = boston, stretch = "lin") +
  geom_sf(data = sites, color = "#e68613", size = 5) +
  geom_sf_label(data = sites, aes(label = name),
                nudge_y = -0.007, fill = "white",
                size = 3.5) +
  coord_sf(
    expand = FALSE
  ) +
  theme_void()
  
harbor_vector_plot <- 
  ggplot() +
  geom_sf(data = boston_vector, fill = "grey") +
  geom_sf(data = sites, color = "#e68613", size = 5) +
  geom_sf_label(data = sites, aes(label = name),
                nudge_y = -0.007, fill = "white",
                size = 3.5) +
  coord_sf(expand = FALSE) +
  theme_void() +
  theme(panel.background = element_rect(fill = "lightblue"))

####
## Combine the maps
####
harbor_plot +
  inset_element(free(seaboard_plot, type = "panel"),
                left   = 1 - 1/2.8,  
                right  = 1+0.1,
                bottom = 0-0.02,
                top    = 1/2.8,
                align_to = "panel", clip = TRUE)
# save out
ggsave("figures/site_map.jpg")


# all rasters
harbor_plot +
  inset_element(free(seaboard_raster_plot, type = "panel"),
                left   = 1 - 1/2.7,  
                right  = 1+0.031,
                bottom = 0-0.02,
                top    = 1/2.7,
                align_to = "panel", clip = TRUE)
# save out
ggsave("figures/site_map_rasters.jpg")

# all vectors
harbor_vector_plot +
  inset_element(free(seaboard_plot, type = "panel"),
                left   = 1 - 1/2.8,  
                right  = 1+0.11,
                bottom = 0-0.02,
                top    = 1/2.8,
                align_to = "panel", clip = TRUE)

ggsave("figures/site_map_vectors.jpg")