Finding 2nd order polygon neighbors (and 3rd, 4th etc.) using ArcGIS for Desktop?

Question

I have a shapefile of contiguous polygons and for each pair of polygons, I need to know how many polygons are between them, following an neighboring/bordering approach. The output is likely to be a non-spatial table.

For example, for this group of polygons:

The output would be something like this:

and so on...

I can use ArcGIS for Desktop for that.

My idea of solution so far:

A model where, for each polygon, I run Polygon Neighbors (Analysis), merge with all the neighbors and run it again, until all there's only one polygon. However, I have little knowledge on coding/scripting, so I'm unsure on how to build this model exactly and if there's a simpler way.

Answer 1

Apologies for a solution using a different set of technologies, but this is something that is very well suited for SQL. A series of self-joins can be used to find nearest neighbors. Example below is using PostgreSQL and PostGIS using New Jersey municipality data.

WITH camden AS (
 SELECT shape 
   FROM municipalities
  WHERE mun = 'CAMDEN CITY'
)
SELECT mun AS Neighbors
  FROM municipalities
  JOIN camden ON ST_Touches(camden.shape, municipalities.shape)
 ;

Produces:

"WOODLYNNE BORO"
"COLLINGSWOOD BORO"
"HADDON TWP"
"GLOUCESTER CITY"
"OAKLYN BORO"
"PENNSAUKEN TWP"

I selected just Camden City using a CTE expression ("with" statement) and then joined that to the municipalities table using a spatial function as the join predicate. ST_Touches returns TRUE if the features touch.

This can be expanded out to all municipalities, using a self join.

SELECT a.mun AS Municipality
 , b.mun AS Neighbors
FROM municipalities a
JOIN municipalities b ON ST_Touches(a.shape, b.shape)
;

Produces the following (sampled two municipalities' results):

"HAMMONTON TOWN"; "HAMILTON TWP"
"HAMMONTON TOWN"; "MULLICA TWP"
"HAMMONTON TOWN"; "WASHINGTON TWP"
"HAMMONTON TOWN"; "SHAMONG TWP"
"HAMMONTON TOWN"; "FOLSOM BORO"
"HAMMONTON TOWN"; "WINSLOW TWP"
"HAMMONTON TOWN"; "WATERFORD TWP"
"RANDOLPH TWP";   "MINE HILL TWP"
"RANDOLPH TWP";   "PARSIPPANY-TROY HILLS TWP"
"RANDOLPH TWP";   "DOVER TOWN"
"RANDOLPH TWP";   "DENVILLE TWP"
"RANDOLPH TWP";   "ROXBURY TWP"
"RANDOLPH TWP";   "ROCKAWAY TWP"
"RANDOLPH TWP";   "CHESTER TWP"
"RANDOLPH TWP";   "MENDHAM TWP"
"RANDOLPH TWP";   "MORRIS TWP"
"RANDOLPH TWP";   "VICTORY GARDENS BORO"

You can go even more once you're working in SQL. Go the further step and calculate directions (Northwest, Southeast, etc) using ST_Azimuth:

 SELECT a.mun AS Municipality
      , b.mun AS Neighbors
      , ST_Azimuth( ST_Centroid(a.shape)
                  , ST_Centroid(b.shape) ) / (2*pi())*360 AS direction
   FROM municipalities a
   JOIN municipalities b ON ST_Touches(a.shape, b.shape)
   ;

Adds the "direction" column:

"CAMDEN CITY"; "COLLINGSWOOD BORO"; 129.847424758369
"CAMDEN CITY"; "OAKLYN BORO"; 149.434213414794
"CAMDEN CITY"; "GLOUCESTER CITY"; 190.159321103297
"CAMDEN CITY"; "PENNSAUKEN TWP"; 54.1779392174728
"CAMDEN CITY"; "HADDON TWP"; 131.032906760376

Here's a presentation I gave that talks about moving away from Desktop GIS and relying more on what your spatial database can do for you.

Hope this helps!

Update: here is using a Recursive CTE to determine neighboring polygons. I am again using NJ municipality data, clipped down to just Camden County, using Camden City as the starting point.

CREATE TABLE public.camden_neighbor AS 
WITH RECURSIVE neighbors (mun_code, shape, depth, path, cycle) AS (
 SELECT m.mun_code, m.shape, 0
      , array[m.mun_code::text], false
   FROM municipalities m
  WHERE m.mun = 'CAMDEN CITY'
  UNION ALL
  SELECT DISTINCT m.mun_code, m.shape, nm.depth+1
      , path || m.mun_code::text
      , m.mun_code = ANY(path)
   FROM municipalities m, neighbors nm
  WHERE ST_Touches(m.shape, nm.shape)
    AND NOT m.mun_code = ANY(path)
    AND nm.depth < 5
)
SELECT DISTINCT mun_code, MIN(depth) as depth
 from neighbors
 group by mun_code
ORDER BY 2, 1
;

Cribbed from the PostgreSQL documentation on Recursive CTEs, I use a depth, path, and cycle column to determine the Order and if the currently evaluated record has been seen before. From there, the results are stuffed into a new, non-spatial table. I then linked the table up to the spatial data in GIS using the "mun_code" field as a key.

Answer 2

Next code:

from math import sqrt
import itertools

def  neighborsInBetween(idx1, idx2):
    if idx1 > idx2:

        tmp = idx1
        idx1 = idx2
        idx2 = tmp

    ver = 0
    nib = 0  #neighbors in between

    while(ver == 0):

        distances = []

        for item in table:
            if item[0] == idx1 and item[1] == idx2:

                if item[3] == False:
                   nib += 1
                else:
                   ver = 1

        new_idx = []

        for item in table:
            if item[0] == idx1 and item[3] == True:
                new_idx.append(item[1])

        for item in table:
            for value in new_idx:
                if value > idx2:

                    tmp = value
                    value = idx2
                    idx2 = tmp
                if item[0] == value and item[1] == idx2:

                    distances.append(item[2])

        try:
            min_d = min ( distances )

        except ValueError:
            min_d = 0

        for item in table:
            if (min_d in item) is True:
                a = item

        idx1 = a[0]

    print nib #neighbors in between

#code starts here

layer = iface.activeLayer()

features = [feature for feature in layer.getFeatures() ]

centroids = [ feature.geometry().centroid().asPoint() for feature in features ]

n = len(centroids)

list = range(n)

table = [ [i, j, sqrt(centroids[i].sqrDist(centroids[j])), features[i].geometry().touches(features[j].geometry())] 
        for i,j in itertools.combinations(list, 2) ]

print "polygon 1   polygon2  neighborsInBetween"

for i, j in itertools.combinations(list, 2):
    print ' {:>3} {:>12}             '.format(i, j),
    neighborsInBetween(i,j)

it was run with a clone of your shapefile (see next image) at the Python Console of QGIS:

and it produced all expected results:

polygon 1   polygon2  neighborsInBetween
   0            1              0
   0            2              1
   0            3              1
   0            4              0
   0            5              0
   0            6              1
   0            7              1
   0            8              2
   0            9              3
   0           10              2
   1            2              0
   1            3              1
   1            4              0
   1            5              1
   1            6              2
   1            7              1
   1            8              2
   1            9              3
   1           10              2
   2            3              0
   2            4              0
   2            5              1
   2            6              2
   2            7              1
   2            8              1
   2            9              2
   2           10              2
   3            4              0
   3            5              1
   3            6              1
   3            7              0
   3            8              0
   3            9              1
   3           10              1
   4            5              0
   4            6              1
   4            7              0
   4            8              1
   4            9              2
   4           10              1
   5            6              0
   5            7              0
   5            8              1
   5            9              2
   5           10              1
   6            7              0
   6            8              1
   6            9              1
   6           10              0
   7            8              0
   7            9              1
   7           10              0
   8            9              0
   8           10              0
   9           10              0