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

At this point, I'd prefer NOT to use PostGIS because I might end up giving the tool away to ArcGIS users. I think 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.

GDAL

Here is a page where I talk about Building GDAL on Linux.

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!

1. Read from a shapefile, because that is what I have handy.

Explore data from command line using ogrinfo.

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

Code up some python.

2. Find the crossings. I am not worried about whether they are really road intersections yet.

3. Write to a feature class in file geodatabase

# Does the geodatabase exist?