Raster-Karte der Lebensraumtypen zufällig ändern?

12

Ich habe ein Raster von Lebensraumtypen für ein bestimmtes Gebiet in Schottland. Ich muss zukünftige Lebensraumszenarien mit Änderungen im Lebensraum erstellen, um die Lebensfähigkeit der Population einer Vogelart zu beurteilen.

Zum Beispiel könnte es in Zukunft 10% mehr Forstwirtschaft in der Region geben. Ich möchte die aktuelle Karte ändern, indem ich die Forstwirtschaft in Blöcken einer bestimmten Größe zufällig hinzufüge. Bisher denke ich nach dem Prinzip, zufällige Punkte aus einem Raster auszuwählen, das Bereiche identifiziert, in denen die Forstwirtschaft stattfinden könnte, und die Blöcke mit der richtigen Größe mithilfe einer Art zellularer Automaten zu züchten.

Scheint dies der beste Weg zu sein, dies zu tun? Gibt es eine bessere Methode?

Wenn dies der beste verfügbare Weg ist, wie könnte ich dies vorzugsweise in R tun? (Ich sehe mir gerade die rpoints-Funktion in "spatstat" zusammen mit dem CellularAutomata-Paket an.)

Ich habe auch Zugriff auf GRASS, QGis und ArcMap 10, wenn es einfachere Möglichkeiten gibt.

Matt Geary
quelle
Hast du dir das rasterPaket schon angesehen? Es gibt viele Werkzeuge, um mit Raster-Daten (noo, rly?) Zu arbeiten.
Roman Luštrik
Danke, Roman. Ja, dies sollte mir die Werkzeuge zum Lesen und Bearbeiten meiner Basiskarte geben.
Matt Geary

Antworten:

18

Haben Sie daran gedacht, eine Markov-Kette zu verwenden ? Dies ist effektiv ein "probabilistischer zellularer Automat", wodurch die gewünschte Zufälligkeit geliefert wird. Anstatt die neue Generation in Bezug auf die lokalen Nachbarn der bestehenden Generation vorzuschreiben, gibt sie eine Wahrscheinlichkeitsverteilung für die neue Generation vor. Diese Verteilung kann beispielsweise aus zeitlichen Abfolgen von Bildern gleicher oder ähnlicher Bereiche geschätzt werden.

Dieses Modell geht intuitiv davon aus, dass eine Zelle nicht unbedingt von einer bewaldeten in eine nicht bewaldete Zelle übergeht (oder umgekehrt ), sondern dass die Wahrscheinlichkeit, dass der Übergang erfolgt, von der unmittelbaren Landbedeckung um sie herum abhängt. Es kann mit mehreren Deckungsklassen, komplexen Konfigurationen von Stadtvierteln umgehen und sogar verallgemeinert werden, um sich an die jüngste Geschichte der Landbedeckungsentwicklung zu erinnern.

Die Übergänge können mithilfe von Map Algebra-Anweisungen implementiert werden, wodurch diese Methode in jedem rasterbasierten GIS praktikabel ist, auch ohne direkten oder schnellen Zugriff auf Daten auf Zellebene. Die Verwendung von R macht es noch einfacher.

Betrachten Sie beispielsweise diese Startkonfiguration mit nur zwei Klassen, Weiß und Schwarz:

Landbedeckungsgitter

Um zu veranschaulichen, was passieren kann, habe ich ein parametrisiertes Modell erstellt (das auf keinen Daten basiert), in dem der Übergang zu Schwarz mit der Wahrscheinlichkeit 1 - q ^ k erfolgt, wobei k die durchschnittliche Anzahl schwarzer Zellen in der 3 × 3-Nachbarschaft ist (k =) 0, 1/9, 2/9, ..., 1). Wenn entweder q klein ist oder der größte Teil der Nachbarschaft bereits schwarz ist, ist die neue Zelle schwarz. Hier sind vier unabhängige Simulationen der zehnten Generation für fünf Werte von q im Bereich von 0,25 bis 0,05:

Ergebnistabelle

Offensichtlich weist dieses Modell viele der Merkmale einer Zertifizierungsstelle auf, es enthält jedoch auch einen Zufallseffekt, der für die Untersuchung alternativer Ergebnisse nützlich ist.


Code

Im Folgenden wird die Simulation in implementiert R.

#
# Make a transition from state `x` using a kernel having `k.ft` as
# its Fourier transform.
#
transition <- function(x, k.ft, q=0.1) {
  k <- zapsmall(Re(fft(k.ft * fft(x), inverse=TRUE))) / length(x)
  matrix(runif(length(k)) > q^k, nrow=nrow(k))
}
#
# Create the zeroth generation and the fft of a transition kernel.
#
n.row <- 2^7 # FFT is best with powers of 2
n.col <- 2^7
kernel <- matrix(0, nrow=n.row, ncol=n.col)
kernel[1:3, 1:3] <- 1/9
kernel.f <- fft(kernel)

set.seed(17)
x <- matrix(sample(c(0,1), n.row*n.col, replace=TRUE, prob=c(599, 1)), n.row)
#
# Prepare to run multiple simulations.
#
y.list <- list()
parameters <- c(.25, .2, .15, .1, .05)
#
# Perform and benchmark the simulations.
#
i <- 0
system.time({
  for (q in parameters) {
    y <- x
    for (generation in 1:10) {
      y <- transition(y, kernel.f, q)
    }
    y.list[[i <- i+1]] <- y
  }
})
#
# Display the results.
#    
par(mfrow=c(1,length(parameters)))
invisible(sapply(1:length(parameters), 
       function(i) image(y.list[[i]], 
                         col=c("White", "Black"),
                         main=parameters[i])))
whuber
quelle
+1 Sehr interessant. Wenn Sie historische Landbedeckungsdaten für ein bestimmtes Gebiet hätten, wäre es möglich, q und / oder k abzuleiten?
Kirk Kuykendall
2
@Kirk Ja, aber ich würde es nicht empfehlen: Das Modell, das ich zur Illustration verwendet habe, ist zu simpel. Sie können jedoch etwas Besseres ableiten: Wenn Sie die empirischen Häufigkeiten von Übergängen aus jeder aufgetretenen Nachbarschaftskonfiguration betrachten, können Sie Modelle der zukünftigen Evolution erstellen, deren Übergänge statistisch die vergangene Evolution nachahmen. Wenn die Übergangsfrequenzen räumlich homogen sind und die Zukunft sich weiterhin wie in der Vergangenheit verhält, können einige dieser Simulationen ein klares Bild davon vermitteln, wie die Zukunft aussehen wird.
whuber
Danke, das scheint genau das zu tun, was ich brauche. Wäre es möglich, den Anteil der Fläche, der sich ändert, zu begrenzen?
Matt Geary
@Matt Ja, zumindest im wahrscheinlichkeitstheoretischen Sinne. Die Theorie beschreibt, wie Markov-Ketten eine asymptotisch stabile Mischung von Anteilen jedes Zustands erreichen. Dies ist ein dynamisches Gleichgewicht: Bei jeder Generation können sich viele Zellen ändern, aber das Nettoergebnis ist, dass ihre Proportionen innerhalb des Gitters gleich bleiben (bis auf kleine zufällige Abweichungen).
Whuber
1
Ich bin ein schrecklicher R-Programmierer. Ich kann den von mir verwendeten Mathematica- Code weitergeben. Mit den Apply-Funktionen von R sollte es sich gut portieren lassen. Sie benötigen einen Kernel, eine Übergangsregel und eine Prozedur, um sie auf ein 2D 0/1-Array anzuwenden. Also: kernel = ConstantArray[1/3^2, {3,3}]für den Kernel; transitionRule [k_] := With[{q = 0.1}, Boole[RandomReal[{0, 1}] > q^k]]für die Regel; und next[a_, kernel_, f_] := Map[f, ListConvolve[kernel, a, {1, 1}, 0], {2}]um sie auf ein Array anzuwenden a . Verwenden Sie beispielsweise, um vier Generationen von Anfang an zu zeichnen ArrayPlot /@ NestList[next[#, kernel, transitionRule] &, start, 3].
whuber