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 want to calculate the angle between two lines where they intersect in PostGIS.

The starting point for angle calculations in PostGIS seems to be ST_Azimuth - but that takes points as input. My first thought was to take the endpoints of the intersecting lines and performing an Azimuth calculation on those. That is not good enough, because most of the line features are not straight, and I am interested in the angle at intersection. So what I came up with is a nested operation that goes through the following steps:

  1. Identify all the intersections between the two line feature tables.
  2. Create a very small buffer around the intersection point
  3. Identify the points where the line features intersect the buffer exterior (taking the first point if there are more than one - I'm really only interested in whether the angle is close to 0, 90 or 180 degrees)
  4. Calculate ST_Azimuth for those two points.

The full SQL is kind of long to post here, but I gisted it here if you're interested. (By the way, is there a better way than to carry over all the fields going down the WITH statements?)

The results don't look right, so I'm clearly doing something wrong:

output example 1 output example 2

EDIT I redid the calculations in EPSG:3785 and the results are a little different but still not right:

output in 3785 #1 output in 3785 #2

My question is where the flaws are in this process. Am I misunderstanding what ST_Azimuth does? Is there a CRS issue? Something else altogether? Or maybe there's a much, much simpler way to do this?

share|improve this question
1  
What was the original CRS? Angle calculations should be made with a conformal projection—not with unprojected lat/long (SRID=4326). –  Mike T May 9 '12 at 20:30
    
It was EPSG:4326 coordinates originally, I included the ST_Translate just to be 100% sure all processing would be done in the same CRS. I will try a conformal projection, thanks. –  mvexel May 9 '12 at 20:39
    
I redid the calculations is EPSG:3785 and it does make a difference - I'll amend the question to show the new results - but the result still does not reflect the actual angle. –  mvexel May 9 '12 at 21:11

2 Answers 2

up vote 6 down vote accepted

I had the epiphany. It is rather mundane. I was leaving out one essential piece of information for PostGIS to calculate the right angle.

What I was calculating was the angle between only the two points intersecting the small buffer exterior. To calculate the angle of the intersection, I need to calculate both angles between both points on the buffer exterior and the intersection point of the two line features and subtract them.

I updated the full SQL, but here's the salient bit:

SELECT
    ...
    abs
    (
        round
        (
            degrees
            (
            ST_Azimuth
            (
                points.point2,
                points.intersection
            )
            -
            ST_Azimuth
            (
                points.point1,
                points.intersection
            )
        )::decimal % 180.0
        ,2
    )
)
AS angle
...
FROM
points 
share|improve this answer
1  
I was thinking about the angle of the buffered point w.r.t. the instersection, but I haven't any time to go in detail. Another aspect is the angular units. You need to multiply the result in radians from ST_Azimuth by 180.0/pi() to get results in degrees. –  Mike T May 9 '12 at 23:52
    
Yep thanks, I use the PostgreSQL degrees() function for that. –  mvexel May 9 '12 at 23:55
    
Cleaver. (I didn't even know there was a degrees function until now.) It would be nice to wrap all this logic up in a function call, but I'm having difficulty in conceptualizing how it would work, i.e. ST_IntersectionAngle(...? –  Mike T May 10 '12 at 0:24
    
I was actually surprised that it is not a PostGIS function. Thanks for your feedback on this. –  mvexel May 10 '12 at 1:48

I recently had to calculate the same thing, but decided on a simpler and likely faster approach.

To find the extra points for the azimuth calculation, I just check a permyriad of the length behind the intersection (or after in the rare case that it happens on the very start of the line) using ST_Line_Locate_Point and ST_Line_Interpolate_Point:

abs(degrees( 
  ST_Azimuth (
    intersection, 
    ST_Line_Interpolate_Point(
      line1, 
      abs(ST_Line_Locate_Point(line1, intersection) - 0.0001)
    )
  )
  -
  ST_Azimuth (
    intersection, 
    ST_Line_Interpolate_Point(
      line2, 
      abs(ST_Line_Locate_Point(line2, intersection) - 0.0001)
    )
  )
))

The permyriad was arbitrary and for more consistent results it would be better to use an absolute offset. To for example check 20m beforehand, you'd change 0.0001 to 20/ST_Length(line1) and 20/ST_Length(line2) respectively.

share|improve this answer

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.