R Studio’s Spatial Capabilities going 3D!

Purpose

This tutorial prominently centres on the spatial abilities of R Studio; an open software program with its own scripting language that is often used for data analysis, statistics and graphing of very large datasets. For university students, R Studio is generally used for its ability to analyze, graph and model datasets to determine relevant correlations and outliers.

However, the software contains libraries capable of manipulating spatial data, making it an asset for geography students looking to model and manipulate large quantities of data!

The spatial data package offers a wide range of tools that can take data and create custom maps with a myriad of customization. By no means is this package an equivalent to ArcMap or other such dedicated mapping platforms but R can do many of the same things with their own advantages.

For more about the spatial data package follow this link: https://cran.r-project.org/web/packages/sp/sp.pdf

This tutorial demonstrates how to create 3D maps and images using the R language – if you're not familiar with R or some of its other mapping capabilities the following tutorials are provided below:

Julia Riddick -> R Studio's Spatial Capabilities at http://gracilis.carleton.ca/CUOSGwiki/index.php/R_Studio%27s_Spatial_Capabilities

Uzayr Siddiqui -> Introducing Geoprocessing Capabilities of SAGA in R Environment using RSAGA (Saga + Rstudio) at http://gracilis.carleton.ca/CUOSGwiki/index.php/Introducing_Geoprocessing_Capabilities_of_SAGA_in_R_Environment_using_RSAGA_(Saga_%2B_Rstudio)

I'm quite fond of R & RStudio and saw how it was under represented within the tutorials over the years. The software isn't as visually 'flashy' as other platforms but much more fun to handle and manipulate at a coding level.

Objective

The objective of this tutorial is to develop the users R Studio skills, competence in the R language and to demonstrate some of the various ways that 3D maps and images can be created.

The goal is to get users thinking on how R can do more than just manipulate statistics. This tutorial is aimed to get users to think beyond tables and lists and start manipulating data in 3 dimensions.

Software Introduction

Rstudio

This is an open software program that creates a user-friendly graphic user interface for accessing the R statistical analysis and scripting language. R is often used for data analysis, statistics and graphing of very large datasets. It also has the ability to interface with the SAGA spatial analysis program, through the 69 modules shown below; users are able to use spatial data and manipulate and perform geoprocessing via special commands within this software. The more you explore and use the commands the more comfortable you will get. Initially everyone may face difficulty, however, if the basics are understood well, it eases the application of major geoprocessing tools. R studio is known for handling large datasets with efficiency and helps in optimizing results quicker in order to perform them multiple times at a time without the needs to run through the steps again and again. Also, R not only processes data but if users may want to visualize their spatial output of their data, R has a package called sp which has that feature.

To download RStudio click the link --> https://www.rstudio.com/products/rstudio/download/

There's a great introduction on the R Studio layout and configuration found on Julia Riddick's R Studio's Spatial Capabilities at http://gracilis.carleton.ca/CUOSGwiki/index.php/R_Studio%27s_Spatial_Capabilities, so take a look and become familiar with the software.

Throughout this tutorial we'll endeavour to use datasets that are pre-installed in R, saving users time in tracking down data as well as for simplicity. R studio has a great volcano dataset that represents Maunga Whau (Mt Eden). This mountain is one of about 50 volcanos in the Auckland volcanic field. This data set represents the topographic information for Maunga Whau on a 10x10 meter grid.

Photograph of Mount Eden, Source: Antonia Baker

The scripts that this tutorial will be using will be straight forward and simple, users will need to do some thinking on their own to wrap their heads around how to apply techniques elsewhere.

3D Mapping in a Myriad of Ways

Before doing any coding in R you'll need to think about what 'tools' you'll need. What packages do you need to install to run any specialized functions?

Through this tutorial we'll focus on four specific packages;

plot3D for more info look here -> https://cran.r-project.org/web/packages/plot3D/plot3D.pdf

plotly for more info look here -> https://cran.r-project.org/web/packages/plotly/plotly.pdf

rgl for more info look here -> https://cran.r-project.org/web/packages/rgl/rgl.pdf

lattice for more info look here -> https://cran.r-project.org/web/packages/lattice/lattice.pdf

The following code will install each and 'unpack' them for use.

#You can click on the packages in the package pane (bottom right) for a full description of their functions
#
#install packages utilized within this tutorial
install.packages("plotly")
install.packages("lattice")
install.packages("plot3D")
install.packages("rgl")
#impliment packages within the script
library(plotly)
library(lattice)
library(plot3D)
library(rgl)

The following examples will also illustrate how to install packages through the install window.

Using Basic Contours

The easiest way of map visualization is the use of a simple contour map.

Install the plotly package and load it using the command library(plotly).

#install the plotly package
install.packages("plotly")
#load the plotly library
library(plotly)

The following codes are examples of two basic ways of plotting contour maps using the Mount Eden dataset.

A contour line map is created through the use of the contour( ) function.

#The contour line creates a basic contour map with lines in blue
contour(volcano, main = "Topographic map of Mount Eden", col = "blue")

A color gradient representation of the contour line map is achieved through the use of the filled.contour( ) function.

#Creates a map that represents elevation change with a color gradient
filled.contour(volcano, color.palette = terrain.colors, main = "Topographic map of Mount Eden")

By utilizing the plot3D package we can take both of these contour map representations and merge them into a single map that demonstrates the elevation with a color gradient and contour lines simultaneously.

We start by installing the plot3D package and installing it

#install the plot3D package
install.packages("plot3D")
#load the plot3D library
library(plot3D)

Next we will be making use of the image2D( ) function, which will generate a color gradient for the elevation data and overlay the contour lines.

#takes a matrix as the first argument and displays it as a color gradient
#overlays contour lines and labels to identify the elevation 
image2D(Volcano, contour = list(col = "black", labcex = 1, lwd = 3, alpha = .6))

This example demonstrates the variability and flexibility of R Studio with respect to visualizing topographic data.

Using a Lattice (Wireframe)

A straight forward method of creating a 3D image on R can be done using the Lattice package. This creates powerful and elegant graphics with minimum tuning.

At its most basic level what this package will do is create a 3D scatter plot (using x, y, z) and draw wires across all the points creating a 3D image.

In the context of mapping it involves using variables like the latitude, longitude and elevation to create a 3D image, akin to a Digital Elevation Model. Within this tutorial we are going to use a stored dataset volcano to illustrate the lattice method.

Before writing any code the lattice package needs to be installed. Clicking on the Tools tab and clicking Install Packages brings up the Install Packages Window. In the blank Packages box type in the desired package and click Install.

This will take a few seconds to install but will allow provide full access to the functions within the package.

A metaphor for what you're doing by installing this package is that you've taken the lattice tool out of your 'tool shed' (R's various packages) and put them on your 'tool belt' for access.

Here's what the code looks like.

#volcano data is accessed and assigned to the variable "cano"
data("volcano")
cano <- (volcano)
#The wireframe function displays the surface of Mount Eden
wireframe(unclass(cano),shade = T, aspect=c(1,1), main = "Surface structure of Mount Eden",
          xlab = "",ylab = "",zlab = "Height", scales = list(draw = FALSE, arrows = FALSE))

The library(lattice); line is telling R that I'll be using some of the lattice tools from my 'belt' in this script.

The following two lines involve bring in the volcano dataset and putting the data into a matrix called cano.

The last line will do the work of creating the fictional scatter plot (you won't see it) and drap wires across the points to create the 3D imaged.

The parameters that are set will use an unclassified matrix called cano, shading the wires, telling it the aspect you want to look at it from, whether there's a title, labeling the 3 axes, and how you want to scales depicted.

This code will create the following image:

Try other parameters and see what happens. Maybe change shade = F or draw = TRUE? Changing the aspect will also change the image.

Using image2D and persp3D

Using the image2D and persp3D extension from R's image and persp function you can create 3D segments of maps, but you won't be using the scatter plot method as seen in utilizing a wireframe. First a 2D image will get generated using a default dataset and then you'll create a 3D segment from this 2D image.

Ensure that you have the necessary packages installed. Make sure you have plot3D and misc3d as they'll allow you to plot with multi-dimensional data. You can click on them in the package pane (bottom left) to get a better description.

Let's create our basic map. For this you'll use the hypsometry data that R has built in. Hypsometry is the measurement of land elevation relative to sea level. So along that line of thought the 2D image you'll create will illustrate areas above and below sea level.

#Create a 2 dimensional map of land elevation
image2D(Hypsometry, xlab = "Longitude", ylab = "Latitude",
        contour = list(levels = 0, col = "black", lwd = 2),
        shade = 0.1, main = "Hypsometry dataset", clab = "m")
#Determine the Area of Interest and draw an overlay
rect(-50, 10, -20, 40, lwd = 3)

Just like creating a graph you set the parameters and labels as you see fit. The Hypsometry data is first, then we set the axis labels,explained how the contours would be set, gave it a title and a simple bar explaining the differences in the colors.

This results in the following:

2D map depicting land elevation, with designated AOI

This should look a bit familiar as most elevation maps have a similar color gradient. Try googling some of the areas that are red and see if they correspond with geographical phenomenon.

From here we'll create a 3D image using a segment of the Hypsometric data (the black square). Creating a 3D image of all the globe would use up a lot of visual space and you'd be sacrificing detail for coverage.

#Create the Area of Interest within the bounds of the previously stated rectangle
#
#alonglong stores data points that are within the longitudinal bounds of the area of interest
alongLong <- which(Hypsometry$x > -50 & Hypsometry$x < -20)
#alonglat stores data points that are within the latitudinal bounds of the area of interest
alongLat <- which(Hypsometry$y > 10 & Hypsometry$y < 40)
#seabed creates a range of elevation values that will be accepted
seabed <- c(-10000,0)
#plot the surface condition based on the area of interest and elevation scale we pre defined
persp3D(z = Hypsometry$z[alongLong,alongLat], xlab = "Longitude",bty = "bl2",
        ylab = "Latitude", zlab = "Depth", main = "Looking below the sea", clab = "Depth, m",
        expand = 0.5, d = 2, phi = 20, theta = 30, resfac = 2,
        contour = list(col = "grey", side = c("zmin", "z")),
        zlim = seabed,colkey = list(side = 1, length = 0.5))

Looking at the alongLat & alongLong variables you can see that hypsometry$_ coincide with the Latitude and Longitude columns of Hypsometry.

From this you're segmenting your AOI from the larger 2D image.

The variable Seabed becomes the upper and lower limits that will be depicted on the scale.

This 3D image illustrates a ridge/trench that runs along the sea. Where else could you apply this kind of 3D image?

Creating 3D interactive maps with the RGL package

Utilizing the rgl package will allow you to create 3D maps that can be rotated 360 degrees outside of the RStudio.

RStudio doesn't have the interface to move the map's perspective on the fly (with your mouse) so it creates an external window to display the model.

A basic script that looks at the volcano dataset looks like this:

#Create a 3D model of Mount Eden
#
data(volcano)
#scale the x,y and z variables in order to exaggerate the difference between the
#highest and lowest values
z <- 3*volcano
y <- 10*(1:nrow(z))
x <- 10*(1:ncol(z))
#set the elevation range
zlim <- range(z)
zlen <- zlim[2]-zlim[1]+1
#Assign color values
colorlut <- terrain.colors(zlen, alpha = 0)
col <- colorlut[z-zlim[1]+1]
#Generate the 3D model of surface conditions
open3d()
rgl.surface(x, y, z, color = col, alpha = 1, back = "lines") 

This code utilizes the volcano dataset but exaggerates the z values in order to produce a model that better demonstrates the change in elevation.

Running this code will generate a 3D model of the selected AOI within another window.

This model can be manipulated by clicking and dragging the cursor across the screen.

From above you can see a color gradient shifting from low on the periphery and high near the peak. Looking from the underside you can see how the 'mountain' was rendered by R.

#Create a 2D model of Mount Eden
#
#Apply a different color scheme to contrast the 3D model
colorlut <- heat.colors(zlen, alpha = 1)
col <- colorlut[z-zlim[1]+1]
#Generate a 2D model of the surface conditions
rgl.surface(x, y, matrix(1, nrow(z), ncol(z)),color = col, back = "fill")

This code creates a 2D model of the same AOI which also depicts the elevation change.

The color scheme is modified to contrast the 3D model.

Running both blocks of code together allows the user to visualize the change in elevation in 2D and 3D simultaneously.

If you look closely you can see how the 2D images corresponds with the changes in elevation. The brighter areas representing areas with higher elevation than those along the edge.

You can now create 3D images that can be interacted with outside of RStudio.

Creating Composites

I've only covered a few ways of mapping in 3D but there are many more. What's really interesting is that you can bring some of the different ways of mapping together and create composites.

When you created some of the previous images/maps you also used a contour to help identify shifts in the sea level/elevation change corresponding with the changing colors.

The next big step is to bring all these different ways of mapping together to create what's called a composite.

The following code will create a composite map utilizing some of the mapping techniques you've used already.

#Create a composite map consisting of a 3D and a 2D contour representation of Mount Eden
#
#Plot the surface based on elevation data
persp3D(z = volcano, zlim = c(-60, 200), phi = 20,
        colkey = list(length = 0.2, width = 0.4, shift = 0.15,
                      cex.axis = 0.8, cex.clab = 0.85), lighting = TRUE, lphi = 90,
        clab = c("", "Height", "m"), bty = "f", plot = FALSE)
#Create a gradient variable for longitude
canoData <- volcano[-1, ] - volcano[-nrow(volcano), ]
#Create a contour map base with the same color scheme as the 3D model
image3D(z = -60, colvar = canoData/10, add = TRUE,
        colkey = list(length = 0.2, width = 0.4, shift = -0.15,
                      cex.axis = 0.8, cex.clab = 0.85),
        clab = c("", "Gradient", "m/m"), plot = FALSE)
#Add contour lines
contour3D(z = -60 + 0.01, colvar = canoData/10, add = TRUE,
          col = "black", plot = TRUE)

You can see how we're using persp3D, image3D and contour3D.

The end result should look like this.

A definition of the parameters used in the persp3D( ), image3D( ) and contour3D( ) functions are provided:

persp3D

• z: matrix containing the elevation data of the volcano dataset
• zlim: limits set to create a range of accepted values
• phi: defines the viewing angle of the plot
• length: Length of the scale
• width: Width of the scale
• shift: Offset of the scale
• cex.axis: defines the text size of the scale components
• cex.clab: defines the text size of the scale title
• lighting: When "True" natural light is imposed on the image creating depth
• lphi: lighting direction
• bty: generates a white background for the plot
• plot: when "TRUE" the matrix is displayed

plot3D

• z: height of the image drawn within the plot
• colvar: color scheme used for the image
• add: when "TRUE" this function draws its image within the same plot as the other function

Try changing some of the parameters to produce a different plot.

Conclusion

As demonstrated RStudio offers more than just stats but can be used of mapping and 3D modelling.

This tutorial is a basic introduction in how to conceptualize and do some basic 3D mapping. Users are encouraged to experiment and explore the software to get a better understanding on the full capabilities that RStudio has to offer.

Using RStudio isn't a replacement for ArcMap, as there are many capabilities that Arc can do much easier and smoother. RStudio has the advantage that it gets users to write code and thinking through steps on how to construct their map/images.