Measuring elevation at points on DTM surface using R

Question

Based on answers to Preventing reflective surface in rasterVis? I have been able to create natural looking digital terrain models using a combination of raster, rasterVis and rgl working on data held in ESRI ASCII grid files.

Next I want to measure the the height (elevation) of points on the surface of the digital terrain model. Ultimately I am interested in establishing the height of certain points on the hills compared to the floor of the valley nearby. Obviously, I could look at the data in the ASCII file - this is raw elevation data, after all - but it is impossible to relate the points in the plot to the data in the file. I would like to be able to visually inspect the plot, click on a point and read off the elevation, then click on another point and read that elevation.

While plotting the data in R has been very easy, I can't think of any obvious way to measure heights using an R-based solution.

I seek suggestions on how to approach this using R. Image showing what I'm after and the code used to produce it follow below. (NB I do not have permission to reproduce the original DTM data, hence the mathematical surface.)

surface

# Create an ESRI ASCII raster grid and plot it.

out.file <- "your_file_name_here.asc"

num.cols <- 50
num.rows <- 50

# Mathematical surface taken from ?persp example
x <- seq(-10, 10, length.out = num.cols)
y <- x
rotsinc <- function(x,y) {
    sinc <- function(x) { y <- sin(x)/x ; y[is.na(y)] <- 1; y }
    10 * sinc( sqrt(x^2+y^2) )
}

my.mat <- outer(x, y, rotsinc)

# Create header as per format specs
header.df <- data.frame(variable = c('ncols','nrows','xllcorner','yllcorner','cellsize','nodata_value'),
                        value = c(num.cols,num.rows,378923,4072345,2,-9999))

fileConn <- file(out.file)
write.table(header.df,
            quote = FALSE,
            row.names = FALSE,
            col.names = FALSE,
            sep = " ",
            file = fileConn)

# Append matrix to file
fileConn <- file(out.file, open = "a")
write.table(my.mat,
            row.names = FALSE,
            col.names = FALSE,
            file = fileConn)

close(fileConn)

# Now we can plot
library(raster)
library(rgdal)
library(rasterVis)

my.raster <- raster(readGDAL(out.file))
plot3D(my.raster,
       smooth = FALSE,
       specular = "black")

score 2 · Answer 1 · edited Apr 13 '17 at 12:34

plot3d and the rgl package are just ways of presenting your raster layer. This does not change the underlying basics of value extraction. What you need is just to extract the z-value of a raster-layer within a given point or location. Simply have a look at the extract() function of the raster package:

or in R:

# Just take a random point somewhere on your 
r.p <- data.frame(x=378980,y=4072400)
extract(my.raster,r.p)

And I predict your next question will be how to detect hill tops and valleys of a DEM raster layer. But please search on this website beforehand. Keyword: Focal maxima/minima. This is not a trivial matters and requires some reflections on your side on what you truly want. Whuber explained the theoretics here pretty well.

Measuring elevation at points on DTM surface using R

1 Answers1