Skip to content

Instantly share code, notes, and snippets.

@mikemahoney218
Last active December 1, 2020 18:01
Show Gist options
  • Save mikemahoney218/ffc33fa8ef455abeaff64306263852fb to your computer and use it in GitHub Desktop.
Save mikemahoney218/ffc33fa8ef455abeaff64306263852fb to your computer and use it in GitHub Desktop.
Code to create a joyplot of Mt. Washington for the 30 Day Map Challenge
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment