1

I have a SpatialPolygonsDataFrame and I would like to create an adjacency matrix. However, I would like to only count two polygons as adjacent if they share more than just vertices (i.e. they share an edge, not just a point). gTouches says polygons are adjacent if they share at least one point (e.g. just one vertex).

My thought was to run:

AdjMat <- gTouches(map, byid=TRUE)

then go look at the intersections for the adjacent polygons that identifies to see if the intersection is a line or a point.

test <- gIntersection(map[x,],map[y,])
exists("test@lines")

Playing around I find "test" is "Formal class SpatialLines" if x and y share an edge, and "Formal class SpatialPoints" if they just share a vertex (and "NULL (empty)" if they aren't adjacent at all. The problem I'm having is that even when it is "Formal class SpatialLines" and "test@lines" is actually in my environment, exists("test@lines") returns FALSE. (When it is "Formal class SpatialPoints, "test@lines" is not available, but "test@coords" is, but exists("test@coords") also returns FALSE.)

How do I get an adjacency matrix that doesn't count polygons as adjacent if their shared boundary/intersection is only one-dimensional?

  • Could you provide a reproducible example for us to play around with? – MichaelChirico Apr 24 '17 at 2:43
  • @MichaelChirico I'm using [link]legis.wisconsin.gov/ltsb/gisdocs/ElectionData/GIS/… if you do map<-readOGR(dsn="WISELR_Wards_WTM8391_041712", layer="WISELR_Wards_WTM8391_041712") then AdjMat<-gTouches(map, byid=TRUE) you will for example get an adjacency matrix showing that, e.g., wards 2163 and 3220 are adjacent, but your answer below shows they only share a point in common. Your answer solves it. I just needed class instead of exists, thanks! – Scott Apr 25 '17 at 13:10
0

It's a double loop, but it gets the job done -- you were definitely on the right track. You could do some initial work minimization by using gTouches first.

Here's an example:

library(sp)
library(rgeos)
square = cbind(c(0, 1, 1, 0),
                c(0, 0, 1, 1))
ex = 
  SpatialPolygons(list(
    Polygons(list(Polygon(square)), 'a'),
    Polygons(list(Polygon(sweep(square, 2L, c(1, 0), '+'))), 'b'),
    Polygons(list(Polygon(sweep(square, 2L, c(1, 1), '+'))), 'c'),
    Polygons(list(Polygon(sweep(square, 2L, c(0, 1), '+'))), 'd'),
    Polygons(list(Polygon(sweep(square, 2L, c(5, 5), '+'))), 'e'),
    #gTouches returns FALSE for this polygon, even though it _overlaps_ a through d
    Polygons(list(Polygon(sweep(square, 2L, c(.5, .5), '+'))), 'f')))

We just loop through pairs of polygons and check if the class of their intersection is SpatialLines -- if so, they intersect by your definition, otherwise (SpatialPolygons, SpatialPoints, or NULL), they don't:

polyids = seq_along(ex)
adjMat = matrix(FALSE, ncol = length(ex), nrow = length(ex))
for (ii in polyids) {
  for (jj in setdiff(polyids, seq_len(ii))) {
    adjMat[ii, jj] = 
      class(gIntersection(ex[ii, ], ex[jj, ])) == 'SpatialLines'
  }
}
#use symmetry for the other half
adjMat[lower.tri(adjMat)] = t(adjMat)[lower.tri(adjMat)]

adjMat
#       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
# [1,] FALSE  TRUE FALSE  TRUE FALSE FALSE
# [2,]  TRUE FALSE  TRUE FALSE FALSE FALSE
# [3,] FALSE  TRUE FALSE  TRUE FALSE FALSE
# [4,]  TRUE FALSE  TRUE FALSE FALSE FALSE
# [5,] FALSE FALSE FALSE FALSE FALSE FALSE
# [6,] FALSE FALSE FALSE FALSE FALSE FALSE

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.