# 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 :

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))) {
class(gIntersection(ex[ii, ], ex[jj, ])) == 'SpatialLines'
}
}
#use symmetry for the other half