Note: if x is a raster, as.matrix(x) is equal to the transpose of matrix(values(x),nrow=nrow(x)). We want to access values directly, so we’ll stay away from as.matrix().
We’ll start with a very simple example raster:
set.seed(42)
habitat <- raster(xmn=0, xmx=5, ymn=0, ymx=5,
resolution=1,
crs="+proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
values(habitat) <- sample( c(1,2), length(habitat), replace=TRUE )
habitat.NA <- habitat
values(habitat.NA)[c(1,12)] <- NA
layout(t(1:2))
plot(habitat,main="habitat",zlim=c(0,2))
plot(habitat.NA,main="habitat, with missings",zlim=c(0,2))
And, we’ll need this:
na_to_zero <- function (x) {
x[is.na(x)] <- 0
return(x)
}Migration is just smoothing. Let’s make sure we understand what it does with boundaries, and with missing values.
We’ll use a weighting matrix that assigns weight \(1/2\) to the center and weight \(1/16\) to the surrounding 8 cells.
w <- matrix(1/16,nrow=3,ncol=3)
w[2,2] <- 1/2
w## [,1] [,2] [,3]
## [1,] 0.0625 0.0625 0.0625
## [2,] 0.0625 0.5000 0.0625
## [3,] 0.0625 0.0625 0.0625We would like to have two types of boundary: absorbing boundaries where migrants try to go, but are killed, and non-boundaries where migrants just don’t try to go. (The latter is a particular sort of reflecting boundary) We may also want to treat the internal boundaries (NA cells) and the external boundaries (border of the raster) differently.
We will often want to zero out the NAs before computation, but need to keep them around in some form to distinguish temporarily empty cells from inhospitable ones.
migration_matrix(): When precomputing a migration matrix, set absorbing boundary elements to accessible, but do not include them in from (and to). Mark reflecting (non-)boundaries as not accessible. External boundaries will be reflecting, so if they should be absorbing, you need to extend the raster and set values appropriately (see examples below).
migrate_raster(): Currently, migration with migrate_raster(), which uses raster::focal(), only implements absorbing boundaries. Reflecting boundaries using raster-based methods requires using the much-slower weighted.sum() function; see below for how to do this.
First, suppose we want migrants to exit at all boundaries. focal() allows this by padding the matrix with zeros, and using na.rm=TRUE.
hf1 <- focal( habitat, w=w, pad=TRUE, na.rm=TRUE, padValue=0 )
hf1.NA <- focal( habitat.NA, w=w, pad=TRUE, na.rm=TRUE, padValue=0 )
layout(t(1:2))
plot(hf1,main="habitat",zlim=c(0,2))
plot(hf1.NA,main="habitat, with missings",zlim=c(0,2))
Here’s the same computation:
hmat <- matrix(values(habitat),nrow=nrow(habitat),ncol=ncol(habitat))
hmat.NA <- matrix(values(habitat.NA),nrow=nrow(habitat),ncol=ncol(habitat))
expect_equivalent(as.numeric(hmat),values(habitat))
expect_equivalent(as.numeric(hmat.NA),values(habitat.NA))
hmat.pad <- rbind( 0, cbind( 0, hmat, 0 ), 0 )
hmat.NA.pad <- rbind( 0, cbind( 0, hmat.NA, 0 ), 0 )
hmf.NA <- hmf <- 0*hmat
for (dx in c(-1,0,1)) {
for (dy in c(-1,0,1)) {
hmf <- hmf + hmat.pad[dx+1+(1:nrow(hmf)),dy+1+(1:ncol(hmf))]*w[2+dx,2+dy]
add.these <- hmat.NA.pad[dx+1+(1:nrow(hmf)),dy+1+(1:ncol(hmf))]
hmf.NA <- hmf.NA + ifelse(is.na(add.these),0,add.these)*w[2+dx,2+dy]
}
}
expect_equivalent(as.numeric(hmf),values(hf1))
expect_equivalent(as.numeric(hmf.NA),values(hf1.NA))Now, let’s do it with migration_matrix. So that migrants exit at the boundaries, we need to not normalize the matrix (normalize=NULL), and compute the matrix for all cells, not just the ones that aren’t missing (from=1:length(habitat), or just construct the matrix from the habitat without missing values). Note that if normalize=NULL then migration_matrix does not normalize the kernel to sum to a constant value, but does multiply the entries by area/sigma^2, where area is the area of a cell. In this case, this factor is equal to 1.
M <- migration_matrix( habitat, kern=function(x) { ifelse(x>0,1/16,1/2) }, sigma=1, radius=1, normalize=NULL )
matrix(M[3,],nrow=nrow(habitat))## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.0000 0.0000 0 0 0
## [2,] 0.0625 0.0625 0 0 0
## [3,] 0.5000 0.0625 0 0 0
## [4,] 0.0625 0.0625 0 0 0
## [5,] 0.0000 0.0000 0 0 0expect_equal( hmf, matrix(M%*%as.numeric(hmat),nrow=nrow(hmf)) )
expect_equal( hmf.NA, matrix(M%*%na_to_zero(as.numeric(hmat.NA)),nrow=nrow(hmf)) )Now, suppose we want migrants to exit at the external boundaries, but not internal ones. To make the internal boundaries reflecting, in focal() we need to pad the matrix with zeros, and use fun=weighted.mean instead of fun=sum, with na.rm=TRUE. (This will be less efficient than with the default, fun=sum.) Note that focal returns locally fun(w*x), where x is the neighborhood of values, so we need to replace w by the matrix of 1s, and use weighted.mean rather than mean.
wm.fun <- function (x, na.rm) { weighted.mean(x,w=w,na.rm=na.rm) }
hf2 <- focal( habitat, w=(w>0), pad=TRUE, fun=wm.fun, na.rm=TRUE, padValue=0 )
hf2.NA <- focal( habitat.NA, w=(w>0), pad=TRUE, fun=wm.fun, na.rm=TRUE, padValue=0 )
layout(t(1:2))
plot(hf2,main="habitat",zlim=c(0,2))
plot(hf2.NA,main="habitat, with missings",zlim=c(0,2))
Now, let’s do it with migration_matrix. So that migrants do not exit at the internal boundaries, we need to normalize the matrix (normalize=1). So that they can exit at the external boundaries, we need to compute the migration matrix on a padded layer, and then use the function subset_migration that restricts the matrix to the part we want.
pad.extent <- extent(habitat)+c(-1,1,-1,1)*res(habitat) # extends by one cell
pad.habitat <- extend(habitat,pad.extent,value=0)
M <- migration_matrix( pad.habitat, kern=function(x) { ifelse(x>0,1/16,1/2) }, sigma=1, radius=1, normalize=1 )
M <- subset_migration( M, old=pad.habitat, new=habitat )
matrix(M[3,],nrow=nrow(habitat))## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.0000 0.0000 0 0 0
## [2,] 0.0625 0.0625 0 0 0
## [3,] 0.5000 0.0625 0 0 0
## [4,] 0.0625 0.0625 0 0 0
## [5,] 0.0000 0.0000 0 0 0expect_equal( values(hf2), as.numeric(M%*%values(habitat)) )
pad.habitat.NA <- extend(habitat.NA,pad.extent,value=0)
M.NA <- migration_matrix( pad.habitat.NA, kern=function(x) { ifelse(x>0,1/16,1/2) }, sigma=1, radius=1, normalize=1 )
M.NA <- subset_migration( M.NA, old=pad.habitat.NA, new=habitat.NA )
expect_equal( values(hf2.NA)[!is.na(values(habitat.NA))], as.numeric(M.NA%*%(values(habitat.NA)[!is.na(values(habitat.NA))])) )Now, suppose we want migration to be conservative everywhere. Now, in focal() we need to pad the matrix with NAs, and use fun=weighted.mean instead of fun=sum, with na.rm=TRUE.
wm.fun <- function (x, na.rm) { weighted.mean(x,w=w,na.rm=na.rm) }
hf3 <- focal( habitat, w=(w>0), pad=TRUE, fun=wm.fun, na.rm=TRUE, padValue=NA )
hf3.NA <- focal( habitat.NA, w=(w>0), pad=TRUE, fun=wm.fun, na.rm=TRUE, padValue=NA )
layout(t(1:2))
plot(hf3,main="habitat",zlim=c(0,2))
plot(hf3.NA,main="habitat, with missings",zlim=c(0,2))
Doing this with migration_matrix is as easy as setting normalize=1.
M <- migration_matrix( habitat, kern=function(x) { ifelse(x>0,1/16,1/2) }, sigma=1, radius=1, normalize=1 )
matrix(M[3,],nrow=nrow(habitat))## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.00000000 0.00000000 0 0 0
## [2,] 0.07692308 0.07692308 0 0 0
## [3,] 0.61538462 0.07692308 0 0 0
## [4,] 0.07692308 0.07692308 0 0 0
## [5,] 0.00000000 0.00000000 0 0 0expect_equal( values(hf3), as.numeric(M%*%values(habitat)) )
M.NA <- migration_matrix( habitat.NA, kern=function(x) { ifelse(x>0,1/16,1/2) }, sigma=1, radius=1, normalize=1 )
expect_equal( values(hf3.NA)[!is.na(values(habitat.NA))], as.numeric(M.NA%*%(values(habitat.NA)[!is.na(values(habitat.NA))])) )Here’s a simple habitat we’ll work with. Note that spatial units are in meters; the values are carrying capacities.
big.habitat <- raster(xmn=-1000, xmx=1000, ymn=-1000, ymx=1000,
resolution=100,
crs="+proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
values(big.habitat) <- sample( 100*c(1,2,5,NA), length(big.habitat), replace=TRUE )
plot(big.habitat)
We can do migration in two ways: either by matrix multiplication, or by doing smoothing operations directly on the layers.
The easiest migration matrix to create is a nearest-neighbor migration:
nn.mig <- layer_adjacency(big.habitat)
image(nn.mig)
Alternatively, consider a Gaussian kernel:
gauss.mig <- migration_matrix(big.habitat,sigma=100,radius=500,kern="gaussian")
image(gauss.mig)
Or, a two-dimensional, symmetric Cauchy kernel:
cauchy.mig <- migration_matrix(big.habitat,sigma=100,radius=800,kern="cauchy")
image(cauchy.mig)
The raster package does smoothing operations directly, with raster::focal. For instance, we can assign random values to the non-NA values in the layer and then smooth them with a Gaussian kernel:
N <- big.habitat
values(N)[!is.na(values(N))] <- rpois(sum(!is.na(values(N))),values(big.habitat)[!is.na(values(N))])
plot(N)
gauss.N <- migration_matrix(N,radius=500,sigma=100,kern="gaussian")
plot(gauss.N)