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!


Leave a Reply

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

You are commenting using your 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