Creates a shaded relief of NEUS

shaded_neus.R

# devtools::install_github("tidyverse/ggplot2")
library(ggplot2)
# library(dplyr)
library(RCurl)
# library(rnaturalearth)
library(sf)
# library(purrr)
library(data.table)
# library(reshape2)
# library(formula.tools)


script <- getURL("https://gist.githubusercontent.com/slarge/1642c59602308e50d97a29895cef5c3b/raw/af08f9e7c88f2a1f943eb4d3e106be2b7f68dd41/geom_relief.R", ssl.verifypeer = FALSE)
eval(parse(text = script))

xmin = -77
xmax = -65
ymin = 35
ymax = 45

xlims <- c(xmin, xmax)
ylims <- c(ymin, ymax)
crs <-  "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 

## North America layer
ne_countries <- rnaturalearth::ne_countries(scale = 10,
                                            continent = "North America",
                                            returnclass = "sf") %>% 
  sf::st_transform(crs = crs)

# 200 m isobath layer
nesbath <- marmap::getNOAA.bathy(lon1 = xmin, lon2 = xmax, 
                                 lat1 = ymin, lat2 = ymax,
                                 resolution = 1,
                                 keep = TRUE)

bathy_df  <- marmap::fortify.bathy(nesbath)
bathy_df <- as.data.table(bathy_df)

bathy_df[, c("light", "dark") := .(ifelse(z > 0, "lemonchiffon1", "azure1"),
                                    ifelse(z > 0, "gray30", "dodgerblue4"))] 
pp <- ggplot() +
  geom_relief(data = bathy_df, aes(x = x, y = y, z = z, light = light, dark = dark),
              raster = TRUE, interpolate = TRUE, sun.angle = 60) +
  # geom_sf(data = ne_countries, color = "black", fill = "transparent", size = 0.25) +
  coord_sf(xlim = c(xmin, xmax),
           ylim = c(ymin, ymax),
           expand = F) +
  theme_void()

ggsave(pp, filename = "shaded_NEUS.png", width = 8, height = 8, units = "in", dpi = 320)

lines_nm <- USAboundaries::us_states(states = "New Mexico")

nmbath <-  marmap::getNOAA.bathy(lon1 = -109.5, lon2 = -102.5, 
                                 lat1 = 31, lat2 = 37.5,
                                 resolution = 1,
                                 keep = TRUE)

bathy_nm  <- marmap::fortify.bathy(nmbath)
clip_nm <- sf::st_as_sf(bathy_nm, coords = c("x","y"),
                         crs = crs)
# bathy_nm$z <- bathy_empty$z
clip_nm <- sf::st_intersection(sf::st_union(clip_nm), sf::st_union(lines_nm))
clip_nm <- sf::st_intersection(clip_nm, lines_nm)
clip_nm <- fortify(clip_nm)
bathy_nm <- as.data.table(bathy_nm)

bathy_nm[, c("light", "dark") := .(ifelse(z > 0, "lemonchiffon1", "lemonchiffon1"),
                                   ifelse(z > 0, "gray30", "gray30"))] 

pd <- ggplot() +
  # geom_relief(data = bathy_nm, aes(x = x, y = y, z = z, light = light, dark = dark),
              raster = TRUE, interpolate = TRUE, sun.angle = 60) +
  geom_sf(data = clip_nm, color = "black", fill = "transparent", size = 0.25) +
  coord_sf(xlim = c(-102.5, -109.5),
           ylim = c(31, 37.5),
           expand = F) +
  theme_void()
pd