估计未定义表面的梯度 [英] 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)

我将所有内容转换为栅格并使用rasterVis::vectorplot进行绘制.一切看起来都很好.矢量场指示最大变化量的方向，阴影类似于我的等高线图...

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)

但是，我需要一个斜率值矩阵.据我了解vectorplot，它使用了raster :: terrain函数:

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

但是，坡度似乎真的很高...(90度是垂直的，对吧?！)

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

我的想法:

1. Vectorplot必须平滑坡度，但是如何?
2. 我相当确定raster::terrain正在使用流动窗口方法来计算斜率.也许窗口太小...可以扩展吗?
3. 我会以不合适的方式进行此操作吗?我还能如何估计未定义表面的斜率?

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?

推荐答案

我使用来自raster的函数用您的数据构建RasterLayer:

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

让我们用levelplot显示该对象:

levelplot(test.rast)

和带有vectorplot的梯度矢量场:

vectorplot(test.rast)

如果只需要坡度，可以通过terrain来获得它:

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

为了了解vectorplot绘制箭头的方式， 在这里，我展示了其(修改后的)代码的一部分，其中水平 计算箭头的垂直分量:

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

由于原始RasterLayer的结构， 水平分量几乎是恒定的，所以让我们画一下 注意垂直组件.下一个代码覆盖了 向量字段在此垂直分量上的箭头.

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)

最后，如果需要坡度和高宽比的值，请使用 getValues:

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

saVals <- getValues(sa)

这篇关于估计未定义表面的梯度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！