Regional locator inset of the Western US

Rendered Map Plate

Map 3 Western US

Map 3: Western United States regional elevation model with whole degree graticule and populated places

Technical Summary

  • Objective: Serve as the widest regional locator panel (Inset 1 of the Map 7 tryptic), establishing continental context for Southern California.
  • Key Packages: terra, tidyverse, tidyterra, geojsonsf
  • Data Sources:
    • Western United States 90-meter Digital Elevation Model (source_data/dem90_hf/dem90_hf.tif)
    • Global shaded relief hillshade (source_data/global_raster/GRAY_HR_SR_OB.tif)
    • California Populated Places vector (source_data/tl_2023_06_place/)
    • California State Boundary (source_data/california_boundary.geojson)
    • Southern California locator box (source_data/socal_aoi.geojson)
  • Styling Decisions:
    • Portrait 3 wide x 4 tall aspect ratio (width = 3, height = 4).
    • Axes labeled in whole degrees (e.g. 125°W, 120°W, 115°W and 32°N, 36°N, 40°N, 44°N).
    • Populated places vector overlay rendered with subtle transparency.
    • Red outline delineates the next zoom extent (Map 4).

Step-by-Step Cartographic Workflow & Intermediate Outputs

The script scripts/map_03.r sets up the continental zoom level for the multi-scale locator sequence.

Step 1: Downsampling Regional DEM & Ocean Shaded Relief Cropping

The high-resolution 90m Western US elevation model is aggregated by a factor of 5 to optimize rendering speeds. The global Natural Earth ocean hillshade (GRAY_HR_SR_OB.tif) is cropped to the California overview bounding box (cali_overview.geojson), and coordinate reference systems are unified.

# Load and downsample DEM for fast regional plotting
zoom_1_dem <- rast("source_data/dem90_hf/dem90_hf.tif")
zoom_1_dem <- aggregate(zoom_1_dem, fact = 5)

# Crop global shaded relief hillshade to West Coast extent
zoom_1_extent <- geojson_sf("source_data/cali_overview.geojson")
zoom_1_hillshade <- crop(world, zoom_1_extent)
zoom_1_dem <- project(zoom_1_dem, crs(zoom_1_hillshade))
zoom_1_dem <- crop(zoom_1_dem, zoom_1_hillshade)

Step 2: Intermediate Plot — Base Terrain & Next-Zoom Inset Box

An alpha-blended hillshade raster is paired with a viridis terrain gradient. A red bounding box is added to indicate the sub-regional scope of Map 4 (The Bight of California).

# Initial raster overlay test with bounding box
ggplot() +
  geom_raster(data = zoom_1_dem_df, aes(x = x, y = y, fill = dem90_hf)) +
  scale_fill_viridis_c() +
  geom_raster(data = zoom_1_hillshade_df, aes(x = x, y = y, alpha = GRAY_HR_SR_OB)) +
  scale_alpha(range = c(0.3, 0.6), guide = "none") +
  geom_spatvector(data = zoom_2_crop_extent, color = "red", lwd = 1.5, fill = NA) +
  my_theme +
  coord_sf()

Map 3.1 Intermediate Output

Intermediate Output 3.1: Raw terrain elevation and ocean hillshade

Map 3.2 Graticule Test

Intermediate Output 3.2: Adding subtle graticule guidelines

Step 3: Vector Context & Whole-Degree Geographic Coordinate Labels

To orient the viewer, California Census populated places (tl_2023_06_place.shp) are simplified (simplifyGeom) and rendered in semi-transparent light gray (#EEEEEE50). The state perimeter is traced in crisp black (linewidth = 0.7).

Finally, custom coordinate tick formatters convert raw decimal degrees into clear geographical markers: * Longitude: seq(-125, -115, by = 5) → 125°W, 120°W, 115°W * Latitude: seq(32, 44, by = 4) → 32°N, 36°N, 40°N, 44°N

# Finalizing coordinate breaks and subtle census polygons
zoom_1_final_plot <- ggplot() +
  geom_raster(data = zoom_1_dem_df, aes(x = x, y = y, fill = dem90_hf)) +
  scale_fill_viridis_c(guide = "none") +
  geom_raster(data = zoom_1_hillshade_df, aes(x = x, y = y, alpha = GRAY_HR_SR_OB)) +
  scale_alpha(range = c(0.05, 0.3), guide = "none") +
  geom_spatvector(data = cal_poly, fill = NA, color = "black", linewidth = 0.7) +
  geom_spatvector(data = places_simple, fill = NA, color = "#EEEEEE50", linewidth = 0.2) +
  geom_spatvector(data = zoom_2_crop_extent, color = "red", lwd = 1, fill = NA) +
  scale_x_continuous(breaks = seq(-125, -115, by = 5), labels = function(x) paste0(abs(x), "°W")) +
  scale_y_continuous(breaks = seq(32, 44, by = 4), labels = function(y) paste0(y, "°N")) +
  coord_sf(expand = FALSE) +
  labs(
    title = "Pacific Coast of California",
    subtitle = "Western United States"
  )

Map 3 Final Inset

Intermediate Output: Final portrait rendering for Map 7 top tryptic

Full R Script (scripts/map_03.r)

Show R Code 
# map 3-4-5 is a tryptic,
# that runs horizontally on top of map 7
# zooming in to campus
# ie: the zoom to Cali locator sheet

# each of the 3 maps is a grayscale alpha hillshade
# over (some color-scheme) DEMs,
# -- just as in the overlay episode --
#  -- which overlay episode Jonnie?     --


# clean the environment and hidden objects
rm(list = ls())

library(terra)
library(geojsonsf)
# library(sf)
library(ggplot2)
library(tidyterra)
library(dplyr)
# library(ggpubr)
# library(raster)

# make sure output window is 1x1
# because you muck with it a lot
par(mfrow = c(1, 1))

# We'll need a grayscale palette later
grays <- colorRampPalette(c("black", "white"))(255)

# set map number
current_sheet <- 3
# set ggplot counter
current_ggplot <- 0

# our auto ggtitle maker
gg_labelmaker <- function(plot_num) {
  gg_title <- c("Map:", current_sheet, " ggplot:", plot_num)
  plot_text <- paste(gg_title, collapse = " ")
  print(plot_text)
  current_ggplot <<- plot_num
  return(plot_text)
}
# every ggtitle should be:
# ggtitle(gg_labelmaker(current_ggplot+1))
# end automagic ggtitle           #######

# we will use this to make the graticule and axes
my_theme <- theme(
  axis.title.x = element_blank(),
  axis.title.y = element_blank(),
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  legend.position = "none",
  panel.ontop = TRUE,
  panel.grid.major = element_line(color = "#FFFFFF33"),
  panel.background = element_blank()
)

# set up a local CRS to use throughout
campus_DEM <- rast("source_data/campus_DEM.tif")
campus_crs <- crs(campus_DEM)
str(campus_crs)


# #####################
# Map 3
# Zoom 1: west US overview

# hillshade:
world <- rast("source_data/global_raster/GRAY_HR_SR_OB.tif")
plot(world)

# dem:
zoom_1_dem <- rast("source_data/dem90_hf/dem90_hf.tif")
plot(zoom_1_dem)
# this is a bit too big later on. let's downsample it:
zoom_1_dem <- aggregate(zoom_1_dem, fact = 5)
plot(zoom_1_dem)

# and get our clipping extent
zoom_1_extent <- geojson_sf("source_data/cali_overview.geojson")
zoom_1_extent_vect <- vect(zoom_1_extent)
zoom_1_extent_vect
zoom_1_extent


# clip the world hillshade to zoom_1
zoom_1_hillshade <- crop(x = world, y = zoom_1_extent)
plot(zoom_1_hillshade)

crs(zoom_1_hillshade)
zoom_1_extent_vect <- project(x = zoom_1_extent_vect, y = zoom_1_hillshade)

# clip the dem
crs(zoom_1_dem) == crs(zoom_1_hillshade)
zoom_1_dem <- project(zoom_1_dem, crs(zoom_1_hillshade))
crs(zoom_1_dem) == crs(zoom_1_hillshade)

zoom_1_dem <- crop(x = zoom_1_dem, y = zoom_1_hillshade)

plot(zoom_1_dem)
plot(zoom_1_hillshade)


# we'll need this locator box socal_aoi.geojson
# is the next crop extent.
# it came from a Planet harvest that we did
zoom_2_crop_extent <- geojson_sf("source_data/socal_aoi.geojson")
zoom_2_crop_extent <- vect(zoom_2_crop_extent)


# now let's set up for ggplot
# and alphas
#################################################

# first check that projections match
crs(zoom_1_hillshade) == crs(zoom_1_dem)
crs(zoom_1_dem) == crs(zoom_2_crop_extent)
crs(zoom_1_hillshade) == crs(zoom_2_crop_extent)
# the rasters don't match the vector.
# but we still seem to be able to overlay.


# hillshade as ggplot
str(zoom_1_hillshade)
zoom_1_hillshade_df <- as.data.frame(zoom_1_hillshade, xy = TRUE)
colnames(zoom_1_hillshade_df)

zoom_1_hillshade_plot <- ggplot() +
  geom_raster(
    data = zoom_1_hillshade_df,
    aes(x = x, y = y, fill = GRAY_HR_SR_OB)
  ) +
  scale_fill_continuous() +
  coord_sf() +
  ggtitle("Zoom 1 Hillshade", subtitle = gg_labelmaker(current_ggplot + 1))

# zoom_1_hillshade_plot

# zoom 1 DEM as ggplot
str(zoom_1_dem)
zoom_1_dem_df <- as.data.frame(zoom_1_dem, xy = TRUE)
colnames(zoom_1_dem_df)

zoom_1_dem_plot <- ggplot() +
  geom_raster(
    data = zoom_1_dem_df,
    aes(x = x, y = y, fill = dem90_hf)
  ) +
  scale_fill_viridis_c() +
  coord_sf() +
  ggtitle("Zoom1 DEM", subtitle = gg_labelmaker(current_ggplot + 1))

# zoom_1_dem_plot


# hillshade as alpha
zoom_1_hillshade_plot <- ggplot() +
  geom_raster(
    data = zoom_1_hillshade_df,
    aes(x = x, y = y, alpha = GRAY_HR_SR_OB)
  ) +
  scale_alpha(range = c(0.3, 0.6), guide = "none") +
  coord_sf() +
  ggtitle("Zoom 1 Hillshade Alpha", subtitle = gg_labelmaker(current_ggplot + 1))

# zoom_1_hillshade_plot


# now overlay them
zoom_1_overlay_plot <- ggplot() +
  geom_raster(
    data = zoom_1_dem_df,
    aes(x = x, y = y, fill = dem90_hf)
  ) +
  scale_fill_viridis_c() +
  geom_raster(
    data = zoom_1_hillshade_df,
    aes(x = x, y = y, alpha = GRAY_HR_SR_OB)
  ) +
  scale_alpha(range = c(0.3, 0.6), guide = "none") +
  geom_spatvector(data = zoom_2_crop_extent, color = "red", lwd = 1.5, fill = NA) +
  coord_sf() +
  my_theme +
  ggtitle("Western US Fancy Overlay", subtitle = gg_labelmaker(current_ggplot + 1))
# zoom_1_overlay_plot


# the hillshade for the water is very subtle.
# zoom_1_overlay_plot

# add the graticule
zoom_1_overlay_plot <- ggplot() +
  geom_raster(
    data = zoom_1_dem_df,
    aes(x = x, y = y, fill = dem90_hf)
  ) +
  scale_fill_viridis_c() +
  geom_raster(
    data = zoom_1_hillshade_df,
    aes(x = x, y = y, alpha = GRAY_HR_SR_OB)
  ) +
  scale_alpha(range = c(0.3, 0.6), guide = "none") +
  geom_spatvector(data = zoom_2_crop_extent, color = "red", lwd = 1.5, fill = NA) +
  my_theme +
  coord_sf() +
  ggtitle("Western US Fancy Overlay", subtitle = gg_labelmaker(current_ggplot + 1))

# zoom_1_overlay_plot
# ggsave("images/map3.1.png", width = 3, height = 4, plot = last_plot(), bg = "white")

# now seize control of labels
# zoom_1_overlay_plot
# ggsave("images/map3.2.png", width = 3, height = 4, plot = last_plot(), bg = "white")


# now let's add
# 'california populated places'
# which is census data
# for the sake of a nice vizualization
places <- vect("source_data/tl_2023_06_place/tl_2023_06_place.shp")
places <- project(places, crs(zoom_1_dem))
places <- crop(places, zoom_1_hillshade)
places_simple <- simplifyGeom(places, tolerance = 0.005)

# California state boundary polygon
cal_poly <- vect("source_data/california_boundary.geojson")
cal_poly <- project(cal_poly, crs(zoom_1_dem))

# now one to save
zoom_1_final_plot <- ggplot() +
  geom_raster(
    data = zoom_1_dem_df,
    aes(x = x, y = y, fill = dem90_hf)
  ) +
  scale_fill_viridis_c(guide = "none") +
  geom_raster(
    data = zoom_1_hillshade_df,
    aes(x = x, y = y, alpha = GRAY_HR_SR_OB)
  ) +
  scale_alpha(range = c(0.05, 0.3), guide = "none") +
  geom_spatvector(data = cal_poly, fill = NA, color = "black", linewidth = 0.7) +
  geom_spatvector(data = places_simple, fill = "NA", color = "#EEEEEE50", linewidth = 0.2) +
  geom_spatvector(data = zoom_2_crop_extent, color = "red", lwd = 1, fill = NA) +
  scale_x_continuous(breaks = seq(-125, -115, by = 5), labels = function(x) paste0(abs(x), "°W")) +
  scale_y_continuous(breaks = seq(32, 44, by = 4), labels = function(y) paste0(y, "°N")) +
  theme(
    plot.title = element_text(size = 9.5, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 8, hjust = 0.5),
    axis.title = element_blank(),
    axis.text = element_text(size = 8, color = "gray20"),
    axis.ticks = element_line(color = "gray40", linewidth = 0.3),
    legend.position = "none",
    panel.ontop = TRUE,
    panel.grid.major = element_line(color = "#FFFFFF44", linetype = "dashed"),
    panel.background = element_blank()
  ) +
  coord_sf(expand = FALSE) +
  labs(
    title = "Pacific Coast of California",
    subtitle = "Western United States"
  )

zoom_1_final_plot

#######################################################################################


ggsave("images/map3.png", width = 3, height = 4, plot = zoom_1_final_plot, bg = "white")
ggsave("final_output/map_03.png", width = 3, height = 4, plot = zoom_1_final_plot, bg = "white")