localpeaks.R

library(SpatialPosition)
library(raster)
data(hospital)
# Compute Stewart potentials from known points (hospital) on a
# grid defined by its resolution
pot <- stewart(knownpts = hospital, varname = "capacity",
               typefct = "exponential", span = 750, beta = 3,
               resolution = 100, mask = paris)
# Create a raster of potentials values
ras <- rasterStewart(x = pot)


plot(ras)

## Find the maximum value within the 9-cell neighborhood of each cell
f <- function(X) max(X, na.rm = TRUE)
ww <- matrix(1, nrow=7, ncol=7) ## Weight matrix for cells in moving window
localmax <- focal(ras, fun=f, w=ww, pad=TRUE, padValue=0)
## Does each cell have the maximum value in its neighborhood?
r2 <- ras==localmax
## No max==0 allowed
r2[localmax<50] <- 0
## Get x-y coordinates of those cells that are local maxima
maxXY <- xyFromCell(r2, Which(r2==1, cells=TRUE))


plot(ras)
points(maxXY[,1], maxXY[,2], col = "red", pch = 20)