Here is the script:
# Packages
pkgs <- c("sf", "sp", "terra", "raster", "elevatr", "rayshader")
to_install <- setdiff(pkgs, rownames(installed.packages()))
if (length(to_install)) install.packages(to_install, quiet = TRUE)
invisible(lapply(pkgs, library, character.only = TRUE))
# --- Mount Fuji AOI (25 km buffer) ---
# Summit coords (lon, lat)
fuji_ll <- sf::st_sfc(sf::st_point(c(138.7274, 35.3606)), crs = 4326)
crs_fuji <- "EPSG:32654"
fuji_utm <- sf::st_transform(fuji_ll, crs_fuji) # UTM Zone 54N (meters)
aoi_utm <- sf::st_buffer(fuji_utm, dist = 25000) # 25 km radius
aoi_ll <- sf::st_transform(aoi_utm, 4326) |> sf::st_as_sf()
# --- Fetch DEM with elevatr, already projected to UTM (meters) ---
dem_raster <- elevatr::get_elev_raster(
locations = aoi_ll, # pass polygon in UTM
z = 12, # zoom; raise to 13 for more detail (larger download)
clip = "bbox")
# Convert to terra and clip/mask to AOI for a clean edge
dem <- terra::rast(dem_raster)
dem <- terra::project(dem, crs_fuji)
terra::plot(dem)
# Height matrix for rayshader
elmat <- rayshader::raster_to_matrix(raster::raster(dem))
# Set zscale to horizontal resolution in meters (good physical default)
zscale <- mean(terra::res(dem)) # meters per pixel
zscale
ray <- rayshader::ray_shade(
elmat,
sunaltitude = 40, # elevation angle of the sun
sunangle = 315, # azimuth (NW light)
zscale = zscale,
multicore = TRUE
)
amb <- rayshader::ambient_shade(
elmat,
zscale = zscale,
multicore = TRUE
)
lamb <- rayshader::lamb_shade(
elmat,
sunaltitude = 45,
sunangle = 315, # WNW
zscale = zscale
)
png("fuji_desert_base.png", width = 3425, height = 3400, res = 300)
par(mar = rep(0, 4))
elmat |>
rayshader::sphere_shade(texture = "desert") |>
# rayshader::add_shadow(ray, 0.6) |>
rayshader::plot_map()
dev.off()
png("fuji_desert_ray.png", width = 3425, height = 3400, res = 300)
par(mar = rep(0, 4))
elmat |>
rayshader::sphere_shade(texture = "desert") |>
rayshader::add_shadow(ray, 0.6) |>
rayshader::plot_map()
dev.off()
png("fuji_desert_lamb.png", width = 3425, height = 3400, res = 300)
par(mar = rep(0, 4))
elmat |>
rayshader::sphere_shade(texture = "desert") |>
rayshader::add_shadow(lamb, 0.6) |>
rayshader::plot_map()
dev.off()
png("fuji_desert_amb.png", width = 3425, height = 3400, res = 300)
par(mar = rep(0, 4))
elmat |>
rayshader::sphere_shade(texture = "desert") |>
rayshader::add_shadow(amb, 0.6) |>
rayshader::plot_map()
dev.off()