Finding Intersections with Python

Where roads meet

We need to be able to build a table of all the places roads meet.

• This is not a page about the general case of geometric intersections. Road intersections.
• And we don't want to find the places that roads CROSS (over- or under-crossings). Just where they join.

Both of these points are salient because they mean dealing with dirty input data. Always the case.

• Roads should always meet at a vertex
• Roads that meet should always have a common vertex, not just be "close".
• Roads that do not meet (crossings) should NOT meet these conditions.
  • Any gaps should indicate no connection.
  • Any place roads touch but not at a vertex should indicate a crossing.

But none of this is ever true! Alas, we always get data that looks good on a map but is not routable.

To this end, I want a tool that:

• Generates a point layer of intersections
• Generates a point layer of error cases

The tools

I have a version of this program already written in C# using ArcObjects, and it's a BIG MESS. Worst of all, every time I do a build in a new release of ArcGIS, it breaks!

I want something easy to maintain, easy to update, easy to run, and reasonably fast.

To accomplish this I am reading Python Geospatial Development by Erik Westra. Packt Pubs

I've been working a lot with PostGIS lately and it would be ideal for this project. HOWEVER I need to be able to give it to some ESRI users and I don't want to teach them how to use PostGIS as part of that project. This means I will be using Shapely.

I want to be able to use File Geodatabases, so I will be using a custom built GDAL with Python bindings on Linux and the OSGeo4W version on Windows.

GDAL

Here is a page where I talk about Building GDAL on Linux. I guess I should also do a write up on Using OSGeo4W but I have not done that yet.

Shapely

Here is the home of Shapely, where it is described as "PostGIS-ish operations outside a database context for Pythoneers and Pythonistas." Sounds like what I need.

The steps

Starting off with the simplest version, ignoring all the data cleaning.

1. Read centerline data with OGR (part of GDAL).
2. Find all the crossings with Shapely.
3. Write points with appropriate attributes (FENAME1, FENAME2)

With the tools at our disposal, this should be so easy! It turns out that it is!.

I was able to write a script in Python and develop it on a Linux system, and the first time I ran it on Windows, it ran successfully to completion with NO CHANGES. Awesome

The next part now will be getting it to run in ESRI Model Builder.

Explore data from command line using ogrinfo.

ogrinfo -so cntrline.shp cntrline
Gives lots of info I won't include here

Read shapes into Shapely geometries

Time to code up some python. I want to be able to read ESRI File Geodatabases too. And write them.

class InputFeatures:
    
    # Input fields
    input_fieldname = "FULLNAME"

    # Stubby little streets less than this are ignored.
    MinLineLength = 3.0
    
    FeatureList = [] # tuple containing (hashed name, ogr feature, shapely geom)
    SpatialReference = None
    
    BadNameCount = 0
    EmptyNameCount = 0
    BadGeometryCount = 0
    
    DebugFeatureCount = 0 # Set to some value to load a few features only -- good for testing
    
    _d_ignored_names = {} # ALLEY, UNKNOWN...

    def __init__(self, input_fieldname, ignored_path, st_types_path):
        self.input_fieldname = input_fieldname
        
        icount = self._load_ignored_names(ignored_path)
        if icount == 0:
            print "Could not load ignored names file '%s'" % ignored_path

        self.condenser = streetnames.Condense()
        ccount = self.condenser.LoadStreetTypes(st_types_path)    
        if ccount == 0:
            print "Could not load street types file '%s'" % st_types_path
            
        return

    #####
    def LoadFeatures(self, fcname):
        """ Load everything from the input feature layer into memory.
            Well.. not EVERYTHING. Ignores meaningless features. """
        feature_count = 0
        (ogr_ds, ogr_layer) = self._open_layer(fcname)
        if not ogr_layer: 
            print "Can't open input layer '%s'" % fcname
        else:
            feature_count = self._loadFeatures(ogr_layer)  
            self.SpatialReference = ogr_layer.GetSpatialRef()
            ogr_ds.Destroy()
        return feature_count

    def _open_layer(self, fcname):
        (workspacename, layername) = fc_to_ogrname(fcname)
        ds = None
        layer = None
        try:
            ogr.RegisterAll()
            ds = ogr.Open(workspacename);
            layer = None
            if ds:
                layer = ds.GetLayerByName(layername)
        except Exception:
            print "Could not open input feature class '%s" % fcname
        return (ds, layer)

    def _loadFeatures(self, layer):
        max_feature_count = layer.GetFeatureCount()
        debug_count = self.DebugFeatureCount
        if (debug_count > 0) and (debug_count < max_feature_count):
            max_feature_count = debug_count

        for i in range(max_feature_count):
            if progress and not i % 1000: print i, "..",
                
            feature = layer.GetFeature(i)
            if not feature:
                self.BadGeometryCount += 1
                continue
       
            name = feature.GetField(self.input_fieldname)
            if not name:
                self.EmptyNameCount += 1
                continue

            # Ignore all streets with names like ALLEY and UNKNOWN
            if self._d_ignored_names.has_key(name):
                self._d_ignored_names[name] += 1 # count each one to help debugging
                self.BadNameCount += 1;
                continue

            try:
                # Convert OGR geometry to Shapely geometry
                geometry = feature.GetGeometryRef()
                shgeom = shapely.wkt.loads(geometry.ExportToWkt())
            except:
                print "%d : %s (EXCEPTION)" % (i, name)
                self.BadGeometryCount += 1
                continue
    
            msg = self._featureIsValid(shgeom)
            if msg :
                print "%d : %s : %s" % (i, name, msg)
                self.BadGeometryCount += 1
                continue
    
            # Geometry is good and name is good, so add it to our list
       
           # Make a condensed name from the FULLNAME field.

            # I wonder what would happen if I actually used hash() here
            hashed = self.condenser.HashStreetName(name)
            self.FeatureList.append( (hashed, feature, shgeom) )

        return len(self.FeatureList)
    
    def _featureIsValid(self, shgeom):
        """ Return True if this feature passes validity tests. """
        
        if not shgeom.is_valid:
            return "Geometry is not valid."
        
        if shgeom.is_empty:
            return "Geometry is empty."
        
        mp = list(shgeom.boundary.geoms)
        if len(mp) == 0:
            return "Geometry is empty!"
        
        if shgeom.length < self.MinLineLength:
            return "Short line ignored."
      
        return None

def isshapefilename(fcname):
    if fcname[-4:].lower() == '.shp':
        return True
    return False

def fc_to_ogrname(fcname):
    """ Convert an ESRI feature class name into an OGR workspace and layername """
    
    try:
        if isshapefilename(fcname):
            # We're working with a shapefile.
            # Don't test for existence because we might be making a new shapefile.
            workspacename = fcname
            layername = os.path.basename(fcname)[:-4]
            return (workspacename, layername)
    except Exception:
        pass
    
    try:
        directory = os.path.dirname(fcname)
        if directory[-4:].lower() == '.gdb':
            # We're working with a file geodatabase
            if not os.path.exists(directory):
                print "File geodatabase '%s' does not exist" % directory
                return (None, None)
            workspacename = directory
            layername = os.path.basename(fcname)
            return (workspacename, layername)
    except Exception:
        pass

    print "Something is wrong with feature class '%s'" % fcname
    return (None,None)

2. Use Shapely to find the intersections.

First I find where lines cross.

Then I look at the endpoints of each line segment.

class GenerateCrossPoints:
    
    # Units depend on projection of data source, usually in feet.
    BufferDistance = 3.0

    # Trim pre-directionals from output street names when set true
    trim_predirectionals = True 

    # Input 
    input_fieldname = "FULLNAME"
    featureList = None

    # Dictionaries for fast look up tables.
    d_street_types = {} # ST AV BLVD
    d_street_translations = {}  # stuff like first => 1

    # OUTPUT
    # Has to be a list because of dupe names + modifier (1st and MAIN 1, 1st and MAIN 2,...)
    IntersectionList = [] # Tuple (geom, (hash1, hash2), fename1, fename2, modifer, rating)
    ErrorList = []        # Tuple (geom, fename1, fename2, errdesc)
    # used to differentiate between intersections with same name
    _d_modifier = {}      # [hashed names] => counter count instances of each name

    # rating, lower is better
    # 0 = vertex at line crossing
    # 1 = endpoint touches a vertex
    # 2 = crossing
    # 3 = endpoint touches a line
    # 4 = line is near vertex
    # 5 = endpoint is near a line

    # Regular expressions, compiled ahead of need

    # Used for trimming predirectionals
    re_predir  = re.compile(r'^(N|S|W|W|NE|NW|SE|SW)\s(.*)$', re.IGNORECASE)

    #####################################################

    def __init__(self, input_fieldname):
        self.input_fieldname = input_fieldname
        pass   
    
    #########
    def FindIntersections(self, featureList):
        self.featureList = featureList
        self._findLineIntersections()      # where lines cross
        self._findEndpointIntersections()  # where endpoints meet at vertices
        return

    def _getVertices(self, geom, points):
        """ Given a shapely geometry (line or multi)
            Return a list of the points in that object """
    
        if geom.type == "MultiLineString":
            for line in geom.geoms:
                self._getVertices(line, points)
        else:
            for p in geom.coords:
                points.append(shapely.geometry.Point(p))
        return
    
    def _adjustName(self, name):
        """ Adjust the fename according to trim flag. """
        if self.trim_predirectionals:
            mo = self.re_predir.search(name)
            if mo:
                return mo.group(2)
        return name
    
    def _proximate(self, keytuple, pt):
        """ Return True if the named pt is close to an existing intersection with the same name. """

        for (ex_pt, ex_keytuple, ex_fename1, ex_fename2, modifier, rating) in self.IntersectionList:
            if keytuple == ex_keytuple:
                # Has same name, is it nearby?
                buffered = pt.buffer(self.BufferDistance)
                if ex_pt.within(buffered):
                    if verbose:
                        print "Nearby + same name = ignored '%s x %s'" % keytuple
                    return True
        return False
    
    def _add_intersection(self, pt, keytuple, fename1, fename2, rating):
        """ Given geometry for a point, condensed name pair, full name pair, and rating,
        check to see if this is already in our list of intersections
        if not, then add it. """
        
        # Don't add if it already exists
    
        # 1ST x MAIN
        if (keytuple in self._d_modifier):
            if self._proximate(keytuple, pt): return;
    
        swapped = (keytuple[1], keytuple[0])  # MAIN x 1ST
        if (swapped in self._d_modifier):
            if self._proximate(swapped, pt): return;
            #print "Swapping names to remove duplication"
            t = fename1
            fename1 = fename2
            fename2 = t
            keytuple = swapped
    
        # Is there already an intersection with this name?
        if self._d_modifier.has_key(keytuple):
            modifier = self._d_modifier[keytuple]
            self._d_modifier[keytuple] += 1
        else:
            self._d_modifier[keytuple] = 1
            modifier = None
    
        # This should be the only place we need to trim the 2 street names, just before storing them.
        self.IntersectionList.append((pt, keytuple, self._adjustName(fename1), self._adjustName(fename2), modifier, rating))
        
        return
    
    #####
    def _findLineIntersections(self):
   
        print "Looking for line intersections"
        i = test_count = hit_count = 0
        for (condensedname, feature, shgeom) in self.featureList:
            for (othercondensedname, otherfeature, othergeom) in self.featureList:
                if progress and not i % 1000:
                    print i,
                    if not i % 10000: print
                i += 1
    
                test_count += 1
                if condensedname == othercondensedname:
                    # no self intersections allowed
                    continue
    
                if shgeom.intersects(othergeom) :
                    pt = shgeom.intersection(othergeom)
                    hit_count += 1
    
                    # two lines cross
                    rating = 2
                    
                    # is there a vertex here?
                    # or is there a vertex NEAR here?
                    vertices = []
                    self._getVertices(othergeom, vertices)
                    buffered = shgeom.buffer(self.BufferDistance)
                    
                    for vertex in vertices:
                        if shgeom.equals(vertex):
                            rating = 0 # vertex at line crossing
                            break
                        if vertex.within(buffered):
                            rating = 4 # line is near vertex
                            break;
    
                    fename1 = feature.GetField(self.input_fieldname) 
                    fename2 = otherfeature.GetField(self.input_fieldname)
                    keytuple = (condensedname, othercondensedname)
                    self._add_intersection(pt, keytuple, fename1, fename2, rating)
    
                    if verbose: print "Hit! '%s x %s'" % keytuple
            
        if progress: print
        return hit_count
    
    #####
    def _findEndpointIntersections(self):
        """ This should be run after generating line intersections.
            It checks the endpoints of every line. If there is an intersection
            there already, fine. Otherwise see if there is a line nearby and
            if so, add it. And endpoint is near but not touching a line. """
    
        print "Looking at endpoints.."  
        i = 0
        for (condensedname, feature, geom) in self.featureList:
            i += 1
            #if progress: print "Checking endpoints for %s" % (condensedname)
            
            fename1 = feature.GetField(self.input_fieldname)
            endpoints = list(geom.boundary.geoms)
    
            self._process_linestring_endpoint(condensedname, fename1, endpoints[0])
            self._process_linestring_endpoint(condensedname, fename1, endpoints[1])
        
        return
    
    #####
    def _process_linestring_endpoint(self, condensedname, fename1, endpoint) :
    
        # If there is already an intersection here,
        for (i_point, keytuple, ex_fename1, ex_fename2, modifier, rating) in self.IntersectionList:
            if endpoint.equals(i_point):
                # We're good already. Not an error
    #            if verbose: print "Endpoint already matched '%s'" % hashkey
                return
        
        # else
        #   is there a line nearby?
        buffered = endpoint.buffer(self.BufferDistance)
        test_count = hit_count = 0
        for (othercondensedname, otherfeature, othergeom) in self.featureList:
            test_count += 1
            if condensedname == othercondensedname:
                # no self intersection allowed
                pass
    
            # Are we touching
            elif endpoint.intersects(othergeom) :
                wkt = shapely.wkt.dumps(endpoint)
                hit_count += 1
    
                fename2 = otherfeature.GetField(self.input_fieldname)
                rating = 3 # endpoint actually touches a line
                
                keytuple = (condensedname, othercondensedname)
                self._add_intersection(endpoint, keytuple, fename1, fename2, rating)            
                if verbose: print "Endpoint touches line '%s x %s'" % keytuple
            
            # Are we not quite touching
            elif buffered.intersects(othergeom) :
                wkt = shapely.wkt.dumps(endpoint)
                hit_count += 1
    
                fename2 = otherfeature.GetField(self.input_fieldname)
                rating = 5 # endpoint is near a line
    
                keytuple = (condensedname, othercondensedname)
                self._add_intersection(endpoint, keytuple, fename1, fename2, rating)            
                if verbose: print "Endpoint near line '%s x %s'" % keytuple
    
        return hit_count

# end of class definition

3. Write intersection points to a feature class

This code can write to either a shapefile or a geodatabase. It looks at the string you pass in and decides. I use the ESRI concept of a "feature class" to make it easier to run this script from ESRI Model Builder.

• Shapefiles end in .shp for example "foo/bar/intersections.shp"
• Geodatabases have a layername at the end and a GDB in them, for example "goo/gump/features.gdb/intersections".

#####
def createOutputLayer(fcname, sref):
    
    (ogr_ds, ogr_layer) = create_layer(fcname, sref)
        # Add the attribute fields

    fieldDefn = ogr.FieldDefn(fename1_fieldname, ogr.OFTString);
    fieldDefn.SetWidth(40)
    ogr_layer.CreateField(fieldDefn)
    
    fieldDefn = ogr.FieldDefn(fename2_fieldname, ogr.OFTString);
    fieldDefn.SetWidth(40)
    ogr_layer.CreateField(fieldDefn)
    
    fieldDefn = ogr.FieldDefn(modifier_fieldname, ogr.OFTString);
    fieldDefn.SetWidth(40)
    ogr_layer.CreateField(fieldDefn)
    
    fieldDefn = ogr.FieldDefn(rating_fieldname, ogr.OFTInteger);
    ogr_layer.CreateField(fieldDefn)
    
    return (ogr_ds, ogr_layer)

def create_layer(fcname, sref):
    """ Create a layer for output """
    (ds_name, layer_name) = fc_to_ogrname(fcname)
    
    if isshapefilename(fcname):
        driver = ogr.GetDriverByName("ESRI Shapefile")
        if os.path.exists(ds_name):
            driver.DeleteDataSource(ds_name)
        ogr_ds = driver.CreateDataSource(ds_name)
    else:
        driver = ogr.GetDriverByName("FileGDB")
        if os.path.exists(ds_name):
            ogr_ds = driver.Open(ds_name)
            rval = ogr_ds.DeleteLayer(layer_name)
        else:
            ogr_ds = driver.CreateDataSource(ds_name)
    if ogr_ds is None:
        print "Create layer '%s' failed" % ds_name
        return (None, None)
    ogr_layer = ogr_ds.CreateLayer(layer_name, sref, ogr.wkbPoint)
    fieldDefn = ogr.FieldDefn(ID_FIELDNAME, ogr.OFTInteger)
    ogr_layer.CreateField(fieldDefn)

    return (ogr_ds, ogr_layer)

def saveIntersections(ogr_layer, intersectionList):
    global fename1_fieldname, fename2_fieldname, modifier_fieldname
        
    featureDefn = ogr_layer.GetLayerDefn()

    not_a_point = 0
    fid = 0
    for (p, keytuple, fename1, fename2, modifier, rating) in intersectionList:

        if p.type == "Point":
            wkt = shapely.wkt.dumps(p)
            feature = ogr.Feature(featureDefn)
            feature.SetGeometry(ogr.CreateGeometryFromWkt(wkt))

            fid += 1
            feature.SetField("id", fid);
            feature.SetField(fename1_fieldname, fename1)
            feature.SetField(fename2_fieldname, fename2)
            feature.SetField(rating_fieldname, rating)
            feature.SetField(modifier_fieldname, modifier) 
            ogr_layer.CreateFeature(feature)
                
            # Now create a second feature with the streetnames reversed.
            fid += 1
            feature.SetField("id", fid);
            feature.SetField(fename1_fieldname, fename2)
            feature.SetField(fename2_fieldname, fename1)
            feature.SetField(modifier_fieldname, modifier)
            ogr_layer.CreateFeature(feature)
            
            feature.Destroy()
        else:
            not_a_point += 1
    
    return fid