Clump raster cells into patches, with optional neighbourhood distance.

patchify.R

patchify <- function(x, distance, p4s, givedist=TRUE) {
  # x: a binary Raster layer (0 or NA for background, and 1 for areas to be clumped)
  # distance: the neighbourhood distance. Patches that occur within this distance of 
  #           one another will be clumped. This should be in the units of the CRS given
  #           in p4s. If this is zero, the function will identify patches defined by 
  #           8-connectivity (i.e. queen's case).
  # p4s: an equal-area projection appropriate for the region. This will be used when 
  #      calculating inter-patch distances. The returned objects will be in the 
  #      original coordinate system.
  # givedist: should the distance matrix be returned? (logical). Distances are in the 
  #           units of p4s, and are shortest cartesian distances between patches.
  require(raster)
  require(rgeos)
  require(SDMTools)
  if(!is(x, 'Raster')) x <- raster(x)
  if(!is(p4s, 'CRS')) p4s <- CRS(p4s)
  if(is.na(proj4string(x))) stop(substitute(x), ' lacks a CRS.')
  x[x == 0] <- NA
  cc <- ConnCompLabel(x)
  p <- polygonizer(cc) # https://gist.github.com/johnbaums/26e8091f082f2b3dd279
  suppressWarnings(proj4string(p) <- proj4string(x))
  pproj <- spTransform(p, p4s)
  bproj <- gBuffer(pproj, width=distance)
  pproj$patch <- over(pproj, disaggregate(bproj))
  b <- spTransform(pproj, CRS(proj4string(x)))
  pout <- aggregate(b, by='patch')
  pout$patch <- as.factor(pout$patch)
  rout <- rasterize(pout, x)
  out <- list(patchrast=rout, patchpoly=pout)
  if(isTRUE(givedist)) {
    poutproj <- spTransform(pout, p4s)
    d <- gDistance(poutproj, poutproj, byid=TRUE)
    out <- c(out, list(distance=d))
  } 
  return(out)
}

Example:

First, simulate some spatially-correlated data

library(rasterVis)
library(raster)
set.seed(1)
xy <- expand.grid(x=seq(145, 150, 0.1), y=seq(-40, -35, 0.1))
Dd <- as.matrix(dist(xy))
w <- exp(-1/nrow(xy) * Dd)
Ww <- chol(w)
xy$z <- t(Ww) %*% rnorm(nrow(xy), 0, 0.1) 
coordinates(xy) <- ~x+y
r <- rasterize(xy, raster(points2grid(xy)), 'z')
r2 <- raster(r)
res(r2) <- 0.01
r2 <- resample(r, r2)
proj4string(r2) <- '+init=epsg:4283'
rthr <- r2 > quantile(r2, 0.9)
levelplot(rthr, margin=FALSE, col.regions=0:1, colorkey=FALSE)

Source the required gists

Here we load the polygonizer (see here) and patchify gists.

devtools::source_gist('26e8091f082f2b3dd279', filename='polygonizer.R')
devtools::source_gist('6d155cfca28e02b05ad5', filename='patchify.R')

Now identify patches

Here we set the neighbourhood distance to 1000 m, and specify Australian Albers (i.e., EPSG:3577) as the projection that will be used when calculating inter-patch distances. (NB: The returned raster and polygon objects will revert to the original coordinate system.)

foo <- patchify(rthr, distance=1000, '+init=epsg:3577')

The function returns the patches as a Raster object and as a SpatialPolygonsDataFrame. Below, we plot each.

# raster
spplot(foo$patchrast, col.regions=rainbow)

# poly
spplot(foo$patchpoly, col.regions=rainbow(length(foo$patchpoly)))

What happens if we increase the neighbourhood distance to 10 kilometres?

bar <- patchify(rthr, distance=5000, '+init=epsg:3577')
spplot(bar$patchpoly, col.regions=rainbow(length(bar$patchpoly)), 
    scales=list(draw=TRUE, tck=1:0))