library(latex2exp)
#### Define area
lat_range <- c(-5.0, -6.0)
lon_range <- c(15, 20)
#lon_center <- -5.4
#lat_center <- 36.75
# And the window around this centroid
lat_radius <- 0.05
lon_radius <- 0.1
gen_random_table <- function(world_num=1, return_data=FALSE, rand_seed=NULL) {
if (!is.null(rand_seed)) {
set.seed(rand_seed)
} else {
rand_seed <- sample(1:9999, size=1)
set.seed(rand_seed)
}
lat_center <- runif(1, min=min(lat_range), max=max(lat_range))
lon_center <- runif(1, min=min(lon_range), max=max(lon_range))
lon_lower <- lon_center - lon_radius
lon_upper <- lon_center + lon_radius
lat_lower <- lat_center - lat_radius
lat_upper <- lat_center + lat_radius
coords <- data.frame(x = c(lon_lower, lon_lower, lon_upper, lon_upper),
y = c(lat_lower, lat_upper, lat_lower, lat_upper))
coords_sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326)
# Using geodata
elev <- geodata::elevation_3s(lon = lon_center, lat = lat_center, res=2.5, path = getwd())
elev <- terra::crop(elev, coords_sf)
return(elev)
}
compute_hillshade <- function(elev) {
# Calculate hillshade
slopes <- terra::terrain(elev, "slope", unit = "radians")
aspect <- terra::terrain(elev, "aspect", unit = "radians")
hs <- terra::shade(slopes, aspect)
return(hs)
}
plot_hillshade <- function(elev, hs, title=NULL) {
## Plot hillshading as basemap
# (here using terra::plot, but could use tmap)
# lat_upper_trimmed <- lat_upper - 0.001*lat_upper
# lat_lower_trimmed <- lat_lower - 0.001*lat_lower
base_plot <- terra::plot(
hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE, mar=c(0,0,1,0), grid=TRUE
)
# xlim = c(elev_xmin, elev_xmax), ylim = c(elev_ymin, elev_ymax))
# overlay with elevation
color_vec <- terrain.colors(25)
plot(elev, col = color_vec, alpha = 0.5, legend = FALSE, axes = FALSE, add = TRUE)
# add contour lines
terra::contour(elev, col = "grey30", add = TRUE)
if (!is.null(title)) {
world_name <- TeX(sprintf(r"(World $\omega_{%d} \in \Omega$)", world_num))
title(world_name)
}
}
elev_r1 <- gen_random_table(rand_seed=6810)
hs_r1 <- compute_hillshade(elev_r1)
plot_hillshade(elev_r1, hs_r1)