## S3 method for class 'RasterLayer':
focal(x, w=3, fun, filename='', na.rm=FALSE, pad=FALSE, padValue=NA, NAonly=FALSE, ...)
w=3
na.rm
argument (or ignore it, e.g. as one of the 'dots' arguments. For example, TRUE
, NA
will be removed from focal computations. The result will only be NA
if all focal cells are NA
. Except for some special cases (weights of 1, functions like min, max, mean), using TRUE
, additional 'virtual' rows and columns are padded to x
such that there are no edge effects. This can be useful when a function needs to have access to the central cell of the filterTRUE
, only cell values that are NA
are replaced with the computed focal valueswriteRaster
focal
uses a matrix of weights for the neighborhood of the focal cells. The default function is sum
. It is computationally much more efficient to adjust the weights-matrix than to use another function through the fun
argument. Thus while the following two statements are equivalent (if there are no NA
values), the first one is faster than the second one:
a <- focal(x, w=matrix(1/9, nc=3, nc=3))
b <- focal(x, w=3, fun=mean)
There is, however, a difference if NA
values are considered. One can use the na.rm=TRUE
option which may make sense when using a function like mean
. However, the results would be wrong when using a weights matrix.
Laplacian filter: filter=matrix(c(0,1,0,1,-4,1,0,1,0), nrow=3)
Sobel filter: filter=matrix(c(1,2,1,0,0,0,-1,-2,-1) / 4, nrow=3)
r <- raster(ncols=36, nrows=18, xmn=0)
r[] <- runif(ncell(r))
# 3x3 mean filter
r3 <- focal(r, w=matrix(1/9,nrow=3,ncol=3))
# 5x5 mean filter
r5 <- focal(r, w=matrix(1/25,nrow=5,ncol=5))
# Gaussian filter for square cells
fgauss <- function(sigma, n=5) {
m <- matrix(nc=n, nr=n)
col <- rep(1:n, n)
row <- rep(1:n, each=n)
x <- col - ceiling(n/2)
y <- row - ceiling(n/2)
# according to http://en.wikipedia.org/wiki/Gaussian_filter
m[cbind(row, col)] <- 1/(2*pi*sigma^2) * exp(-(x^2+y^2)/(2*sigma^2))
# sum of weights should add up to 1
m / sum(m)
}
gf=fgauss(1.5)
rg <- focal(r, w=gf)
# The max value for the lower-rigth corner of a 3x3 matrix around a focal cell
f = matrix(c(0,0,0,0,1,1,0,1,1), nrow=3)
f
rm <- focal(r, w=f, fun=max)
# global lon/lat data: no 'edge effect' for the columns
xmin(r) <- -180
r3g <- focal(r, w=matrix(1/9,nrow=3,ncol=3))
## focal can be used to create a cellular automaton
# Conway's Game of Life
w <- matrix(c(1,1,1,1,0,1,1,1,1), nr=3,nc=3)
gameOfLife <- function(x) {
f <- focal(x, w=w, pad=TRUE, padValue=0)
# cells with less than two or more than three live neighbours die
x[f<2 | f>3] <- 0
# cells with three live neighbours become alive
x[f==3] <- 1
x
}
# simulation function
sim <- function(x, fun, n=100, pause=0.25) {
for (i in 1:n) {
x <- fun(x)
plot(x, legend=FALSE, asp=NA, main=i)
dev.flush()
Sys.sleep(pause)
}
invisible(x)
}
# Gosper glider gun
m <- matrix(0, nc=48, nr=34)
m[c(40, 41, 74, 75, 380, 381, 382, 413, 417, 446, 452, 480,
486, 517, 549, 553, 584, 585, 586, 619, 718, 719, 720, 752,
753, 754, 785, 789, 852, 853, 857, 858, 1194, 1195, 1228, 1229)] <- 1
init <- raster(m)
# run the model
sim(init, gameOfLife, n=150, pause=0.05)
Run the code above in your browser using DataLab