Losowo zmieniasz mapę rastrową typów siedlisk?

12

Mam raster typów siedlisk dla określonego obszaru w Szkocji. Muszę stworzyć przyszłe scenariusze siedlisk ze zmianami w siedlisku, aby ocenić żywotność populacji gatunku ptaka.

Na przykład w przyszłości na tym obszarze może być o 10% więcej leśnictwa. Chciałbym zmienić obecną mapę, losowo dodając leśnictwo w blokach o określonym rozmiarze. Do tej pory myślę zgodnie z zasadami wybierania losowych punktów z rastra, który identyfikuje obszary, w których może wystąpić leśnictwo, i uprawy bloków o odpowiednich rozmiarach za pomocą jakiegoś automatu komórkowego.

Czy to wydaje się być najlepszym sposobem na to? Czy istnieje lepsza metoda?

Jeśli jest to najlepszy dostępny sposób, jak mogę to zrobić, najlepiej w R? (Obecnie patrzę na funkcję rpoints w „spatstat” wraz z pakietem CellularAutomata)

Mam również dostęp do GRASS, QGis i ArcMap 10, jeśli są prostsze sposoby na którymkolwiek z nich.

Matt Geary
źródło
Czy obejrzałeś już rasterpaczkę? Ma wiele narzędzi do pracy z danymi rastrowymi (noo, rly?).
Roman Luštrik
Dzięki, Roman. Tak, to powinno dać mi narzędzia do czytania i manipulowania moją mapą podstawową.
Matt Geary,

Odpowiedzi:

18

Czy myślałeś o użyciu łańcucha Markowa ? Jest to faktycznie „probabilistyczny automat komórkowy”, zapewniając w ten sposób pożądaną losowość. Zamiast zalecać nowe pokolenie w kategoriach lokalnych sąsiadów istniejącego pokolenia, określa rozkład prawdopodobieństwa dla nowego pokolenia. Rozkład ten można oszacować, powiedzmy, na podstawie sekwencji czasowych obrazów tego samego lub podobnego obszaru.

Intuicyjnie ten model mówi, że komórka niekoniecznie dokona przejścia z zalesionej na nieleśną (lub odwrotnie ), ale szanse, że sprawi, że przejście będzie zależeć od pokrycia terenu bezpośrednio wokół niego. Może obsługiwać wiele klas pokrycia, złożone konfiguracje dzielnic, a nawet zostać uogólniony, aby „zapamiętać” najnowszą historię ewolucji pokrycia terenu.

Przejścia można zaimplementować za pomocą instrukcji Map Algebra, co czyni tę metodę możliwą do zastosowania w dowolnym systemie GIS opartym na rastrze, nawet bez bezpośredniego lub szybkiego dostępu do danych na poziomie komórki. Używanie R sprawia, że ​​jest to jeszcze łatwiejsze.

Rozważmy na przykład tę konfigurację początkową zawierającą tylko dwie klasy, białą i czarną:

Siatka pokrycia terenu

Aby zilustrować, co może się zdarzyć, stworzyłem sparametryzowany model (nie oparty na żadnych danych), w którym przejście do czerni występuje z prawdopodobieństwem 1 - q ^ k, gdzie k jest średnią liczbą czarnych komórek w sąsiedztwie 3 na 3 (k = 0, 1/9, 2/9, ..., 1). Gdy q jest małe lub większość sąsiedztwa jest już czarna, nowa komórka będzie czarna. Oto cztery niezależne symulacje dziesiątej generacji dla pięciu wartości q w zakresie od 0,25 do 0,05:

Tabela wyników

Oczywiście model ten ma wiele cech charakterystycznych CA, ale zawiera także efekt losowy przydatny do badania alternatywnych wyników.


Kod

Poniżej zaimplementowano symulację w 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
źródło
+1 Bardzo interesujące. Gdybyś miał historyczne dane dotyczące terenu dla konkretnego obszaru, czy byłoby możliwe uzyskanie q i / lub k?
Kirk Kuykendall
2
@Kirk Tak, ale nie poleciłbym tego: model użyty do ilustracji jest zbyt uproszczony. Ale możesz uzyskać coś lepszego: patrząc na empiryczne częstotliwości przejść z każdej konfiguracji sąsiedztwa, która miała miejsce, możesz stworzyć modele przyszłej ewolucji, których przejścia statystycznie naśladują ewolucję z przeszłości. Jeśli częstotliwości przejściowe są przestrzennie jednorodne i jeśli przyszłość nadal będzie zachowywać się jak przeszłość, uruchomienie kilku z tych symulacji może dać jasny obraz tego, co przyniesie przyszłość.
whuber
Dzięki, to wydaje się robić dokładnie to, czego potrzebuję. Czy można by ustalić limit proporcji zmieniającego się obszaru?
Matt Geary,
@Matt Tak, przynajmniej w sensie probabilistycznym. Teoria opisuje, w jaki sposób łańcuchy Markowa osiągają asymptotycznie stabilną mieszaninę proporcji każdego stanu. Jest to równowaga dynamiczna: przy każdym pokoleniu wiele komórek może się zmieniać, ale wynikiem netto jest utrzymanie ich proporcji w obrębie siatki (do niewielkich odchyleń szans).
whuber
1
Jestem okropnym programistą R. Mogę udostępnić kod Mathematica, którego użyłem; z funkcjami stosującymi R, powinien dobrze się przenosić. Potrzebujesz jądra, reguły przejścia i procedury, aby zastosować je do tablicy 2D 0/1. Zatem: kernel = ConstantArray[1/3^2, {3,3}]dla jądra; transitionRule [k_] := With[{q = 0.1}, Boole[RandomReal[{0, 1}] > q^k]]dla reguły; i next[a_, kernel_, f_] := Map[f, ListConvolve[kernel, a, {1, 1}, 0], {2}]zastosować je do tablicy a . Np. Aby wykreślić cztery pokolenia od początku , użyj ArrayPlot /@ NestList[next[#, kernel, transitionRule] &, start, 3].
whuber