In ArcGIS 10.2, I'm trying to identify the neighbors AND neighbors of neighbors for counties in the United States.

Although I can identify which counties are neighbors to each county using either a distance table or the Polygon neighbors tool, I'm having trouble identifying the neighbors of neighbors while also excluding those neighbors of neighbors that are neighbors of the county.

Here is a picture of what I'm trying to identify. I'd like to identify the light blue features for each red feature. The dark blue counties are the neighbors of the red county; I want to identify the light blue one (i.e., neighbors OF neighbors of a feature)...

New York

  • 1
    What do you see the output being (eg., table multiple relates)? I believe that the reason that this was closed is that it is not clear as to the output. I perform this type of Nth order neighbor analysis frequently in R using the spdep package. If you provide some clarification on what type of results you would like to see and request reopening the question, I would be happy to share code. – Jeffrey Evans Jan 26 '17 at 18:34

To do this I would use ArcPy and Python to write a script from this pseudo-code.

  1. Use the Polygon Neighbors tool to write a table which has one row for each pair of polygons that are neighbours.
  2. From that table make a Python list of all the polygons that are neighbors to the polygon of interest. This could be done with list comprehension on a Search Cursor with a where clause to isolate rows with that polygon of interest.
  3. Iterate through the list of polygon neighbors from 2. to get their neighbors appended into a single list of "neighbors of neighbors". The code inside that loop will be similar to that in 2.
  4. From the list of "neigbors of neighbors" subtract the list of neighbors and then also subtract the original polygon of interest.

You should now have a list of "neigbors of neighbors", but excluding the original polygon of interest and its neighbors.


Sorry, I cant really offer you a solution in ArcGIS. If you think through the problem you will quickly realize that a one-to-one match is not possible in a flat attribute table. Each polygon will have it's own unique number of 1st, 2nd, ... order neighbors. You can do a many-to-one relate but need to generate the neighbor contingencies first.

Here is an example of creating Nth order Polygon contingencies using R and the spdep package.

columbus <- readOGR(system.file("shapes/columbus.shp", package="spData")[1])

Using the columbus polygon data we create first-order neighbors (those touching a source polygon) using the spdep::poly2nb function and second-order neighbors (those touching the polygons touching a source polygon) using the spdep::nblag function. These neighbor objects are then coerced into sparse matrices.

nb.1st <- nb2listw(poly2nb(columbus), style = "B", zero.policy = TRUE)
  W1 <- as(nb.1st, "CsparseMatrix")

nb.2nd = nblag(poly2nb(columbus), 2)
  W2 <- as(nb2listw(nb.2nd[[2]], style = "B", zero.policy = TRUE), "CsparseMatrix")

Now we can index a polygon that we are interested (i) in and use the neighbor matrices to find the associated Nth order neighbors. The vectors created bellow are tracking the row.names of the original polygon object. As such, we are creating an index based on matching row.names between the polygon and matrix object given [W>0].

i = 5  # polygon index to check neighbors 

nb.1st.idx <- colnames(W1)[which(W1[i,] > 0 )]
  nb.1st.idx <- which(row.names(columbus) %in% nb.1st.idx)

nb.2nd.idx <- colnames(W2)[which(W2[i,] > 0 )] 
  nb.2nd.idx <- which(row.names(columbus) %in% nb.2nd.idx)

We now have the corresponding 1st and 2nd order polygons for polygon 5 and can plot the results.

plot(columbus, border="grey")
  plot(columbus[nb.2nd.idx,], col="cyan", add=TRUE)
    plot(columbus[nb.1st.idx,], col="blue", add=TRUE)
      plot(columbus[i,], col="red", add=TRUE)
    title("First and Second Order Neighbors for polygon 5")       
    legend("topleft", legend=c("source polygon","1st order","2nd order"),

Keep in mind that this is only for one polygon. To do this for all polygons you would have to iterate through all polygons and store the Nth order polygon indexes in a list and then figure out how to relate them to the original data.

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.