包中的focal()
函数raster
是为这样的计算而设计的。以下代码返回所有局部最大值的坐标,包括边缘上的坐标和“高原”部分的坐标。
library(raster)
## Construct an example matrix
set.seed(444)
msize <- 10
x <- matrix(sample(seq_len(msize), msize^2, replace=TRUE), ncol=msize)
## Convert it to a raster object
r <- raster(x)
extent(r) <- extent(c(0, msize, 0, msize) + 0.5)
## Find the maximum value within the 9-cell neighborhood of each cell
f <- function(X) max(X, na.rm=TRUE)
ww <- matrix(1, nrow=3, ncol=3) ## Weight matrix for cells in moving window
localmax <- focal(r, fun=f, w=ww, pad=TRUE, padValue=NA)
## Does each cell have the maximum value in its neighborhood?
r2 <- r==localmax
## Get x-y coordinates of those cells that are local maxima
maxXY <- xyFromCell(r2, Which(r2==1, cells=TRUE))
head(maxXY)
# x y
# [1,] 8 10
# [2,] 10 10
# [3,] 3 9
# [4,] 4 9
# [5,] 1 8
# [6,] 6 8
# Visually inspect the data and the calculated local maxima
plot(r) ## Plot of heights
windows() ## Open a second plotting device
plot(r2) ## Plot showing local maxima