14

For every one of the 208,781 Census block groups, I'd like to retrieve the FIPS IDs of all of its 1st order neighbors. I have all the TIGER boundaries downloaded and merged into a single 1GB shapefile.

I tried an ArcPython script that uses SelectLayerByLocation for BOUNDARY_TOUCHES at its core, but it takes over 1 second for each block group which is slower than I'd like. This is even after I limit the SelectLayerByLocation search to block groups in the same state. I found this script, but it also uses SelectLayerByLocation internally so it's not any faster.

The solution doesn't have to be Arc-based--I'm open to other packages, though I'm most comfortable coding with Python.

5
  • 2
    Since version 9.3, there have been tools in the Spatial Statistics toolbox to do this. Starting at 10.0, they are very efficient. I recall running a similar operation on a shapefile of comparable size (all blocks within one state) and it completed in 30 minutes, 15 of that just for disk I/O--and this was two years ago on a much slower machine. The Python source code is accessible, too.
    – whuber
    Dec 1, 2011 at 16:24
  • Which geoprocessing tool in Spatial Statistics did you use?
    – dmahr
    Dec 1, 2011 at 16:35
  • 1
    I forget its name; it is specifically for creating a table of polygon neighbor relationships. The help system encourages you to create this table before running any of the neighbor-based spatial stats tools, so that the tools don't have to recompute this information on the fly each time they run. A significant limitation, at least in the 9.x version, was that the output was in .dbf format. For a large input shapefile that won't work, in which case you either have to break the operation into pieces or hack the Python code to output in a better format.
    – whuber
    Dec 1, 2011 at 16:39
  • 1
    Is it Generate Spatial Weights Matrix?
    – dmahr
    Dec 1, 2011 at 16:42
  • Yes, that's it. The Python code fully exploits internal ArcGIS capabilities (which use spatial indexes), making the algorithm quite fast.
    – whuber
    Dec 1, 2011 at 16:49

5 Answers 5

8

For a solution avoiding ArcGIS, use pysal. You could get the weights directly from shapefiles using:

w = pysal.rook_from_shapefile("../pysal/examples/columbus.shp")

or

w = pysal.queen_from_shapefile("../pysal/examples/columbus.shp")

Head for the docs for more info.

3

If you have access to ArcGIS 10.2 for Desktop, or possibly earlier, then I think the Polygon Neighbors (Analysis) tool which:

Creates a table with statistics based on polygon contiguity (overlaps, coincident edges, or nodes).

Polygon neighbors

may make this task much easier now.

3
  • Thanks, PolyGeo. I've updated the accepted answer so the Polygon Neighbors tool gets a bit more exposure. It's definitely more robust than my manual Python-based method, though I'm not sure how the scalability with large datasets compares.
    – dmahr
    Nov 3, 2014 at 14:11
  • I am currently using 10.3, and it fails on my shapefile with ~300K census blocks.
    – soandos
    May 17, 2015 at 7:30
  • @soandos That sounds like it may be worth researching/asking as a new question.
    – PolyGeo
    May 17, 2015 at 8:18
3

Just an update. After following Whuber's advice, I found that the Generate Spatial Weights Matrix simply uses Python loops and dictionaries to determine neighbors. I reproduced the process below.

The first part loops through every vertex of every block group. It creates a dictionary with vertex coordinates as the keys and a list of block group IDs that have a vertex at that coordinate as the value. Note that this requires a topologically neat dataset, as only perfect vertex/vertex overlap will register as a neighbor relationship. Fortunately the Census Bureau's TIGER block group shapefiles are OK in this regard.

The second part loops through every vertex of every block group again. It creates a dictionary with block group IDs as the keys and that block group's neighbor IDs as the values.

# Create dictionary of vertex coordinate : [...,IDs,...]
BlockGroupVertexDictionary = {}
BlockGroupCursor = arcpy.SearchCursor(BlockGroups.shp)
BlockGroupDescription = arcpy.Describe(BlockGroups.shp)
BlockGroupShapeFieldName = BlockGroupsDescription.ShapeFieldName
#For every block group...
for BlockGroupItem in BlockGroupCursor :
    BlockGroupID = BlockGroupItem.getValue("BKGPIDFP00")
    BlockGroupFeature = BlockGroupItem.getValue(BlockGroupShapeFieldName)
    for BlockGroupPart in BlockGroupFeature:
        #For every vertex...
        for BlockGroupPoint in BlockGroupPart:
            #If it exists (and isnt empty interior hole signifier)...
            if BlockGroupPoint:
                #Create string version of coordinate
                PointText = str(BlockGroupPoint.X)+str(BlockGroupPoint.Y)
                #If coordinate is already in dictionary, append this BG's ID
                if PointText in BlockGroupVertexDictionary:
                    BlockGroupVertexDictionary[PointText].append(BlockGroupID)
                #If coordinate is not already in dictionary, create new list with this BG's ID
                else:
                    BlockGroupVertexDictionary[PointText] = [BlockGroupID]
del BlockGroupItem
del BlockGroupCursor


#Create dictionary of ID : [...,neighbors,...]
BlockGroupNeighborDictionary = {}
BlockGroupCursor = arcpy.SearchCursor(BlockGroups.shp)
BlockGroupDescription = arcpy.Describe(BlockGroups.shp)
BlockGroupShapeFieldName = BlockGroupDescription.ShapeFieldName
#For every block group
for BlockGroupItem in BlockGroupCursor:
    ListOfBlockGroupNeighbors = []
    BlockGroupID = BlockGroupItem.getValue("BKGPIDFP00")
    BlockGroupFeature = BlockGroupItem.getValue(BlockGroupShapeFieldName)
    for BlockGroupPart in BlockGroupFeature:
        #For every vertex
        for BlockGroupPoint in BlockGroupPart:
            #If it exists (and isnt interior hole signifier)...
            if BlockGroupPoint:
                #Create string version of coordinate
                PointText = str(BlockGroupPoint.X)+str(BlockGroupPoint.Y)
                if PointText in BlockGroupVertexDictionary:
                    #Get list of block groups that have this point as a vertex
                    NeighborIDList = BlockGroupVertexDictionary[PointText]
                    for NeighborID in NeighborIDList:
                        #Don't add if this BG already in list of neighbors
                        if NeighborID in ListOfBGNeighbors:
                            pass
                        #Add to list of neighbors (as long as its not itself)
                        elif NeighborID != BlockGroupID:
                            ListOfBGNeighbors.append(NeighborID)
    #Store list of neighbors in blockgroup object in dictionary
    BlockGroupNeighborDictionary[BlockGroupID] = ListOfBGNeighbors

del BlockGroupItem
del BlockGroupCursor
del BlockGroupVertexDictionary

In hindsight I realize I could have used a different method for the second part that didn't require looping through the shapefile again. But this is what I used, and it works pretty well even for 1000s of block groups at a time. I haven't tried doing it with the whole USA, but it can execute for an entire state.

2

An alternative might be to use PostgreSQL and PostGIS. I've asked a few questions on how to perform similar calculations on this site:

I found there to be a steep learning curve to figure out how the various pieces of the software fit together, but I've found it wonderful for doing calculations on large vector layers. I've run some nearest neighbor calculations on millions of polygons and it's been quick compared to ArcGIS.

1

Just some comments... the esri/ArcGIS method currently uses dictionaries to hold the information but the core calculations are done in C++ using the Polygon Neighbors Tool. This tool generates a Table that contains the contiguity information as well as optional attrs like the length of shared boundary. You can use the Generate Spatial Weights Matrix Tool if you want to store and subsequently re-use the info over and over again. You can also use this function in WeightsUtilities to generate a dictionary [random access] w/ the contiguity info:

contDict = polygonNeighborDict(inputFC, masterField, contiguityType = "ROOK")

where inputFC = any type of polygon feature class, masterField is the "unique ID" field of integers and contiguityType in {"ROOK", "QUEEN"}.

There are efforts at esri to skip the tabular aspect for Python users and go straight to an iterator which would make many use cases far faster. PySAL and the spdep package in R are fantastic alternatives [see radek's answer]. I think you are required to use shapefiles as the data format in these packages which is in tune w/ this threads input format. Not sure how they deal with overlapping polygons as well as polygons within polygons. Generate SWM as well as the function I described will count those spatial relationships as "ROOK" AND "QUEEN" Neighbors.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.