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)
print(zmat)
     [,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!

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s