3D maps in R

The purpose of this post is to show how make a 3D map using rayshader. The idea is to use a DEM obtained from the SRTM mission (worldwide cover) or an available DEM from other National Institutions, such as Mexico’s National Institute of Statistics and Geography to obtain a 3D relief of the area of interest; and an RGB composite to overlay it on the relief. The cloudless RGB mosaic I am using in this post was obtained from Landsat-8 images and was produced using Google Earth Engine. Finally, a shapefile with the locations of two urban settlements in the area is loaded to show their location.

First load the libraries we are going to use.

library(rayshader)
library(raster)
library(sf)
library(dplyr)

Then, read the images.

dem <- raster("dem_more.tif")
locs <- st_read("Morelia.shp")
rgb <- stack("rgb.tif")

Then, clip the rgb image values to the lower and upper 3 % quantiles so the image can be better appreciated.

quants <- quantile(as.vector(values(rgb)), c(0.03, 0.97))

# Clip values to lower and upper 3 %
rgb[rgb < quants[1]] <- quants[1]
rgb[rgb > quants[2]] <- quants[2]
rgb <- ((rgb - quants[1])* 1/ quants[2]) + 0

Next, let’s transform the rgb into an array and the dem into matrix, so they can be rendered in 3D.

rgb_array = as.array(rgb)
dem_mat <-  raster_to_matrix(dem)

Extract the coordinates of the two human settlements used as reference (as numbers) and drop the geometry.

locs <- locs |>
  mutate(x = st_coordinates(locs)[,1],
         y = st_coordinates(locs)[,2]) |>
  st_drop_geometry()

Afterward, do the 3d rendering with the rgb and dem layers and set additional parameters for the visualization. Subsequently, add the labels of the two localities.

rgb_array |>
plot_3d(dem_mat, 
        zscale = 12, 
        fov = 0, 
        theta = 20, 
        zoom = 0.65, 
        # azimut
        phi = 45, 
        windowsize = c(1000, 800)) 
  render_label(dem_mat,
               long = locs$x[1],
               lat = locs$y[1],
               # altitude = 120000,
               zscale = 19,
               extent = raster::extent(dem),
               text = locs$Nombre[1],
               linecolor = "white",
               textcolor = "white")
  render_label(dem_mat,
               long = locs$x[2],
               lat = locs$y[2],
               # altitude = 120000,
               zscale = 14,
               extent = raster::extent(dem),
               text = locs$Nombre[2],
               linecolor = "white",
               textcolor = "white")

Finally, save a snapshot of the 3d model.

render_snapshot(filename = "Morelia3D.png",
                gravity = "North")

The result:

3D map of Patzcuaro and Morelia surroundings in Michoacán, Mexico.

If you are familiar with the surrounding of Morelia, Michoacán, Mexico, you will immediatly recognize Patzcuaro and Cuitzeo lakes, as well as some hills, such as the Quinceo.