Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

This question already has an answer here:

I have some points that I need to divide into groups based on clusters of points:

x <- c(-1.482156, -1.482318, -1.482129, -1.482880, -1.485735, -1.485770, -1.485913, -1.484275, -1.485866)
y <- c(54.90083, 54.90078, 54.90077, 54.90011, 54.89936, 54.89935, 54.89935, 54.89879, 54.89902)

xy <- SpatialPoints(matrix(c(x,y), ncol=2))
proj4string(xy) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
xyTransformed <- spTransform(xy, CRS("+init=epsg:27700 +datum=WGS84"))


enter image description here

Each group should be defined using this rule: all points in the group should be within 40m of each other.

Based on the points in code above, I believe I should end up with 3-4 groups. Does anybody know how to do this in R?

share|improve this question

marked as duplicate by whuber Jun 25 '13 at 18:39

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

There are 45 posts related to clustering points in R. I am sure you will find many solutions by searching our site for them. – whuber Jun 25 '13 at 18:39
Thanks Whuber. Indeed I searched the site before I posted. i didnt look at every post, but of the post I did check, I couldn't find one that contained reproducible R code – luciano Jun 25 '13 at 19:01
up vote 4 down vote accepted

You can use a hierarchical clustering approach. By applying hclust and cutree you can derive clusters that are within a specified distance. Another way is to use the spdep package and calculate a distance matrix using dnearneigh. If you would like code for the dnearneigh approach let me know and I will post it.

d=40  # Distance threshold         

# Create example data and transform into projected coordinate system
x <- c(-1.482156, -1.482318, -1.482129, -1.482880, -1.485735, -1.485770, -1.485913, -1.484275, -1.485866)
y <- c(54.90083, 54.90078, 54.90077, 54.90011, 54.89936, 54.89935, 54.89935, 54.89879, 54.89902)

xy <- SpatialPointsDataFrame(matrix(c(x,y), ncol=2), data.frame(ID=seq(1:length(x))),
                               proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))

xy <- spTransform(xy, CRS("+init=epsg:27700 +datum=WGS84"))

chc <- hclust(dist(data.frame(rownames=rownames(xy@data), x=coordinates(xy)[,1],
              y=coordinates(xy)[,2])), method="complete")

# Distance with a 40m threshold  
chc.d40 <- cutree(chc, h=d) 

# Join results to meuse sp points
xy@data <- data.frame(xy@data, Clust=chc.d40)

# Plot results
plot(xy, col=factor(xy@data$Clust), pch=19)
        legend("topleft", legend=paste("Cluster", 1:4,sep=""),
                   col=palette()[1:4], pch=rep(19,4), bg="white")
share|improve this answer
does h=200 in the cutree function definitely refer to 200m? – luciano Jun 25 '13 at 19:03
yes, as stated in the cutree help "h is a numeric scalar or vector with heights where the tree should be cut". – Jeffrey Evans Jun 25 '13 at 19:06
if I replace coordinates(meuse) with x and y, as defined in my opening post, I need to use h=0.0009 to get 4 groups. So that is a 0.09cm distance threshold. So it doesn't like h=200 is in metres – luciano Jun 25 '13 at 19:29
I made my example specific to your data. Your transform was not working. Even though you can call spTransform from "sp" it is actually an "rgdal" function. For some reason "sp" is not adding it as a dependence so you have to explicitly add the "rgdal" package. The above code works with the 40m threshold and returns 4 clusters at this distance. – Jeffrey Evans Jun 25 '13 at 21:45

Not the answer you're looking for? Browse other questions tagged or ask your own question.