Estimate the gradient of an undefined surface

问题:

I want to estimate the gradient (slope and aspect) of an undefined surface (i.e., the function is unknown). To test my methods, here is 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 into a raster and plot using rasterVis::vectorplot. Everything looks fine. The vector field indicates the direction of the highest magnitude of change and the shading is similar to my contour plot...

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 matrix of slope values. As 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 seem really high... ( 90 degrees would 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 plot the slope and aspect, they don't seem to agree with the vectorplot graphic.

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

My thoughts:

1. Vectorplot must be smoothing the slope, but how?
2. I am fairly certain that raster::terrain is using a roving-window method to calculate slope. Perhaps the window is too small... can this be expanded?
3. Am I going about this in an inappropriate fashion? How else might I estimate the slope of an undefined surface?

回答1:

I build a 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's display this object with levelplot:

levelplot(test.rast)

and the gradient vector field with vectorplot:

vectorplot(test.rast)

If you only need the slope you get it with terrain:

slope <- terrain(test.rast, unit='degrees')

levelplot(slope, par.settings=BTCTheme())

However, if I understand you right, you really need the gradient, so you should compute both the slope and the aspect:

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

In order to understand the way vectorplot is drawing the arrows, here I show the part of its (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)
    })

Because of the structure of the original RasterLayer, the horizontal component is almost constant, so let's draw our attention on the vertical component. The next code overlays the arrows of the 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 use getValues:

saVals <- getValues(sa)

回答2:

You could also use imager functions for this:

library(imager)

# view (plot) image matrix
image(realy)

# compute gradient
gradient <- imgradient(as.cimg(realy), "xy")

# view x and y gradients
plot(gradient$x)
plot(gradient$y)

# access matrix values
mat.x <- as.matrix(gradient$x)
mat.y <- as.matrix(gradient$y)