Ich möchte die Steigung (Steigung und Ausrichtung) einesnicht definiert Oberfläche (d. h. die Funktion ist unbekannt). Hier sind die Testdaten, um meine Methoden zu testen:

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)

Ich konvertiere alles in ein Raster und zeichne mitrasterVis::vectorplot. Alles sieht gut aus. Das Vektorfeld zeigt die Richtung der höchsten Änderungsgröße an und die Schattierung ähnelt meiner Konturdarstellung ...

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)

Ich benötige jedoch eine Matrix von Steigungswerten. Soweit ich Vektorplot verstehe, wird die Funktion raster :: terrain verwendet:

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

Die Steigung scheint jedoch sehr hoch zu sein ... (90 Grad wären vertikal, oder ?!)

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

und wenn ich die Neigung und den Aspekt zeichne, scheinen sie nicht mit der Vektorgrafik übereinzustimmen.

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

Meine Gedanken:

Vectorplot muss die Steigung glätten, aber wie?Das bin ich mir ziemlich sicherraster::terrain verwendet eine Roving-Window-Methode zur Berechnung der Steigung. Vielleicht ist das Fenster zu klein ... kann dies erweitert werden?Mache ich das auf unangemessene Weise? Wie könnte ich sonst die Neigung einer undefinierten Oberfläche schätzen?