0

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

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

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.