Я хочу оценить градиент (наклон и аспект)не определено поверхность (то есть функция неизвестна). Чтобы проверить мои методы, вот тестовые данные:

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)

Однако мне нужна матрица значений наклона. Как я понимаю vectorplot, он использует функцию 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 использует метод ровничного окна для расчета наклона. Возможно, окно слишком маленькое ... это можно расширить?Я поступаю об этом неуместно? Как еще можно оценить наклон неопределенной поверхности?