Code to create a joyplot of Mt. Washington for the 30 Day Map Challenge

washington_joyplot.R

library(terrainr) # remotes::install.github("mikemahoney218/terrainr")
library(raster)
library(ggplot2)
library(ggridges)
library(sf)

mt_washington <- tmaptools::geocode_OSM("mt washington")
mt_washington <- get_bbox(lat = mt_washington$bbox[c("ymin", "ymax")],
                          lng = mt_washington$bbox[c("xmin", "xmax")])
mt_washington <- set_bbox_side_length(mt_washington, 16000)
mt_washington <- get_tiles(mt_washington)
mt_washington["merged"] <- merge_rasters(mt_washington[[1]])
mt_washington <- raster::raster(mt_washington[["merged"]])
nlines <- 40
line_grid <- expand.grid(
  # this bit is hardcoded but corresponds to raster::extent(mt_washington)
  # extracting these in this format is harder programatically than just 
  # copying and pasting
  x = c(-71.40366, -71.20271),
  y = seq(44.18405, 44.35685, length.out = nlines)
)

outfile <- vector("list", nlines)
for (i in 1:nlines) {
  cur_points <- line_grid[((i * 2) - 1):(i * 2), ]
  linefile <- st_linestring(as.matrix(cur_points))
  linefile <- as(st_sfc(linefile, crs = 4326), "Spatial")
  elevations <- raster::extract(mt_washington, linefile)
  if (length(elevations[[1]]) == 0) {
    next
  } else {
    outfile[[i]] <- data.frame(
    xstart = rep(min(cur_points$x), length(elevations)),
    xend = rep(max(cur_points$x), length(elevations)),
    ystart = rep(min(cur_points$y), length(elevations)),
    yend = rep(max(cur_points$y), length(elevations)),
    elev = as.vector(elevations)[[1]]
    )
  }
}

washingtonelev <- do.call(rbind, outfile)
latvals <- unique(washingtonelev$ystart)
washout <- vector("list", length(latvals))
for (i in seq_along(latvals)) {
  washout[[i]] <- data.frame(
  lon = seq(min(washingtonelev$xstart), 
            max(washingtonelev$xend), 
            length.out = nrow(washingtonelev[washingtonelev$ystart == latvals[[i]], ])),
  lat = rep(latvals[[i]], 
            nrow(washingtonelev[washingtonelev$ystart == latvals[[i]], ])),
  elev = washingtonelev[washingtonelev$ystart == latvals[[i]], ][["elev"]]
)
}
washout <- do.call(rbind, washout)
ggplot(washout, aes(x = lon, y = lat, group = lat, height = elev)) +
  scale_y_discrete(expand = c(0, 0)) +
  scale_x_continuous(expand = c(0, 0), name = "Mt Washington") +
  geom_density_ridges(stat = "identity", scale = 15, fill="black", color = "white") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "black"),
        axis.line = element_blank(),
    axis.text.x=element_blank(),
    plot.background = element_rect(fill = "black"),
    axis.ticks.x=element_blank(),
    axis.title.y=element_blank(),
    axis.text.y=element_blank(),
    axis.ticks.y=element_blank()) + 
  coord_map() + 
  ggsave("30_map.png", width = 8, height = 11, units = "in")