Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I need a graph where each node is a US county and each edge represents a border shared between counties. I do not particularly care about the absolute position or shape of each county (though that would be a plus).

Where can I find that information (for free)? Army Corp of Engineers? US Geo Survey?

Ideally, I could just get a csv list something like

  • FL-Polk, FL-Lake
  • FL-Polk, FL-Orange
  • FL-Polk, FL-Osceola
  • etc for every county, including "kitty-corner" borders, borders across water, or across state borders

I'm a programmer so expect I can handle any standard exchange format. (That could just be hubris though). (I don't have any cool GIS apps. Just humble scripting in perl and python.)

share|improve this question
    
If you have arcgis, this could be derived using ITopologyGraph for a topology built using counties. –  Kirk Kuykendall Sep 13 '11 at 2:56
add comment

5 Answers 5

You could find a topological representation of county boundaries, such as the CTA boundary files (another source would be the OpenStreetMap, but that's not as complete) and then pick out all distinct pairs of county IDs on opposing sides of boundary lines - e.g. in the CTA files, there's these two fields in the link entity (as described in the documentation):

 (10)  county FIPS left    I6      48-53   
 (11)  county FIPS right   I6      54-59

EDIT: After request, here's the details on how to get the data:

  1. To get all county codes, go to the page linked above, click on the 'List of county codes'.

  2. To read the adjacent IDs, there's two options:

    • download the shapefile ("GO" button next to Download SCUL), open the zip and read the .dbf file (there's loads of programs that open it, e.g. Excel, and the file is basically fixed-width-row ASCII table, so you should be able to extract the last two columns easily)

    • or download the native format (i'd first try to see if the state boundaries contains the counties as well) open the .llr file, split it by lines (line ends with CR+LF = 0D 0A) and extract the ids from each line: left FIPS is in chars 48-53; right FIPS is in chars 54-59 (index of the first character in a line is 1)

share|improve this answer
    
+1 Using the topological information is far more efficient than the brute-force comparison of each pair of counties! –  whuber Sep 13 '11 at 13:35
2  
Will this handle the kitty-corner case? –  Kirk Kuykendall Sep 13 '11 at 13:43
    
OSM probably wouldn't be a good data source at this point, because many of the county boundaries have likely been edited and may no longer be topological. I'm surprised no one's mentioned using TIGER (after all, the TI stands for Topologically Integrated!) –  neuhausr Sep 13 '11 at 13:57
    
Couple other things: between Census 2000 and 2010 there were 4 changes in the counties (+2 Skagway and Wrangell, AK; +1 Broomfield, CO; -1 Clifton Forge, VA ). Looks like the CTA file is from 2003, and doesn't appear to include these changes. –  neuhausr Sep 13 '11 at 15:34
    
@Kirk No, it won't handle the kitty-corner. To find those, you'd have to keep track of counties touching a node, which will consume some memory, but it should be pretty straightforward. –  mkadunc Sep 13 '11 at 23:11
show 4 more comments

Since you mention you are a programmer, here's some code that works with arcgis 10. Update: If you don't feel like programming, I've posted a zipped shapefile of the graph here.

public void TestGetNeighbors()
{
    var fLayer = ArcMap.Document.FocusMap.get_Layer(0) as IFeatureLayer;
    var dict = GetNeighborsByName((ITopologyClass)fLayer.FeatureClass, "{0} {1}", "Name", "State_Name");
    foreach (KeyValuePair<string, List<string>> kvp in dict)
    {
        Debug.WriteLine(kvp.Key);
        foreach (string neighbor in kvp.Value)
            Debug.WriteLine("\t" + neighbor);
    }
}
private Dictionary<string, List<string>> GetNeighborsByName(ITopologyClass topoClass, string format, params object[] fldNames)
{
    var fc = topoClass as IFeatureClass;
    if (topoClass.Topology.Cache.BuildExtent == null || topoClass.Topology.Cache.BuildExtent.IsEmpty)
    {
        Debug.WriteLine("building ...");
        topoClass.Topology.Cache.Build(((IGeoDataset)fc).Extent, false);
    }

    // get neighbors by oid
    var oidDict = GetNeighborsByOid(topoClass);

    // use the full names to build the output dictionary
    var nameDict = GetFullNames(fc, format, fldNames);
    var outDict = new Dictionary<string, List<string>>();
    foreach (KeyValuePair<int, List<int>> kvp in oidDict)
    {
        var list = new List<string>();
        foreach (int oid in kvp.Value)
            list.Add(nameDict[oid]);
        outDict.Add(nameDict[kvp.Key], list);
    }
    return outDict;
}

private Dictionary<int, List<int>> GetNeighborsByOid(ITopologyClass topoClass)
{            
    var outDict = new Dictionary<int, List<int>>();
    IFeatureCursor fCur = ((IFeatureClass)topoClass).Search(null, false);
    try
    {
        IFeature feat = null;
        while ((feat = fCur.NextFeature()) != null)
        {
            var neighbors = GetNeighboringFeatureOids(topoClass.Topology.Cache,(IFeatureClass)topoClass , feat.OID);
            outDict.Add(feat.OID, neighbors);
        }
    }
    catch
    {
        throw;
    }
    finally
    {
        if(fCur != null)
            System.Runtime.InteropServices.Marshal.FinalReleaseComObject(fCur);
    }
    return outDict;
}

private List<int> GetNeighboringFeatureOids(ITopologyGraph topoGraph, IFeatureClass fc, int oid)
{
    var outList = new List<int>();
    var enumEdge = topoGraph.GetParentEdges(fc, oid);
    for (int i = 0; i < enumEdge.Count; i++)
    {
        var edge = enumEdge.Next();
        outList.AddRange(GetParentOids(edge.get_LeftParents(true),fc));
        outList.AddRange(GetParentOids(edge.get_RightParents(true),fc));               
    }
    // handle the kitty-corner case
    var enumNode = topoGraph.GetParentNodes(fc, oid);
    for(int i=0; i<enumNode.Count;i++)
    {
        var node = enumNode.Next();
        var list = GetParentOids(node.Parents, fc);
        foreach (int neighborOid in list)
        {
            if (!outList.Contains(neighborOid))
                outList.Add(neighborOid); // kitty corner
        }
    }
    return Reduce(outList,oid);
}

private List<int> Reduce(List<int> inList, int omit)
{
    var outList = new List<int>();
    foreach (int i in inList)
    {
        if (!(outList.Contains(i) || omit == i))
            outList.Add(i);
    }
    return outList;
}

private List<int> GetParentOids(IEnumTopologyParent enumParent, IFeatureClass fc)
{
    var outList = new List<int>();
    for (int i = 0; i < enumParent.Count; i++)
    {
        var parent = enumParent.Next();
        if (parent.m_pFC == fc)
            outList.Add(parent.m_FID);
    }
    return outList;
}

private Dictionary<int, string> GetFullNames(IFeatureClass fc, string format, params object[] fldNames)
{
    var fldList = GetFieldIndexes(fc, fldNames);
    var outDict = new Dictionary<int, string>();
    IFeatureCursor fCur = fc.Search(null, false);
    IFeature feat = null;
    try
    {
        while ((feat = fCur.NextFeature()) != null)
        {
            var valList = new List<string>();
            foreach (int idx in fldList)
            {
                string val = feat.get_Value(idx) is DBNull ? "" : feat.get_Value(idx).ToString();
                valList.Add(val);
            }
            outDict.Add(feat.OID, String.Format(format, valList.ToArray()));
        }
    }
    catch
    {
        throw;
    }
    finally
    {
        if (fCur != null)
            System.Runtime.InteropServices.Marshal.FinalReleaseComObject(fCur);
    }
    return outDict;
}

private static List<int> GetFieldIndexes(IFeatureClass fc, object[] fldNames)
{
    var fldList = new List<int>();
    foreach (string fldName in fldNames)
    {
        int idx = fc.FindField(fldName);
        if (idx == -1)
            throw new Exception(string.Format("field {0} not found on {1}", fldName, fc.AliasName));
        fldList.Add(idx);
    }
    return fldList;
}

It produces output like this:

Yukon-Koyukuk Alaska
    North Slope Alaska
    Northwest Arctic Alaska
    Southeast Fairbanks Alaska
    Nome Alaska
    Fairbanks North Star Alaska
    Denali Alaska
    Wade Hampton Alaska
    Matanuska-Susitna Alaska
    Bethel Alaska
Southeast Fairbanks Alaska
    Yukon-Koyukuk Alaska
    Fairbanks North Star Alaska
    Denali Alaska
    Matanuska-Susitna Alaska
    Valdez-Cordova Alaska ...
share|improve this answer
    
I do not have arcgis or any other specialty geographic tools. Thank you for the zipfile and for providing detailed code. –  Sukotto Sep 19 '11 at 14:41
add comment

UPDATE: Added spatial index to improve performance and brief instructions for using this script on Windows.

#-------------------------------------------------------------------------------
#   This script will build an adjacency table in csv format representing
#   county polygons that "neighbor" each other. This script is intended
#   to illustrate the use of Python, OGR, Shapely, and Rtree.
#
#   County shapefile used in this example: http://cta.ornl.gov/transnet/scuov.zip
#
#   You will need Shapely, GDAL/OGR, and Rtree to run this. The easiest way for
#   Windows users to obtain these dependencies is to download the OSGeo4W
#   installer[1] and run an advanced install to retreive (located under Libs):
#
#       gdal-python (1.8.0)
#       libspatialindex (1.5.0-1, not 1.6.1)
#
#   From the OSGeo4W shell, install setuptools via ez_setup.py[2], and then run
#   these commands:
#
#       python -m easy_install shapely
#       python -m easy_install rtree
#
#   This script can then be run from the OSGeo4W shell.
#
#   [1] http://trac.osgeo.org/osgeo4w
#   [2] http://trac.osgeo.org/osgeo4w/wiki/TracPlugins
#-------------------------------------------------------------------------------
#!/usr/bin/env python

import sys, cPickle
from rtree import Rtree
from osgeo import ogr
from shapely import wkb

class FastRtree(Rtree):
    def dumps(self, obj):
        return cPickle.dumps(obj, -1)

def generator_function(list):
    for i, obj in enumerate(list):
        yield (i, obj[0].bounds, obj)

def main():
    csv = r'C:\county_adjacency.csv' # Path to output csv, modify as needed
    shapefile = r'C:\scuo.shp' # Path to existing shapefile, modify as needed
    fipsindex = 3 # Field index for FIPS code
    statenameindex = 0 # Field index for state name
    countynameindex = 1 # Field index for county name

    driver = ogr.GetDriverByName('ESRI Shapefile')
    dataset = driver.Open(shapefile)

    if dataset is None:
        print 'Open failed.'
        sys.exit(1)

    lyr = dataset.GetLayerByName('scuo') # Modify as needed
    lyr.ResetReading()

    countylist = []

    for feat in lyr:
        geom = feat.GetGeometryRef()
        if geom is not None and (geom.GetGeometryType() == ogr.wkbPolygon or geom.GetGeometryType() == ogr.wkbMultiPolygon):
            countytuple = wkb.loads(geom.ExportToWkb()), \
                          feat.GetFieldAsString(fipsindex), \
                          feat.GetFieldAsString(countynameindex), \
                          feat.GetFieldAsString(statenameindex)
            countylist.append(countytuple)
        else:
            print 'Error reading polygon geometry for: ' + \
                  feat.GetFieldAsString(fipsindex) + ', ' + \
                  feat.GetFieldAsString(countynameindex) + ', ' + \
                  feat.GetFieldAsString(statenameindex)

    idx = FastRtree(generator_function(countylist))

    csvfile = open(csv, 'w')
    csvfile.write('\"fips1\",\"county1\",\"state1\",\"fips2\",\"county2\",\"state2\"' + '\n')
    for i in countylist:
        for j in list(idx.intersection(i[0].bounds, objects='raw')):
            try:
                if i[0].touches(j[0]):
                    csvfile.write('\"' + \
                                  i[1] + '\",\"' + \
                                  i[2] + '\",\"' + \
                                  i[3] + '\",\"' + \
                                  j[1] + '\",\"' + \
                                  j[2] + '\",\"' + \
                                  j[3] + '\"\n')
            except:
                print 'Failed to evaluate: ' + \
                      i[1] + ',' + \
                      i[2] + ',' + \
                      i[3] + \
                      ' with ' + \
                      j[1] + ',' + \
                      j[2] + ',' + \
                      j[3]

    csvfile.close()
    dataset = None

if __name__ == '__main__':
    main()
share|improve this answer
add comment

The US Census Bureau publishes tons of spatial maps. One of these that is interesting is the "Counties and Equivalents".

These can be downloaded here:
http://www.census.gov/cgi-bin/geo/shapefiles2010/main

Once you have that, there is a shapefile importer for SQL Server here:
http://beyondrelational.com/blogs/jason/archive/2010/09/04/import-shapefiles-into-sql-server-and-aggregate-spatial-data-geometry.aspx

Once in SQL server, you can run an "intersects" query (which says that the boundaries of the counties should share at lease 1 point). Something like this:

SELECT c1.name, c2.name FROM
 county c1, county c2 WHERE c1.geog.STIntersects(c2.geog) AND c1.id < c2.id
share|improve this answer
    
But what if a c2.id is less than a c1.id and not the same county? –  Allan Adair Sep 15 '11 at 14:48
    
I think this is a good way to do it, but I would modify the query like so: SELECT c1.name, c2.name FROM county c1 INNER JOIN county c2 ON c1.geom.STTouches(c2.geom) = 1 (note the geometry type, not geography) –  Allan Adair Sep 15 '11 at 14:59
    
@Allan: The point of the c1.id < c2.id is so that every c1 gets matched up with each c2 exactly once. That is, every c1 is also a c2, this is a "self-join". –  John Gietzen Sep 15 '11 at 17:06
    
Ah, yes, I'm quite aware of self joins. I was questioning the use of c1.id < c2.id. I now see that what you are doing provides an answer to a question that I interpreted differently. –  Allan Adair Sep 15 '11 at 18:21
add comment

You could use Yahoo's geoplanet API

http://where.yahooapis.com/v1/place/12587851/neighbors?appid=[yourappidhere]

where 12587851 ist Polk, FL (County)

  • see the documentation at geoplanet API Reference

  • a representation of it with geoplanet explorer

  • all counties of a US state can be retrieved by

     http://where.yahooapis.com/v1/counties/CA?appid=[yourappidhere]
    

    where CA is replaced by the state code (e.g. California)

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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