I have a large-ish set of polygons and am after identifying second-level neighbors of each, that is, the neighbors of the neighbors of each polygon (distinctly, i.e. the 2nd-level neighbors cannot contain self
or 1st-level neighbors).
With a smaller number of polygons this is very easy -- owing to the "traversal" property of an adjacency matrix, we can simply square the neighbors matrix and we'll have the second-level matrix (with some minor touch-ups for the "distinctness" condition), but this doesn't extend easily to a large set of polygons:
library(sp)
library(rgeos)
dim = c(150, 150)
poly = as(GridTopology(c(0, 0), c(1, 1), dim), 'SpatialPolygons')
plot(poly[sample(prod(dim), 100), ])
neighbors = gTouches(poly, byid = TRUE)
neighbors2 = neighbors %*% neighbors
Error: cannot allocate vector of size xxx Gb
(Actually it may well compute if you've got a big-RAM machine, but anyway it will be very slow)
The problem of course is that neighbors
is a huge matrix, and it's quite sparse:
format(object.size(neighbors), 'Gb')
# [1] "1.9 Gb"
mean(neighbors)
# [1] 0.0003520079
This of course is what the returnDense
argument of gTouches
is for, but I'm struggling to use the following output to get second-level neighbors:
neighbors = gTouches(poly, byid = TRUE, returnDense = FALSE)