To begin, we load the necessary packages. These libraries will help us with raster data manipulation, geospatial transformations, 3D rendering, and brightness/contrast adjustments.
# install.packages(c("rayshader", "raster", "sp"))
library(rayshader)
## Warning: package 'rayshader' was built under R version 4.4.2
library(sp)
## Warning: package 'sp' was built under R version 4.4.2
library(raster)
## Warning: package 'raster' was built under R version 4.4.2
library(scales)
library(sf)
## Warning: package 'sf' was built under R version 4.4.2
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(magick)
## Warning: package 'magick' was built under R version 4.4.2
## Linking to ImageMagick 6.9.12.98
## Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fontconfig, x11
library(showtext)
## Warning: package 'showtext' was built under R version 4.4.2
## Loading required package: sysfonts
## Warning: package 'sysfonts' was built under R version 4.4.2
## Loading required package: showtextdb
## Warning: package 'showtextdb' was built under R version 4.4.2
library(grid)
library(ragg)
## Warning: package 'ragg' was built under R version 4.4.2
font_add_google("Merriweather", "merriweather")
font_add_google("Playfair Display", "playfair")
font_add_google("Lora", "lora")
font_add_google("Arvo", "arvo")
showtext_auto()
Here, we load elevation data from two SRTM tiles. These tiles are merged into a single dataset to cover the region of interest. We then visualize the elevation using height shading to check the data.
# Load SRTM elevation tiles
elevation1 = raster::raster("S09E117.SRTMGL1.hgt/S09E117.hgt")
elevation2 = raster::raster("S09E118.SRTMGL1.hgt/S09E118.hgt")
# Merge tiles into one seamless elevation dataset
elevation = raster::merge(elevation1, elevation2)
# Visualize the elevation data
height_shade(raster_to_matrix(elevation)) %>%
plot_map()
Next, we load satellite imagery as separate red, green, and blue bands. These are stacked together to create a single RGB image. Finally, we visualize the raw satellite imagery using plotRGB to ensure it loaded correctly.
# Load RGB bands of satellite imagery
r = raster::raster("LC09_L1TP_115066_20240713_20240713_02_T1/LC09_L1TP_115066_20240713_20240713_02_T1_B2.TIF")
g = raster::raster("LC09_L1TP_115066_20240713_20240713_02_T1/LC09_L1TP_115066_20240713_20240713_02_T1_B3.TIF")
b = raster::raster("LC09_L1TP_115066_20240713_20240713_02_T1/LC09_L1TP_115066_20240713_20240713_02_T1_B4.TIF")
# Stack bands to form an RGB image
imagery = raster::stack(r, g, b)
# Visualize
raster::plotRGB(imagery, scale=255^2)
To improve the appearance of the satellite imagery, we apply a gamma correction. This adjusts the brightness and contrast to better match human perception. We use the square root transformation to achieve this and visualize the corrected result.
# Apply gamma correction to brighten the imagery
imagery_corrected <- sqrt(raster::stack(r, g, b))
# Visualize
raster::plotRGB(imagery_corrected)
Our elevation data uses a geographic CRS (latitude/longitude), while the imagery uses a projected CRS (UTM). To align them, we transform the elevation data to match the CRS of the imagery.
# Check CRS of the imagery
raster::crs(r)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 50N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",117,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16050]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]
# Check CRS of the elevation data
raster::crs(elevation)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]
# Transform elevation data to match imagery's CRS
elevation_utm <- raster::projectRaster(elevation, crs = crs(r), method = "bilinear")
Using the bounding box of the area of interest, we crop both the elevation and imagery datasets. This ensures that all data is aligned spatially and focuses only on the target region.
# Define bounding box for the region of interest (ROI)
top_right <- c(x = 118.1023, y = -8.13384)
bottom_left <- c(x = 117.8759, y = -8.360502)
# Convert bounding box to UTM and crop the datasets
extent_latlong <- sp::SpatialPoints(rbind(bottom_left, top_right), proj4string = sp::CRS("+proj=longlat"))
extent_utm <- sp::spTransform(extent_latlong, raster::crs(elevation_utm))
e <- raster::extent(extent_utm)
# Crop datasets to the ROI
imagery_cropped <- raster::crop(imagery_corrected, e)
elevation_cropped <- raster::crop(elevation_utm, e)
Now, we convert the cropped elevation and imagery datasets into formats suitable for 3D plotting. The elevation data is converted to a matrix, while the RGB imagery is processed into a normalized array.
# Convert elevation to a matrix
el_matrix <- rayshader::raster_to_matrix(elevation_cropped)
# Prepare RGB imagery as a normalized array
r_cropped <- rayshader::raster_to_matrix(imagery_cropped[[1]])
g_cropped <- rayshader::raster_to_matrix(imagery_cropped[[2]])
b_cropped <- rayshader::raster_to_matrix(imagery_cropped[[3]])
imagery_array <- array(0, dim = c(nrow(r_cropped), ncol(r_cropped), 3))
imagery_array[,,1] <- r_cropped / 255
imagery_array[,,2] <- g_cropped / 255
imagery_array[,,3] <- b_cropped / 255
imagery_array = aperm(imagery_array, c(2,1,3))
imagery_contrast = scales::rescale(imagery_array, to=c(0,1))
plot_map(imagery_contrast)
# Prepare the imagery as an RGB array
if (inherits(imagery_contrast, "RasterStack") || inherits(imagery_contrast, "RasterBrick")) {
r <- rayshader::raster_to_matrix(imagery_contrast[[1]])
g <- rayshader::raster_to_matrix(imagery_contrast[[2]])
b <- rayshader::raster_to_matrix(imagery_contrast[[3]])
# Create RGB array for rayshader
imagery_array <- array(0, dim = c(nrow(r), ncol(r), 3))
imagery_array[,,1] <- r / max(r)
imagery_array[,,2] <- g / max(g)
imagery_array[,,3] <- b / max(b)
} else {
imagery_array <- imagery_contrast # Use directly if already in RGB array format
}
We use a custom function to adjust the brightness and contrast of the imagery array. This ensures the 3D visualization has clear and vibrant colors.
# Adjust brightness and contrast for the imagery array
# Define a brightness and contrast adjustment function
adjust_brightness_contrast <- function(img, brightness = 1, contrast = 1) {
img <- img * contrast
img <- img + (brightness - 1)
img <- pmin(pmax(img, 0), 1) # Ensure values stay within [0, 1]
return(img)
}
# Convert the elevation matrix to a matrix if it is not already
if (inherits(el_matrix, "RasterLayer")) {
elevation_matrix <- raster_to_matrix(el_matrix)
} else {
elevation_matrix <- el_matrix # Use directly if already a matrix
}
# Prepare the imagery as an RGB array
if (inherits(imagery_contrast, "RasterStack") || inherits(imagery_contrast, "RasterBrick")) {
r <- rayshader::raster_to_matrix(imagery_contrast[[1]])
g <- rayshader::raster_to_matrix(imagery_contrast[[2]])
b <- rayshader::raster_to_matrix(imagery_contrast[[3]])
# Create RGB array for rayshader
imagery_array <- array(0, dim = c(nrow(r), ncol(r), 3))
imagery_array[,,1] <- r / max(r)
imagery_array[,,2] <- g / max(g)
imagery_array[,,3] <- b / max(b)
} else {
imagery_array <- imagery_contrast # Use directly if already in RGB array format
}
# Adjust brightness and contrast for the imagery array
imagery_array <- adjust_brightness_contrast(imagery_array, brightness = 1, contrast = 2.5)
plot_map(imagery_array)
To enhance the 3D model’s realism, we add shadow and ambient shading layers based on the elevation matrix. These are combined with the imagery to create the final map.
# Generate shadow and ambient layers
shadow_layer <- ray_shade(el_matrix, zscale = 15, lambert = TRUE)
ambient_layer <- ambient_shade(el_matrix, zscale = 15)
# Combine imagery with shading
combined_map <- height_shade(el_matrix) %>%
add_overlay(imagery_array, alphalayer = 0.8) %>%
add_shadow(shadow_layer, 0.25) %>%
add_shadow(ambient_layer, 0.15)
Finally, we create a 3D plot of Mount Tambora using the combined map and elevation data. The render_snapshot function includes a description of the site.
# Plot 3D map
plot_3d(
combined_map, el_matrix,
zscale = 12, fov = 0, theta = -45, phi = 45,
windowsize = c(3000, 2400), zoom = 0.75,
baseshape = 'circle'
)
render_highquality(filename = "output/3d_tambora_base.png", clear = TRUE)
# Load the rendered 3D image
image_3d <- image_read("output/3d_tambora_base.png")
# Export to PNG using ragg
agg_png("output/3D_Tambora_Annotated.png", width = 3000, height = 1835, res = 300)
# Start a new plot
grid.newpage()
grid.rect(gp = gpar(fill = NA, col = NA))
grid.raster(image_3d) #overlay 3d
# Overlay Title
grid.text(
"Tambora Volcano",
x = unit(0.03, "npc"), y = unit(0.95, "npc"),
just = "left", gp = gpar(fontfamily = "merriweather", fontsize = 70, col = "#333333", fontface = "bold")
)
# Overlay Subtitle
grid.text(
"Located on Indonesia's Sumbawa Island, Mount Tambora erupted in 1815, marking the largest volcanic event\nin recorded history. The explosion caused thousands of deaths and disrupted global weather, leading to the\n'Year Without a Summer' in 1816, characterized by widespread crop failures and food shortages worldwide.",
x = unit(0.03, "npc"), y = unit(0.865, "npc"),
just = "left", gp = gpar(fontfamily = "sans", fontsize = 42.5, col = "#555555", lineheight=0.3)
)
# Overlay Data Source
grid.text(
"Data Source: USGS, URS Earthdata NASA \n© Godfried Junio Matahelemual",
x = unit(0.975, "npc"), y = unit(0.05, "npc"),
just = "right", gp = gpar(fontfamily = "sans", fontsize = 25, col = "#777", lineheight=0.3)
)
# Close the device to save the file
dev.off()
## png
## 2
# Load the image
annotated_image <- png::readPNG("output/3D_Tambora_Annotated.png")
# Display the final result
grid.newpage()
grid.raster(annotated_image)