9

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?

Improve this question
2
  • 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
    Commented Jun 25, 2013 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
    Commented Jun 25, 2013 at 19:01

1 Answer 1

Reset to default
12

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")
Improve this answer
12
  • does h=200 in the cutree function definitely refer to 200m?
    – luciano
    Commented Jun 25, 2013 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
    Commented Jun 25, 2013 at 19:06
  • 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
    Commented Apr 10, 2018 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
    Commented Apr 10, 2018 at 21:57
  • 1
    @HermanToothrot this is because d represents a distance in meters and your distance units would be in decimal degrees. You can either project your data to a distance based projection or figure out the appropriate decimal degree value representing your desired distance.
    – Jeffrey Evans
    Commented Dec 31, 2019 at 19:08

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