7

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)

library(sp)
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"))

plot(xyTransformed)

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?

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
7

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.

require(sp)
require(rgdal)
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)
     box(col="black")
       title(main="Clustering")
        legend("topleft", legend=paste("Cluster", 1:4,sep=""),
                   col=palette()[1:4], pch=rep(19,4), bg="white")
  • 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
  • 1
    @val, look at the object resulting from cutree (chc.d40), there are 4 unique values. I am not telling hclust to make 4 clusters, this is just the result when I apply a 40m distance constraint using cutree. If you want an explicit number of clusters with a distance condition then you have to implement a different approach such as kNN or a rule based model using a distance matrix. – Jeffrey Evans Apr 10 '18 at 20:55
  • 1
    @val Yes, but it would be safer to use something like: length(unique(x)) to return number of clusters. – Jeffrey Evans Apr 10 '18 at 21:57

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