fig.dim <- 4 knitr::opts_chunk$set(fig.width=2*fig.dim,fig.height=fig.dim,fig.align='center') library(Matrix) library(raster) library(testthat) library(landsim)
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
We 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 NA
s 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)) expect_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 1
s,
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)) expect_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 NA
s,
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)) expect_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)
The naive method for constructing the migration matrix from a migration kernel $k$ and scale $\sigma$ has the migration rate from a cell centered at $x$ to a cell centered at $y$ equal to $k(\frac{y-x}{\sigma})$. However, if $\sigma$ is of the same size or smaller than $|x-y|$, this is not a good approximation to the corresponding continuous model. The migration rate between squares $S_x$ and $S_y$ is actually proportional to $$ \int_{S_x} \int_{S_y} k\left(\frac{u-v}{\sigma}\right) du dv . $$
Note that even if the initial kernel only depends on distance, the resulting discretized one may not. Let's check what this looks like, summing up migration rates from the cell at zero to other cells for a Gaussian kernel with a $\sigma$ equal to one-half the cell size.
naive.migration <- migration( kern="gaussian", sigma=0.5, radius=5, normalize=NULL ) center.loc <- as.integer(cellFromXY( habitat, c(0,0) )) habitat.fine <- disaggregate(habitat, fact=10) fine.locs <- xyFromCell(habitat.fine,seq_len(prod(dim(habitat.fine)))) fine.in.coarse <- cellFromXY(habitat,fine.locs) M.fine <- migration_matrix( habitat.fine, kern=naive.migration$kern, sigma=naive.migration$sigma, radius=naive.migration$radius, normalize=naive.migration$normalize, from=which(fine.in.coarse==center.loc), to=seq_len(prod(dim(habitat.fine))) ) M.aggr <- aggregate_migration( M.fine, old=habitat.fine, new=habitat, from.old=which(fine.in.coarse==center.loc), to.old=seq_len(prod(dim(habitat.fine))), from.new=center.loc, to.new=seq_len(prod(dim(habitat))) ) M.coarse <- migration_matrix( habitat, kern=naive.migration$kern, sigma=naive.migration$sigma, radius=naive.migration$radius, normalize=naive.migration$normalize, from=center.loc, to=seq_len(prod(dim(habitat))) ) plm <- function (x,...) { z <- habitat; values(z) <- x; plot(z,...) } # the difference in the matrices layout(t(1:2)) plm(M.aggr[1,],main="aggregated",zlim=range(M.aggr@x,M.coarse@x)) plm(M.coarse[1,],main="coarse",zlim=range(M.aggr@x,M.coarse@x)) plm((M.coarse-M.aggr)[1,],main="error") plm((abs(M.coarse-M.aggr)/M.aggr)[1,],main="relative error") # now, as a function of distance distvec <- values( distanceFromPoints( habitat, xyFromCell(habitat,cell=center.loc) ) ) matplot( distvec, cbind( M.aggr[1,], M.coarse[1,] ) ) matplot( distvec, cbind( M.aggr[1,], M.coarse[1,] ), log='y' )
Our strategy for producing an approximately correct migration matrix
without actually integrating over all pairs of cells
will be to construct an approximating function of distance as with the procedure above on the migration stencil
and to use this as before.
This is implemented with the function discretize kernel
:
diskern <- discretize_kernel( get_kernel(naive.migration$kern), res=res(habitat), radius=naive.migration$radius, sigma=naive.migration$sigma ) M.dk <- migration_matrix( habitat, kern=diskern, sigma=naive.migration$sigma, radius=naive.migration$radius, normalize=naive.migration$normalize, from=center.loc, to=seq_len(prod(dim(habitat))) ) matplot( distvec, cbind( M.aggr[1,], M.coarse[1,], M.dk[1,] ), xlab='distance', ylab='migration probability' ) xx <- seq(0,max(distvec),length.out=400) yy <- diskern(xx/naive.migration$sigma) * prod(attr(diskern,"res")/naive.migration$sigma) lines( xx, yy, col=3 ) legend("topright",col=1:3,lty=1:3,legend=c("aggregated","coarse","discretized coarse")) matplot( distvec, cbind( M.aggr[1,], M.coarse[1,], M.dk[1,] ), log='y', xlab='distance', ylab='migration probability' ) lines( xx, yy, col=3 )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.