Estimate surface gradient undefined

Question

I want to evaluate the gradient (slope and aspect) of an undefined surface (i.e. the function is unknown). To test my methods, here are the test data:

require(raster); require(rasterVis) set.seed(123) x <- runif(100, min = 0, max = 1) y <- runif(100, min = 0, max = 1) e <- 0.5 * rnorm(100) test <- expand.grid(sort(x),sort(y)) names(test)<-c('X','Y') z1 <- (5 * test$X^3 + sin(3*pi*test$Y)) realy <- matrix(z1, 100, 100, byrow = F) # And a few plots for demonstration # persp(sort(x), sort(y), realy, xlab = 'X', ylab = "Y", zlab = 'Z', main = 'Real function (3d)', theta = 30, phi = 30, ticktype = "simple", cex=1.4) contour(sort(x), sort(y), realy, xlab = 'X', ylab = "Y", main = 'Real function (contours)', cex=1.4) 

I convert everything to raster and sketch using rasterVis::vectorplot . Everything looks good. The vector field indicates the direction of the highest magnitude of the change, and the shading is similar to my outline graph ...

 test.rast <- raster(t(realy), xmn = 0, xmx = 1, ymn = 0, ymx = 1, crs = CRS("+proj")) vectorplot(test.rast, par.settings=RdBuTheme, narrow = 100, reverse = T) 

However, I need a slope value matrix. Since I understand vectorplot, it uses the raster :: terrain function:

 terr.mast <- list("slope" = matrix(nrow = 100, ncol = 100, terrain(test.rast, opt = "slope", unit = "degrees", reverse = TRUE, neighbors = 8)@ data@values , byrow = T), "aspect" = matrix(nrow = 100, ncol = 100, terrain(test.rast, opt = "aspect", unit = "degrees", reverse = TRUE, neighbors = 8)@ data@values , byrow = T)) 

however, the slope seems really high ... (90 degrees will be vertical, right ?!)

 terr.mast$slope[2:6,2:6] # [,1] [,2] [,3] [,4] [,5] #[1,] 87.96546 87.96546 87.96546 87.96550 87.96551 #[2,] 84.68628 84.68628 84.68627 84.68702 84.68709 #[3,] 84.41349 84.41350 84.41349 84.41436 84.41444 #[4,] 84.71757 84.71757 84.71756 84.71830 84.71837 #[5,] 79.48740 79.48741 79.48735 79.49315 79.49367 

and if I draw the slope and aspect, they don't seem to be consistent with the vectorplot graphics.

 plot(terrain(test.rast, opt = c("slope", "aspect"), unit = "degrees", reverse = TRUE, neighbors = 8)) 

My thoughts:

• Vectorplot should smooth the slope, but how?
• I'm pretty sure raster::terrain uses the roving method to calculate the slope. Perhaps the window is too small ... can it be expanded?
• Am I thinking it inappropriately? How else can I evaluate the surface slope undefined?

Answer 1

I am RasterLayer with your data using functions from raster :

 library(raster) library(rasterVis) test.rast <- raster(ncol=100, nrow=100, xmn = 0, xmx = 1, ymn = 0, ymx = 1) xy <- xyFromCell(test.rast, 1:ncell(test.rast)) test.rast[] <- 5*xy[,1] + sin(3*pi*xy[,2]) 

Let me display this object using levelplot :

 levelplot(test.rast) 

and gradient vector field with vectorplot :

 vectorplot(test.rast) 

If you only need a slope, you will get it with terrain :

 slope <- terrain(test.rast, unit='degrees') levelplot(slope, par.settings=BTCTheme()) 

However, if I understand you correctly, you really need a gradient, so you have to calculate both the slope and the aspect:

 sa <- terrain(test.rast, opt=c('slope', 'aspect')) 

To understand how vectorplot draws arrows, here I show part of my (modified) code where the horizontal and vertical components of the arrows are calculated:

 dXY <- overlay(sa, fun=function(slope, aspect, ...){ dx <- slope*sin(aspect) ##sin due to the angular definition of aspect dy <- slope*cos(aspect) c(dx, dy) }) 

Due to the structure of the original RasterLayer , the horizontal component is almost constant, so let's focus on the vertical component. The following code overlays the arrows of a vector field over this vertical component.

 levelplot(dXY, layers=2, par.settings=RdBuTheme()) + vectorplot(test.rast, region=FALSE) 

Finally, if you need the values ​​of the slope and aspect of using getValues :

 saVals <- getValues(sa)