Second-level neighbors of large polygon (neighbors' neighbors)

Question

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)

Answer 1

The key is figuring out how to kludge the neighbors object into an ngCMatrix object; turns out it's relatively straightforward, with logic basically the same as unnest functions in common data manipulation tools:

library(Matrix)

neighbors_sparse = sparseMatrix(
  i = rep(seq_along(neighbors), lengths(neighbors)),
  j = unlist(neighbors)
)


neighbors_sparse[1:5, 1:5]
# 5 x 5 sparse Matrix of class "ngCMatrix"

# [1,] . | . . .
# [2,] | . | . .
# [3,] . | . | .
# [4,] . . | . |
# [5,] . . . | .

Then we can square this sparse matrix efficiently:

neighbors_2 = neighbors_sparse %*% neighbors_sparse

Then "turn off" the "lower-order" neighbors (also efficient with the [ method for ngCMatrix):

neighbors_2[neighbors_sparse] = FALSE

NN = 1:nrow(neighbors_2)
neighbors_2[sparseMatrix(i = NN, j = NN)] = FALSE

neighbors_2[1:5, 1:5]
# 5 x 5 sparse Matrix of class "ngCMatrix"

# [1,] . . | . .
# [2,] . . . | .
# [3,] | . . . |
# [4,] . | . . .
# [5,] . . | . .