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)

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])

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!