Creating landscapes and simulating species occupation with MetaLandSim

February 23, 2019
By

(This article was first published on R Code – Geekcologist , and kindly contributed to R-bloggers)

MetaLandSim is an R package I’ve developed to simulate metapopulational dynamics on an habitat network. It also computes a lot of landscape connectivity metrics and simulates range expansion (check the manual if you are interested).

I’ve mentioned it here frequently, and sometimes I like to post a simple script that I used to answer a user’s question, or a that I’ve written for any other reason. This is a run through the main functions of the package.

Much more complex uses are possible, and I’ve already published a paper with some advanced uses of the package (with another one on the way – I hope!).

Let’s start!

First of all, as always, we should load the package:

# load the package
library(MetaLandSim)
## Warning: package 'MetaLandSim' was built under R version 3.5.1
## Loading required package: tcltk
## Loading required package: igraph
## Warning: package 'igraph' was built under R version 3.5.1
## ## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union

Now, to create a random landscape of 10 patches with mean area of 0.2 hectares, and with the landscape mosaic being a square with 1000 meters side.

The species mean dispersal ability is 120 meters (this is merely for graphical purposes, in order to connect the patches producing the graph).

A plot with the landscape graph is displayed graphically.

So… let’s create the landscape:

# set up the plotting window
plot.new()
par(mfrow = c(1,2))

set.seed(200)
rl1 <- rland.graph(mapsize = 1000, 
                   dist_m = 80, 
                   areaM = 0.2, 
                   areaSD = 0.115, 
                   Npatch = 10, 
                   disp = 120, 
                   plotG = TRUE)

# look at the node.characteristics
nodes <- rl1$nodes.characteristics
nodes

 

And this data frame shows the characteristics of each patch:

##           x         y      areas    radius cluster    colour nneighbour ID
## 1  533.7724 583.76503 0.18420882 24.214766       1 #FF0000FF   102.9282  1
## 2  454.3649 649.25291 0.10651184 18.412977       1 #FF0000FF   102.9282  2
## 3  589.5783 691.03989 0.04015217 11.305234       2 #FFBF00FF   120.9222  3
## 4  667.3315 839.29374 0.08311528 16.265428       3 #80FF00FF   167.4060  4
## 5  711.6001  96.50122 0.36824823 34.236976       4 #00FF40FF   149.0693  5
## 6  523.8247 235.35054 0.12806475 20.190165       5 #00FFFFFF   112.1129  6
## 7  566.7674 131.78785 0.01301267  6.435886       5 #00FFFFFF   112.1129  7
## 8  153.7271 649.28867 0.28659388 30.203587       6 #0040FFFF   300.6378  8
## 9  383.2137 307.29806 0.15232900 22.019951       7 #8000FFFF   157.9491  9
## 10 922.1776 646.32962 0.20853287 25.763943       8 #FF00BFFF   319.6587 10
# label the node points
text(x = nodes[,'x'],y = nodes[,'y'], pos = 2, offset = .5, labels = as.character(nodes[,'ID']))

# look at information associated with this object
# note the class is a landscape
class(rl1)
## [1] "landscape"
names(rl1)
## [1] "mapsize"               "minimum.distance"      "mean.area"            
## [4] "SD.area"               "number.patches"        "dispersal"            
## [7] "nodes.characteristics"

 

In order to look at the structure of the object ‘rl1’, our random landscape (the graph):

str(rl1)
## List of 7
##  $ mapsize              : num 1000
##  $ minimum.distance     : num 80
##  $ mean.area            : num 0.157
##  $ SD.area              : num 0.11
##  $ number.patches       : int 10
##  $ dispersal            : num 120
##  $ nodes.characteristics:'data.frame':   10 obs. of  8 variables:
##   ..$ x         : num [1:10] 534 454 590 667 712 ...
##   ..$ y         : num [1:10] 583.8 649.3 691 839.3 96.5 ...
##   ..$ areas     : num [1:10] 0.1842 0.1065 0.0402 0.0831 0.3682 ...
##   ..$ radius    : num [1:10] 24.2 18.4 11.3 16.3 34.2 ...
##   ..$ cluster   : int [1:10] 1 1 2 3 4 5 5 6 7 8
##   ..$ colour    : Factor w/ 8 levels "#0040FFFF","#00FF40FF",..: 6 6 8 5 2 3 3 1 4 7
##   ..$ nneighbour: num [1:10] 103 103 121 167 149 ...
##   ..$ ID        : int [1:10] 1 2 3 4 5 6 7 8 9 10
##  - attr(*, "class")= chr "landscape"

 

Obtaining the summary of the landscape:

summary_landscape(rl1)
##                                             Value
## landscape area (hectares)                 100.000
## number of patches                          10.000
## mean patch area (hectares)                  0.157
## SD patch area                               0.110
## mean distance amongst patches (meters)    391.228
## minimum distance amongst patches (meters)  60.300

 

This landscape is really a mathematical graph, where edges connect nodes (points) if they are within the dispersal distance:

help('edge.graph')
e1 <- edge.graph(rl1) 
e1
##   ID node A node B   area ndA  area ndB       XA       YA       XB
## 1  1      2      1 0.10651184 0.1842088 454.3649 649.2529 533.7724
## 2  2      7      6 0.01301267 0.1280648 566.7674 131.7879 523.8247
##         YB distance
## 1 583.7650 102.9282
## 2 235.3505 112.1129

 

And now, let’s populate the patches with a species, considering that red is empty, green is occupied:

help('species.graph')
set.seed(10)
sp_t0 <- species.graph(rl = rl1, 
                       method = "percentage", 
                       parm = 50, 
                       nsew = "none", 
                       plotG = TRUE) 

# label the node points
text(x = nodes[,'x'],y = nodes[,'y'], pos = 2, offset = 0.5, labels = as.character(nodes[,'ID']))

 

This figure shows the empty landscape on the left (rl1) and the occupied on the right (sp_t0):

fig1

Look at the information regarding this object. Note the class is metapopulation also note the addition of element ‘distance.to.neighbors’.

Let’s check the class of object ‘sp_t0’:

class(sp_t0)
## [1] "metapopulation"
names(sp_t0)
## [1] "mapsize"                "minimum.distance"      
## [3] "mean.area"              "SD.area"               
## [5] "number.patches"         "dispersal"             
## [7] "distance.to.neighbours" "nodes.characteristics"

 

And now, looking at the structure of ‘sp_t0’:

str(sp_t0)
## List of 8
##  $ mapsize               : num 1000
##  $ minimum.distance      : num 80
##  $ mean.area             : num 0.157
##  $ SD.area               : num 0.11
##  $ number.patches        : num 10
##  $ dispersal             : num 120
##  $ distance.to.neighbours:'data.frame':  10 obs. of  10 variables:
##   ..$ 1 : num [1:10] 0 103 121 288 519 ...
##   ..$ 2 : num [1:10] 103 0 142 285 610 ...
##   ..$ 3 : num [1:10] 121 142 0 167 607 ...
##   ..$ 4 : num [1:10] 288 285 167 0 744 ...
##   ..$ 5 : num [1:10] 519 610 607 744 0 ...
##   ..$ 6 : num [1:10] 349 420 460 621 234 ...
##   ..$ 7 : num [1:10] 453 530 560 715 149 ...
##   ..$ 8 : num [1:10] 386 301 438 548 785 ...
##   ..$ 9 : num [1:10] 315 349 436 603 390 ...
##   ..$ 10: num [1:10] 393 468 336 320 589 ...
##  $ nodes.characteristics :'data.frame':  10 obs. of  9 variables:
##   ..$ x         : num [1:10] 534 454 590 667 712 ...
##   ..$ y         : num [1:10] 583.8 649.3 691 839.3 96.5 ...
##   ..$ areas     : num [1:10] 0.1842 0.1065 0.0402 0.0831 0.3682 ...
##   ..$ radius    : num [1:10] 24.2 18.4 11.3 16.3 34.2 ...
##   ..$ cluster   : int [1:10] 1 1 2 3 4 5 5 6 7 8
##   ..$ colour    : Factor w/ 8 levels "#0040FFFF","#00FF40FF",..: 6 6 8 5 2 3 3 1 4 7
##   ..$ nneighbour: num [1:10] 103 103 121 167 149 ...
##   ..$ ID        : int [1:10] 1 2 3 4 5 6 7 8 9 10
##   ..$ species   : num [1:10] 1 0 1 1 1 1 0 0 0 0
##  - attr(*, "class")= chr "metapopulation"

 

Look at some elements like distance.to.neighbors:

head(sp_t0$distance.to.neighbours, n = 3)
##          1        2        3        4        5        6        7        8
## 1   0.0000 102.9282 120.9222 288.3278 518.6990 348.5565 453.1799 385.6524
## 2 102.9282   0.0000 141.5232 285.4300 609.6756 419.6902 529.5323 300.6378
## 3 120.9222 141.5232   0.0000 167.4060 606.9312 460.4089 559.7171 437.8463
##          9       10
## 1 314.8046 393.4118
## 2 349.2787 467.8218
## 3 435.7111 335.5909

Note the species column in nodes.characteristics, it contains the information on the presence and absence of the focal species in the habitat network:

sp_t0$nodes.characteristics
##           x         y      areas    radius cluster    colour nneighbour ID
## 1  533.7724 583.76503 0.18420882 24.214766       1 #FF0000FF   102.9282  1
## 2  454.3649 649.25291 0.10651184 18.412977       1 #FF0000FF   102.9282  2
## 3  589.5783 691.03989 0.04015217 11.305234       2 #FFBF00FF   120.9222  3
## 4  667.3315 839.29374 0.08311528 16.265428       3 #80FF00FF   167.4060  4
## 5  711.6001  96.50122 0.36824823 34.236976       4 #00FF40FF   149.0693  5
## 6  523.8247 235.35054 0.12806475 20.190165       5 #00FFFFFF   112.1129  6
## 7  566.7674 131.78785 0.01301267  6.435886       5 #00FFFFFF   112.1129  7
## 8  153.7271 649.28867 0.28659388 30.203587       6 #0040FFFF   300.6378  8
## 9  383.2137 307.29806 0.15232900 22.019951       7 #8000FFFF   157.9491  9
## 10 922.1776 646.32962 0.20853287 25.763943       8 #FF00BFFF   319.6587 10
##    species
## 1        1
## 2        0
## 3        1
## 4        1
## 5        1
## 6        1
## 7        0
## 8        0
## 9        0
## 10       0

 

Summarizing the object ‘sp_t0’:

summary_metapopulation(sp_t0)
##                                             Value
## landscape area (hectares)                 100.000
## number of patches                          10.000
## mean patch area (hectares)                  0.157
## SD patch area                               0.110
## mean distance amongst patches (meters)    391.228
## minimum distance amongst patches (meters)  60.300
## species occurrence - snapshot 1            50.000

 

And now, finally, we can simulate the occupation on the next time step (function ”spom’), after defining the parameters setting up the options for the simulation:

# ultimately we need to provide the 'rules' by which patch dynamics occur; notice there several options
help("spom")

#  establish parameters that control changes in occupancy through time; here use the built-in dataset that contains the key parameters
help('param1')
data(param1)

# view this object; note the four parameters named alpha, x, y, e;  these will be used to control our dynamics 
View(param1)

# to simulate dynamics through time, use spom function
# spom stands for stochastic patch occupancy model
help('spom')

sp_t1 <- spom(sp_t0, kern = "op1", 
              conn = "op1", 
              colnz = "op1", 
              ext = "op1",
              param_df = param1, 
              beta1 = NULL, 
              b = 1, 
              c1 = NULL, 
              c2 = NULL, 
              z = NULL, 
              R = NULL)

class(sp_t1)
## [1] "list"

 

The dynamics are contained in the data frame ‘nodes.characteristic’. It depicts the species presence at time step t (species), time step t+1 (species2) and turnover:

head(sp_t1$nodes.characteristics)
##          x         y      areas   radius cluster    colour nneighbour ID
## 1 533.7724 583.76503 0.18420882 24.21477       1 #FF0000FF   102.9282  1
## 2 454.3649 649.25291 0.10651184 18.41298       1 #FF0000FF   102.9282  2
## 3 589.5783 691.03989 0.04015217 11.30523       2 #FFBF00FF   120.9222  3
## 4 667.3315 839.29374 0.08311528 16.26543       3 #80FF00FF   167.4060  4
## 5 711.6001  96.50122 0.36824823 34.23698       4 #00FF40FF   149.0693  5
## 6 523.8247 235.35054 0.12806475 20.19017       5 #00FFFFFF   112.1129  6
##   species species2 turn
## 1       1        1    0
## 2       0        0    0
## 3       1        1    0
## 4       1        1    0
## 5       1        1    0
## 6       1        1    0

 

To simulate the dynamics across multiple steps and get the results use the ‘iterate.graph’ function, which runs the full simulation:

it1 <- iterate.graph(iter = 2, 
                     mapsize = 10000, 
                     dist_m = 10, 
                     areaM = 0.05, 
                     areaSD = 0.02, 
                     Npatch = 250, 
                     disp = 800,
                     span = 100, 
                     par1 = "hab", 
                     par2 = 2, 
                     par3 = NULL, 
                     par4 = NULL, 
                     par5 = NULL, 
                     method = "percentage", 
                     parm = 50, 
                     nsew = "none", 
                     succ = "none", 
                     param_df = param1,
                     kern = "op1", 
                     conn = "op1", 
                     colnz = "op1", 
                     ext = "op1", 
                     beta1 = NULL, 
                     b = 1, 
                     c1 = NULL, 
                     c2 = NULL, 
                     z = NULL, 
                     R = NULL, 
                     graph = TRUE)
## Completed iteration 1  of  2 
## Completed iteration 2  of  2

 

Which will open the web browser with the following graphs:

2019-02-22 (2).png

And that’s all for now!

Thanks!

 

To leave a comment for the author, please follow the link and comment on their blog: R Code – Geekcologist .

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)