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

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)
``````

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,] . . | . .
``````