Adjacency Matrix Not Including Vertices

by Scott   Last Updated April 23, 2017 20:22 PM

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

So, any advice on how to get an adjacency matrix that doesn't count polygons as adjacent if their shared boundary/intersection is only one-dimensional? Thanks!

Tags : r adjacency


Answers 1


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
MichaelChirico
MichaelChirico
April 24, 2017 03:10 AM

Related Questions



Adjacency matrix from TopoJSON

Updated July 20, 2015 17:09 PM



Adjacency matrix from polygons in GeoJSON

Updated April 15, 2015 23:09 PM