我想估计一个梯度(斜率和坡向)不明确的表面(即函数未知)。为了测试我的方法,这里是测试数据:
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))
我的想法:
- Vectorplot 必须平滑斜率,但是如何平滑呢?
- 我相当确定
raster::terrain
使用移动窗口方法来计算坡度。也许窗户太小了……这个可以扩大吗?
- 我是否以不恰当的方式处理这件事?我还能如何估计未定义表面的斜率?