Hi there!
I was training some ways to simulate animal (or other organisms) movements having into account habitat suitability. To do this, I used my previous eWalk model as the underlying process to simulate random or directional walks. This model is based on Brownian / Ornstein–Uhlenbeck process. You can find
more about eWalk model here!
Today, I will add one more element to this movement simulations. In this case, we will have into account the habitat or environmental preferences of the simulated species, to perform a simulation like this:
First, we will create a raster layer as a random environmental variable, for example tree cover.
library (raster)
library (dismo)
tc <- raster(nrows=100, ncols=100, xmn=0, xmx=100, ymn=0,ymx=100)
tc[] <- runif(10000, -80, 180)
tc <- focal(tc, w=matrix(1, 5, 5), mean)
tc <- focal(tc, w=matrix(1, 5, 5), mean)
plot(tc)
Second, we will define the species class. The species will be defined by their position (coordinates), and their optimum for the environmental variable.
species <- setClass("species", slots=c(x="numeric", y="numeric", opt="numeric"))
Here we will define the Red deer species as a specimen in the coordinates (50,50) and an optimum of 80 for the environmental variable (tree cover). In the same way, we will define the Egyptian mongoose as a specimen in the coordinates (50,50) and an optimum of 30 for the tree cover variable.
Red_deer <- species(x= 50, y =50, opt= 90)
Egyptian_mongoose <- species(x= 50, y =50, opt= 30)
Now, we will load the "
go" function (I do not have a name for it yet). It require a species (
sp), a raster layer with any environmental variable (
env), number of iterations (
n), a Brownian motion parameter (that is, how random is the movement of your species), a geographical optimum (the wanted destination of your species
theta_x and
theta_y), and the attraction strength or "interest" of the species to get this position (
alpha_x and
alpha_y). The syntaxis should be something like this:
path <- go (sp, env, n, sigma, theta_x, alpha_x, theta_y, alpha_y)
Here is the function to load (I will comment the function in a future post):
go <- function (sp, env, n, sigma, theta_x, alpha_x, theta_y, alpha_y) {
track <- data.frame()
track[1,1] <- sp@x
track[1,2] <- sp@y
for (step in 2:n) {
neig <- adjacent(env,
cellFromXY(env, matrix(c(track[step-1,1],
track[step-1,2]), 1,2)),
directions=8, pairs=FALSE )
options <- data.frame()
for (i in 1:length(neig)){
options[i,1]<-neig[i]
options[i,2]<- sp@opt - env[neig[i]]
}
option <- c(options[abs(na.omit(options$V2)) == min(abs(na.omit(options$V2))), 1 ],
options[abs(na.omit(options$V2)) == min(abs(na.omit(options$V2))), 1 ])
new_cell <- sample(option,1)
new_coords <- xyFromCell(env,new_cell)
lon_candidate<--9999
lat_candidate<--9999
while ( is.na(extract(env, matrix(c(lon_candidate,lat_candidate),1,2)))) {
lon_candidate <- new_coords[1]+ (sigma * rnorm(1)) + (alpha_x * ( theta_x - new_coords[1]))
lat_candidate <- new_coords[2]+ (sigma * rnorm(1)) + (alpha_y * ( theta_y - new_coords[2]))
}
track[step,1] <- lon_candidate
track[step,2] <- lat_candidate
}
return(track)
}
Well, now we can perform a simple experiment with our two specimens. We will simulate random movement of these two species having into account their environmental optimums. The "
go" function will return us the track or the path followed by each specimen (coordinates by each step).
deer_simul <- go (Red_deer, tc, 100, 2, 90, 0, 90, 0)
mongoose_simul <- go (Egyptian_mongoose, tc, 100, 2, 90, 0, 90, 0)
We can plot the paths...
plot(tc)
lines(deer_simul, lwd=1.5, col="red")
points(deer_simul, cex=0.3, col="red")
lines(mongoose_simul, lwd=1.5, col="blue")
points(mongoose_simul, cex=0.3, col="blue")
legend("topleft", legend=c("deer","mongoose"), col=c("red","blue"),
lty=c(1,1), lwd=c(2,2))
To test if each species is actually "searching" their environmental optimum, we can extract the environmental values by step for each species and plot their density distributions.
plot(density(extract(tc, deer_simul)),lwd=3, col="red", xlim=c(20,80),
ylim=c(0,max(c(density(extract(tc, deer_simul))$y,
density(extract(tc, mongoose_simul))$y))),
main="locations density distribution", xlab="tree cover")
lines(density(extract(tc, mongoose_simul)),lwd=3, col="blue")
legend("topleft", legend=c("deer","mongoose"), col=c("red","blue"),
lty=c(1,1), lwd=c(3,3))
So, we can see that the deer is actually using patches with a higher value of tree cover than the mongoose... our simulation worked! You can use the code of
this post to perform a GIF like the one above.
That's all I have to say about this for now... In the next posts we will simulate more animal movements and migrations!
PD: Let me show you a nice song about the
mongoose...
