[R-sig-Geo] Need an algorithm to identify spatial clusters

[Bjarke Christensen]

Hi,

Mark Andersen:
>I need something that
>will delineate clusters of contiguous pixels with the same value (i.e.,
>burned in my simulations). Could anyone point me towards such an
>algorithm, or help me get started finding one?

It seems to me that you are basically looking for a cluster analysis
algorithm, so simply omitting 'spatial' from your search would hopefully
lead you towards the Cluster Analysis Task View
(http://cran.r-project.org/web/views/Cluster.html).

I've used agglomerative clustering (functions stats:hclust or
cluster:agnes) to accomplish a similar task before. Here is my suggestion,
assuming that your simulated forest is an equispaced grid of trees: Save
the coordinates of only those trees that have burned. Calculate the
distances between these points using dist. Run agglomerative clustering
using the minimum distance from a point to any point in a cluster as your
agglomeration method. Cut the tree at a minimum distance slightly larger
than your grid cell dimension. Voilá: you have one cluster per patch of
burned trees.

An example:

#Generate random data for an example:
set.seed(45)
pts <- expand.grid(1:10, 1:10)
names(pts) <- c("X", "Y")
pts$burned<- rbinom(nrow(pts), 1, .3)
coordinates(pts) <- c("X", "Y")
gridded(pts) <- T

#We now have a SpatialGridDataFrame, which I assume is similar to what your
data would look like.

image(pts, attr="burned", col=c("white", "gray"))

#Now calculate euclidean distances between those trees that are burned. If
that's not good enough, use lapply with sp:spDistsN1.

#Coordinates of burned trees:
coords <- coordinates(pts)[pts$burned==1,]
#Distances between them:
distances <- dist(coords)
#Agglomerative cluster analysis:
hc <- hclust(distances, method="single")
#Cluster burned trees by adjacency
patches.rook <- cutree(hc, h=1.1)
patches.queen <- cutree(hc, h=1.5)

#How many patches?
max(patches.rook)  #counts the big patch on the right as one
max(patches.queen) #counts the big patch on the right as five
#or
nlevels(factor(patches.rook))
nlevels(factor(patches.queen))

#How big is each (measured in number of trees/pixels)?
table(patches.rook)

Hope this helps,
Bjarke Christensen.