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