1

I'm trying to cluster the geometries in a sf object that touch each other. For the moment I'm using the same approach that it's explained in the answer here, but the problem is that the same idea doesn't work if I consider a bigger number of geometries because the touching matrix becomes too memory-heavy and I'm not using its sparse properties.

A small and (I hope) reproducible example

# packages
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright

# download data
milan_primary_highways <- opq("Milan, Italy") %>%
  add_osm_feature(key = "highway", value = "primary") %>%
  osmdata_sf() %>%
  trim_osmdata(bb_poly = getbb("Milan, Italy", format_out = "polygon")[[1]]) %>%
  magrittr::use_series(osm_lines)

# Determine the non-sparse matrix of  touching geometries
touching_matrix <- st_touches(milan_primary_highways, sparse = FALSE)
#> although coordinates are longitude/latitude, st_touches assumes that they are planar

# Cluster the geometries
hc <- hclust(as.dist(!touching_matrix), method="single")

# Cut the dendrogram
groups <- cutree(hc, h=0.5)

# Results
table(groups)
#> groups
#>    1    2    3    4    5    6    7    8    9   10   11   12 
#> 1444   21   18    8    4    5    8    2    3    5    2    2
plot(st_geometry(milan_primary_highways), col = groups)

Is there an alternative approach that takes advantage of the sparse property of the touching matrix (that I cannot use because as.dist just accepts a numeric matrix). Do you want me to provide an example when this approach does not work? The error is simply "R cannot allocate a vector of size n GB" when I use the function hclust(as.dist(x)).

I don't know if it's important but I tried to explain the reason why I need this clustering here

  • Just to be clear are you saying your example given works just fine but it fails when you scale up to a larger dataset? – Spacedman Jan 30 at 14:09
  • And also for clarity a cluster is a set of geometries that are touching, and using h=0.5 effects that? – Spacedman Jan 30 at 14:10
  • "Just to be clear are you saying your example given works just fine but it fails when you scale up to a larger dataset?" Yes – agila Jan 30 at 14:18
  • "And also for clarity a cluster is a set of geometries that are touching, and using h=0.5 effects that?" I set h = 0.5 since all cluster of geometries that touch each other are at distance d = 0. – agila Jan 30 at 14:20
  • having read the Q you linked to (which I answered!) I think I understand it all now! – Spacedman Jan 30 at 14:20
1

Make an adjacency list:

touching_list = st_touches(milan_primary_highways)

this isn't NxN so won't grow - its a list of length N where each element is only the adjacent elements. This might help it scale up to larger problems.

Now use igraph to create a graph objects and get the connectivity:

library(igraph)
g = graph.adjlist(touching_list)
c = components(g)

The membership element of c is of length N and says which cluster it is in. It summarises the same as your example:

> table(c$membership)

   1    2    3    4    5    6    7    8    9   10   11   12 
1444   21   18    8    4    5    8    2    3    5    2    2 

So if you want to assign it to the data, do:

 milan_primary_highways$groups = c$membership
  • I just tried it on the complete dataset (formed by slightly more than 26000 geometries) and it works perfectly. Thank you very much. – agila Jan 30 at 14:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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