Finding Street Intersections in PostGIS

I recently had a colleague approach me needing a street intersection search added to our enterprise REST web services. I hadn’t really thought of it before, but it didn’t seem too difficult. I figured out some SQL and rammed it in our web service framework. 30 minutes later, a new REST service entered the world.

Here was the basic SQL (streets to look for: RUTH and LAKEDELL; centerline table is ROADS):

select distinct(AsText(intersection(b.the_geom, a.the_geom))) as the_intersection from (select the_geom from roads where streetname = 'RUTH') a, (select the_geom from roads where streetname = 'LAKEDELL') b where a.the_geom && b.the_geom and intersects (b.the_geom, a.the_geom)

Basically I'm getting the line geometry from the roads layer in a way that it acts like two different tables - one for RUTH and one for LAKEDELL. Then it finds the line pairs that intersect and returns the intersection(s) like so:

POINT(1472605.08 549803.579999998)

This is good enough for government work, but I starting thinking about what it was doing. What if you had two lines that cross but didn't intersect, like a road going under a bridge. This call would flag it as an intersection, which might not be what the user wants. So then I tried this:

select distinct(AsText(intersection(b.the_geom, a.the_geom))) as the_intersection from (select the_geom from roads where streetname = 'RUTH') a, (select the_geom from roads where streetname = 'LAKEDELL') b where a.the_geom && b.the_geom and touches(b.the_geom, a.the_geom) and not(ST_Crosses(a.the_geom,b.the_geom))

This query finds the line segments that touch but don't cross, eliminating those segments that cross but don't intersect (i.e. no node). This is better, but still not perfect. The problem here would be a line that intersects a road segment at one point but crosses without intersecting at another point. This query would toss those right out of contention.

Of course, that would very rarely happen. In fact, the first query is still probably good enough for government work. But suppose we’re really anal about this one. Here’s what I came up with:

select AsText(a.the_geom) as the_intersection from (select ST_StartPoint(the_geom) as the_geom from roads where streetname = 'RUTH' union select ST_EndPoint(the_geom) as the_geom from roads where streetname = 'RUTH') a, (select ST_StartPoint(the_geom) as the_geom from roads where streetname = 'LAKEDELL' union select ST_EndPoint(the_geom) as the_geom from roads where streetname = 'LAKEDELL' ) b where a.the_geom && b.the_geom and a.the_geom = b.the_geom

I decided the problem with the first two queries was I was using the wrong geometry type. I really wanted to go point hunting, not line hunting. So I again create two SQL selects for RUTH and LAKEDELL, but this time I create a point geometry from the start and (via union) end points of the line segments. So now instead of two sets of lines, I have two sets of points. Perfect. Now I find the points from one that are identical with the other, and I have real node-intersections with no false positives or negatives.

This query could probably get tuned a bit (right now on 40k record roads layer it runs in ~90ms), but you get the idea. I thought I’d post it here because I thought the thought process was interesting and to show some of the flexibility you have when geoprocessing with PostGIS. For the record, I’m sticking with the first query for our production web service, as with it people want to get to an area and they don’t really care if one line of pavement hits another - they just want the map to pan to RUTH and LAKEDELL. But if your needs are different, it’s a good idea to think about the problem thoroughly and apply your tools with an eye toward what they’re really doing.