Versión en español aquí
Nowadays, it is not only possible but also straightforward to make cool and informative 3d maps in R. Thanks to
rayrender by Tyler Morgan-Wall and
elevatr by Jeff Hollister, we can combine a georeferenced image and a Digital Elevation Model (DEM) to really bring an existing map to life.
Twitter user @flotsam has many great examples, such as this recent one:
A 1972 NRCan (formerly the Department of Energy, Mines, and Resources) map of Montréal, Canada. Spent more time trying to decipher cryptic filenames to find what I want than actually rayshading it. :P#rayshader adventures, an #rstats tale pic.twitter.com/O756ddhs9G— flotsam (@researchremora) August 2, 2021
The overall process for rayshading a georeferenced image is summarized step-by-step in this gist by Tyler Morgan-Wall.
There are many georeferenced maps online, and we can use GIS software to match up any digital map with its corresponding spatial context by creating ground control points. However, I like making my own thematic maps in
ggplot2 and wanted to plot one in 3d with
To make maps in
ggplot2 using popular spatial packages (
stars, etc.), we are already working with spatial data, so there should be no need to assign real-world coordinates to each pixel of the raster in the final map manually. I also kept losing my patience trying to create control points for a polynomial warp using a trackpad or the weird orange nub on my computer with my shaky hand.
This StackOverflow thread explains how to pull the information that describes the maximum and minimum xy values for the panel in a ggplot object and then assign them to the raster image we exported this object as.
If we want to avoid manual georeferencing and use this approach, the key point here (which is kind of a hack) is to make sure that the ggplot object has no margins so that the xy values that describe the panel extent match up with the overall extent of the image we export and then use for
rayshader. For this, we can use a custom theme (
theme_nothing) that comes with the
cowplot package by Claus Wilke, and the rationale is explained here. Basically, this theme is equivalent to setting various plot elements to blank with no expansion of the scales.
Let me demonstrate. This example uses spatial data from the CONABIO Geoinformation Portal and the CONANP spatial Information download area. I downloaded shapefiles for Protected Areas (Áreas Naturales Protegidas ) and Landcover types (Uso del suelo y vegetación, escala 1:250000, serie VI (continuo nacional)) and then focused on the Iztaccíhuatl-Popocatépetl protected area.
These steps are for reading the shp data, subsetting the features for the area of interest, lumping together some of the land cover categories (somewhat arbitrarily) with
case_when and making a simple map of the volcanoes with pretty colors that I chose for each land cover type. Note how
theme_nothing() removes all the margins and axis ticks/labels, and then I added extra arguments to
coord_sf() to make sure everything stays within the plot panel.
It took some trial and error to export the tiff (using a ragg device) and have no blank space around the edges, but once the output looks good we can read this new image with the
Let’s import the tiff and have a look.
The geospatial data we need is in the ggplot object, and we can pull it out with
ggplot_build and some indexing. This gives us a list with ranges for x and y, which go into the
extent function to finally georeference the map (we also need to define the projection). Next, we export the georeferenced raster.
INT1U means the values are bound between 0 and 255, and
PHOTOMETRIC=RGB is the color interpretation tag.
Let’s re-import the geoTiff
After this, the process follows most guides on rayshading images. First, we transform the raster to an array and use it to get elevation data using
elevatr. We will then crop, resize, and transpose the DEM data (because rasters and arrays are oriented differently in R).
A quick glimpse of the elevation data
Now we can compute the shadows and drape the array on the elevation data.
Here’s a snapshot of the rgl window produced by
plot3d, after fiddling with the camera values. If you are reading this in Mexico City, the area should look familiar.
An overhead view produced by
plot_map() looks pretty good too.
Finally, we can add elements with rayrender and render the fancier maps with
render_highquality(). I like this approach to add more atmosphere by using
light = FALSE and instead lighting the scene with a sphere of light placed at a custom location.
Turn off the lights and add your own🤓: #rstats #rayshader #rayrender— Tyler Morgan-Wall (@tylermorganwall) October 27, 2019
render_highquality(light = FALSE, scene_elements = sphere(y = 100, radius = 10, material = diffuse(lightintensity = 250, implicit_sample = TRUE))) pic.twitter.com/g1QDEphZGh
Here are some renders with different camera locations and light placements.
The results came out quite nice, and with this setup we could even make animations or add additional elements. This approach can also work with base R plots, as long as the plot panel goes all the way to the edges of the output raster image. I hope someone finds this helpful or entertaining. If you have any questions or feedback please let me know.