create polygon adjacency matrix using python

Question

I am trying to create an adjacency matrix from a set of polygons. That is, I have a bunch of polygons and I want to identify which polygons have a common edge or "touch" each other. Once I find this information, I want to create an n x n matrix that indicates whether those each polygon either touches or does not touch the other polygon. Something kinda like this below.

      Poly1    Poly2   Poly3

Poly1   0        1       0
Poly2   1        0       1
Poly3   0        1       0

That is the goal anyway. I was wondering if anyone knew of a python package or such that could help me do this. I need to dynamically generate the polygons and then get their adjacency as I go, so using arcgis or qgis not the best solution.

Thanks.

Answer 1

Just figured out a way to do this in R, inspired by the link I posted in the comments (which uses outdated functionalities but the right packages); as an example, I'll use the "Southwesternmost" counties in Michigan (shapefile here--Allegan, Berrien, Cass, Kalamazoo, St. Joseph, and Van Buren counties)

So, ordering counties alphabetically, the (length>0, i.e., ignoring corners) adjacency matrix we expect is:

  A B C K S V
A 0 0 0 1 0 1
B 0 0 1 0 0 1
C 0 1 0 0 1 1
K 1 0 0 0 1 1
S 0 0 1 1 0 0
V 1 1 1 1 0 0

To get this in R, we can do the following:

library(maptools)
library(spdep)

counties<-readShapeSpatial("~/Desktop/test.shp")
adj_mat<-nb2mat(poly2nb(counties,queen=F,row.names=counties$NAME),style="B")

> adj_mat
           [,1] [,2] [,3] [,4] [,5] [,6]
Allegan       0    0    0    1    0    1
Berrien       0    0    1    0    0    1
Cass          0    1    0    0    1    1
Kalamazoo     1    0    0    0    1    1
St. Joseph    0    0    1    1    0    0
Van Buren     1    1    1    1    0    0

For some bells and whistles:

For some purposes, being a neighbor to yourself is meaningful:

diag(adj_mat)<-rep(1,nrow(adj_mat))

May also want to define column names in addition to row names:

colnames(adj_mat)<-rownames(adj_mat)

Answer 2

The simple approach is to loop through each polygon, filter the remaining polygons by the spatial extent of the polygon (plus a little buffer), and then run your adjacency test.

If your polygons are going to be of the "Simple Feature" variety, i.e. you aren't using a topological data model, you'll need to consider how you define adjacency in terms of spatial relationships rather than topological. For example you could consider two polygons adjacent if their distance is close to zero or if they "touch" or if they share a vertex ... they all will yield slightly different results for edge cases.

Here is a simple script using geodjango (should be able to adapt it to shapely or ogr easily enough). It defines adjacency as polygons that intersect after buffering by a tiny positive threshold value. It doesn't output the same format you need but can be easily modified to do so.