5

我想估计未定义表面的梯度（斜率和方面）（即，函数未知）。为了测试我的方法，这里是测试数据：

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. 一切看起来都很好。矢量场表示变化幅度最大的方向，阴影类似于我的等高线图......

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)

但是，我需要一个斜率值矩阵。据我了解向量图，它使用 raster::terrain 函数：

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度是垂直的，对吧？！）

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

如果我绘制斜率和坡向，它们似乎与矢量图形不符。

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

我的想法：

  1. Vectorplot 必须平滑斜率，但如何？
  2. 我相当确定这raster::terrain是使用粗纱窗方法来计算斜率。可能是窗口太小了……这个可以扩大吗？
  3. 我是否以不恰当的方式处理这件事？我还能如何估计未定义曲面的斜率？
4

2 回答 2

11

RasterLayer我使用以下函数使用您的数据构建一个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：

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

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

坡

但是，如果我理解正确，您确实需要梯度，因此您应该计算斜率和坡向：

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

为了理解vectorplot绘制箭头的方式，这里我展示了其（修改后的）代码的一部分，其中计算了箭头的水平和垂直分量：

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，水平分量几乎是不变的，所以让我们把注意力集中在垂直分量上。下一个代码将向量场的箭头覆盖在这个垂直分量上。

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

dY

最后，如果您需要坡度和坡向的值，请使用 getValues：

saVals <- getValues(sa)
于 2013-05-04T07:29:42.317 回答
0

您还可以imager为此使用函数：

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)
于 2018-11-29T19:06:10.850 回答