## Simulated categorical habitat/landscape pattern (simulated habitat distribution in the landscape)
## Implementation of the Modified Random Cluster (MRC) method of Saura and Martínez-Millán (2000) Landscape Ecology 15: 661-678.
## The paper by Saura and Martínez-Millán (2000), the SIMMAP software where the MRC method was initially implemented
## and the SIMMAP user manual are available at http://www2.montes.upm.es/personales/saura/software.html#simmap.
## Read that paper (particularly the description of steps A-D there available) and the SIMMAP user manual for
## better understanding this implementation and the impact of the simulation parameters on the resultant patterns.
## The author of this R implementation is Murray Efford (University of Auckland, New Zealand), 2012-04-09,10,11,12
## Santiago Saura (Universidad Politécnica de Madrid, Spain) has made slight modifications (minor aspects only) 2012-06-06.
## The original R implementation by Murray Efford is named "randomHabitat" and is part of
## the "secr" package available at http://cran.r-project.org/web/packages/secr/index.html
## This implementation Requires 'raster' and 'igraph' packages:
## Robert J. Hijmans & Jacob van Etten (2011). raster: Geographic
## analysis and modeling with raster data. R package version 1.9-33.
## http://CRAN.R-project.org/package=raster
## Csardi G, Nepusz T: The igraph software package for complex network
## research, InterJournal, Complex Systems 1695. 2006.
## http://igraph.sf.net
## Restrictions in this implementation.
## 1) adjacency for step 'B' is based on a 1-step rook move (direction = 4)
## (i.e. 4 neighbours are considered for defining habitat patches)
## 2) Only two classes in the landscape pattern are considered (habitat vs non-habitat)
## It would be easy to adapt the R code to make it more general as related to 1) and 2).
## See the SIMMAP sofware package for generating patters without these restrictions.
## Input parameters (values to be specified below)
## 'Lx' and 'Ly' size of the rectangular pattern as the number of pixels in the vertical and horizontal directions
## 'p' initial probability (controls the degree of habitat fragmentation)
## 'A' relative area (proportion) of the landscape occupied by habitat
## 'plt' logical for whether to plot resultant patterns (both intermediate steps and final pattern)
Lx <- 200
Ly <- 200
p <- 0.4
A <- 0.7
plt <- TRUE
## This is fixed in the current code (adjacency for step 'B' is based on a 1-step rook move (direction = 4))
dirs <- 4
n<- Lx*Ly
## create rasterLayer
require(raster)
layer <- raster(nrows = Lx, ncols = Ly, xmn = 1, xmx = Ly, ymn = 1, ymx = Lx)
## step A. Percolation map generation
values(layer) <- rep(0, Lx*Ly)
values(layer)[sample.int(Lx*Ly, round(Lx*Ly*p))] <- 1 ## as close as possible
if (plt) plot(layer, useRaster = FALSE, col=c("white","black")) ## plots selected cells in black (cells with a value of in the percolation map)
## the cells with a value of 1 in the percolation map will be distributed among non-habitat (class 1) and habitat (class 2) in step C
## step B. Cluster identification (clustering of adjoining pixels)
clumped <- clump(layer, directions=dirs, gaps = FALSE)
## step C. Cluster type assignment
## ncluster is the number of clusters with a value of 1 in the percolation map
ncluster <- max(values(clumped), na.rm = TRUE)
types <- factor(c(0,1)) # non-habitat, habitat
numTypes <- as.numeric(types) # 1(non-habitat),2(habitat)
## this now assigns the clusters either to 1 (non-habitat) or 2 (habitat)
clustertype <- sample(numTypes, ncluster, replace = TRUE, prob = c(1-A,A))
values(clumped) <- clustertype[values(clumped)]
if (plt) plot(clumped, useRaster = FALSE, col=c("red","green")) ## in black the provisional non-habitat (class 1) and in green the provisional habitat (class2)
## this provisional assignment will be extended to the full image (the pixels with 0 in the inital percolation map) in step D.
## step D. Filling in image
cellsUnassigned <- (1:n)[is.na(values(clumped))]
cellsAssigned <- (1:n)[!is.na(values(clumped))]
tempadj <- adjacency(clumped, cellsUnassigned, cellsAssigned, directions = 8)
tempadj <- split(tempadj[,2], tempadj[,1])
fillinType <- function (adjcells) {
type <- values(clumped)[adjcells]
type <- factor(type)
freq <- tabulate(type)
result <- as.numeric(levels(type)[freq == max(freq)])
if (length(result)>1)
result <- sample (result, 1)
result
}
# cells with typed neighbours
filled <- sapply(tempadj, fillinType)
filledCells <- as.numeric(names(filled))
values(clumped)[filledCells] <- filled
# cells with no typed neighbours
notfilledCells <- cellsUnassigned[!(cellsUnassigned %in% filledCells)]
randomType <- sample(numTypes, length(notfilledCells), replace = TRUE, prob = c(1-A,A))
values(clumped)[notfilledCells] <- randomType
if (plt) plot(clumped, useRaster = FALSE, col=c("red","green"))