3D maps
Authors: Iuliia Burdun, Anto Aasa
November 30, 2022
In this example, we will learn about 3-D visualisation in R using the
package rayshader. The code used in this tutorial was
modified from https://www.rayshader.com/. First, we will create a
digital elevation model (DEM) with a 1-m spatial resolution. Further, we
will add an aerial image at a 20-cm spatial resolution above the DEM.
Finally, we will make an animation with 3-D visualisation. The study
area is a Rummu quarry (59.227398, 24.193968) - a man-made heap that
originated due to the mining of Vasalemma limestone and marble. This
place is very popular during the hot season since it is possible to swim
and dive in transparent water (link).
Lat’s start with packages installation. We have to install packages
from GitHub (not from CRAN! ), otherwise, you will face the problem with
plotting. Package remotes is designed to enable downloading
the packages from GitHub repositories.
install.packages("remotes")
library(remotes)
remotes::install_github("tylermorganwall/rayshader")
remotes::install_github("dmurdoch/rgl")
install.packages("magick") # package implements functionality for manipulating
# high-dimensional arrays using efficient vectorised methodsCall the packages we will use:
library(magick)
library(rayshader)
library(raster)
library(rayrender)Download the DEM data from Landboard Estonia. link
Select 1-m resolution GeoTIFF data for map sheet number 63613. After the file is downloaded - place it to your working directory.
Set your working directory. This should be the folder where you can save the outputs.
setwd("/YOUR/WORKING/DIRECTORY")Import downloaded raster:
DEM <- raster::raster( "63613_dem_1m.tif") # It works if file "63613_dem_1m.tif" is located in your working directory
crs(DEM) # see the coordinate system## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=lcc +lat_0=57.5175539305556 +lon_0=24 +lat_1=59.3333333333333
## +lat_2=58 +x_0=500000 +y_0=6375000 +ellps=GRS80 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["Estonian Coordinate System of 1997 (with axis order normalized for visualization)",
## BASEGEOGCRS["EST97",
## DATUM["Estonia 1997",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4180]],
## CONVERSION["Estonian National Grid",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",57.5175539305556,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",24,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",59.3333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",58,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",6375000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (Y)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (X)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]]]Plot the DEM file
plot(DEM,col = terrain.colors(n=50))
DEM covers a big area while we need a small part of it. Rummu quarry is the small area in the bottom left corner of the figure. Create the polygon to clip DEM to the area of interest:
polygon <- as(extent(510800, 511300, 6565300, 6565750), 'SpatialPolygons')
crs(polygon) <- crs(DEM) # assign projection to the polygon
crs(polygon) # polygon and DEM have the same coordinate system## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=lcc +lat_0=57.5175539305556 +lon_0=24 +lat_1=59.3333333333333
## +lat_2=58 +x_0=500000 +y_0=6375000 +ellps=GRS80 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["Estonian Coordinate System of 1997 (with axis order normalized for visualization)",
## BASEGEOGCRS["EST97",
## DATUM["Estonia 1997",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4180]],
## CONVERSION["Estonian National Grid",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",57.5175539305556,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",24,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",59.3333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",58,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",6375000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (Y)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (X)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]]]DEM_crop <- crop(DEM, polygon) # crop DEM by polygon
plot(DEM_crop,col = terrain.colors(n=50))
To use rayshader package, we need to convert DEM from
raster to a matrix
DEM_matrix <- raster_to_matrix(DEM_crop)Let’s try the first image:
DEM_matrix %>%
sphere_shade(texture = "desert") %>% # `imhof1`,`imhof2`,`imhof3`,`imhof4`,`desert`, `bw`
plot_map()