Monthly Archives: April 2014

Plotting DEMs in 3D using persp(): Part 2

Quick Note: In my last post, I introduced some of the basics of the persp() function, and how to shade the facets of a persp() plot using a desired color ramp. I’m going to build off this earlier post to show how you can take a 2D DEM, and plot it in 3 dimensions. I’m going to be using the output from the script from this post, so don’t forget to run that script first before you start running the script below.

First, I’m going to extract some summary information from our final DEM, “interpDEM_mask”, particularly the minimum elevation value, summary.elev[1]. I then create a simple function, setBaseElev, that sets all NA’s in our raster to the minimum elevation value. This allows us to plot a base surface around our DEM, and give our 3D plot some spatial context:

summary.elev <- summary(interpDEM_mask)
setBaseElev <- function(x) {
  x[] <- as.numeric(summary.elev[1])
elev3D <- calc(interpDEM_mask, setBaseElev)

Next, we convert our Raster object to a regular matrix, so that we can extract the nrow and ncol, and eventually use persp() to plot it. You’ll also notice that I’ve trimmed a few rows and columns from my matrix. This isn’t necessary, but I did this because I realized that my bounding box for the raster was waaay larger than the actual DEM:

zData <-round(as.matrix(elev3D),1)
zData <- zData[c(50:540), c(50:300)]

The next lines of code are required so that we can plot our DEM using the graphics package version of persp(). If you read the persp() help page, you’ll find that both the x and y vectors need to be INCREASING in values. Since we are only plotting for the visual effect, and not for analyzing data, we can simply create vectors that have the same length as the dimensions of zData. You’ll notice that the x and y vectors are multipled by 12; this is to match the resolution of the raster we created earlier:

x = (12 * (1:nrow(zData)))    
y = (12 * (1:ncol(zData)))

We want to color our DEM based on the elevation data of our DEM (duh!), so we create a matrix, containing elevation data derived from our DEM data, which will tell R how to color the facets. This matrix MUST BE one less that the number of rows and columns in the elevation matrix (zData). We use the same function we introduced in Part 1 to calculate the elevation values that go into the zfacet matrix:

zfacet <- (zData[-1, -1] + zData[-1, -ncz] + 
           zData[-nrz, -1] + zData[-nrz, -ncz])/4

We then create a vector of colors that we will use to color the DEM facets. Here I am creating a vector with 100 colors, the first being “grey”, and the rest in hexadecimal (hex) notation from the built-in terrain.colors() color ramp. We use the cut() function to divide up the zfacet elevation data in 100 equal ranges. Each of these ranges is then assigned the corresponding color from the color vector, starting with the lowest elevation being “grey”, and so on and so forth:

nbcol <- 99
color <- c("grey", terrain.colors(nbcol))
facetcol <- cut(zfacet, nbcol+1)

Finally, we are ready to plot! I’ve added a bunch of arguements to the persp() function to make the output look pretty, but you can play around with different values to see how they change the final product:

res = persp(x, y, z = zData*2, theta = 90, phi = 30,
            col = color[facetcol],
            scale = FALSE, expand = 0.75, 
            ltheta = 75, shade = 0.75, border = NA,
            box = F, ticktype = "detailed")
3d plot of our DEM using persp().

3d plot of our DEM using persp().

Pretty nice DEM! I’m working on few new posts, so hope to get them up soon. Ciao!


Plotting DEMs in 3D using persp()

A little background… I got the idea to plot DEMs in 3D from a great R website that seemed quite popular a few years ago, but appears to have been taken down. This site helped me tremendously when I was first learning R. Anyway, one script demonstrated how to plot a DEM in 3D using persp(), but I always felt the script was more complicated that it needed to be.

Fast forward a couple years, and with the creation of the great raster package, working with rasters seems to have gotten a lot easier. Raster has a function (also called persp()) that allows you to plot raster* objects in 3D, without having to convert any data. It works well, but seems rather limited. Coloring the cells (or facets as they are referred to the help page) according to the z value seemed particularly difficult.

So, I decided to go back and see if I could simplify plotting DEMs using the graphic package version of persp(). I’ll demonstrate how I’ve plotted DEMs using persp() in an upcoming post, but I thought it’d be good to go back and review how persp() works, since I’ve had some battles getting data to plot correctly in the past.

Here are some key points to remember when using persp() to plot raster (or any type of data). I’ll then demonstrate them using a simple data matrix.

Key Point 1) The z matrix represents the “wire frame” of the plot NOT the cell centers.
Key Point 2) the x and y vectors represent the (x,y) coordinates where each value in the z matrix is plotted. Also, the x and y vectors MUST be increasing. This means that the z matrix is rotated 90 degrees counterclockwise when plotted. The length of the x vector must be the same length as the nrow(z) and the length of the y vector must equal ncol(z).
Key Point 3) When a nr x nc z matrix is plotted, it creates a (nr-1) x (nc-1) matrix of facets (cells) that can be colored using the “col” argument. This matrix is also rotated 90 degrees counterclockwise when plotted.

To demostrate these points, I will create a simple 4×4 matrix, and x and y coordinate vectors that increase:

zmat <- matrix(seq(1,16,1), nc = 4)
x <- seq(1, 4, 1)
y <- seq(1, 4, 1)
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
persp(x, y, zmat, axes= T, nticks = 4,
      ticktype = "detailed", theta = 320)

persp() plot our our simple DEM.

persp() plot our our simple DEM.

See that zmat[1,1] of our persp() surface is plotted at the origin (1,1,1), which is the lower left corner (of the x-y plane). persp() converts the x and y vectors to (x,y) coordinates using something similar to:

n = 1
coords <- list()
for(i in 1:length(x)){
  for(j in 1:length(y)){
    coords[[n]] <- c(x[i], y[j])
    n <- n+1

You can see that zmat[1,1] get plotted at (1,1,zmat[1,1]), zmat[1,2] get plotted at (1,2,zmat[1,2]), …. You’ll also notice in the image above, the wireframe created a 3×3 matrix of facets (cells) that can be colored according to a 3×3 matrix of colors. To borrow some code from for the persp() help page:

# create 3x3 matrix which averages the value of
# the 4 points bounding each facet.
nrzmat <- nrow(zmat)
nczmat <- ncol(zmat)
facetValues <- (zmat[-1, -1] + zmat[-1, -nczmat] + 
          zmat[-nrzmat, -1] + zmat[-nrzmat, -nczmat])/4

#Color Ramp
jet.colors <- colorRampPalette( c("blue", "green") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Recode facetValues into color indices
facetcol <- cut(facetValues, nbcol)

persp(x, y, zmat, axes= T, nticks = 4, ticktype = "detailed",
      theta = 320, col = color[facetcol])
persp() plot of DEM colored according to facetValue matrix.

persp() plot of DEM colored according to facetValue matrix.

The cool trick I learned from the help page example is the equation used to calculate the facetValue matrix. This will come in handy when plotting our own DEM. The equation simply calculates the average value of the four points that make up each facet, and this is what determines the facet color!

I’ve decided to break this post up into two parts. I’ll have part 2 up shortly. Cheers!