6

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:

1]

The output would be something like this:

2] 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.

  • The easiest way would be to step out to the python library PySal. Another option would be to use the spdep R library. – Jeffrey Evans Feb 26 '16 at 4:55
4

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.

Camden Neighbors calculated using SQL, displayed in ArcGIS

  • Looking at your question again, you wanted to see not just neighbors, but n-order neighbors. That too can be solved using SQL. PostgreSQL has recursive CTEs, so you can use those to determine the order of the neighboring polygons. – John Reiser Jan 12 '16 at 14:41
  • Thank you for your answer John! I'm a beginner in SQL, so I'm not sure I'll get the way to do it, but I'll try my best! Thank you very much! :) – Beatriz Viseu Jan 13 '16 at 15:36
1

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:

enter image description here

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

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.