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]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.



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)