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.
- Read centerline data with OGR (part of GDAL).
- Find all the crossings with Shapely.
- 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?
