How to calculate the angle at which two lines intersect in PostGIS?

Question

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:

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

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?

Answer 1

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 

Answer 2

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.

Answer 3

@mvexel's answer is very helpful. However, that method does not preserve the line's direction. My approach is

1. Identify all the intersections between the two line feature tables.
2. Create a very small buffer around the intersection point.
3. Using the buffer to intersect with the otherway. Geting the endpoint of intersected line as point 1.
4. Using the buffer to intersect with the motorway. Geting the endpoint of intersected line as point 2.
5. Calculate ST_Azimuth for those two points

code link

First: get intersection points.

with intersections as (
select ST_Intersection(a.geometry, b.geometry) as intersection,
       a.geometry as otherway_geom,
       b.geometry as motorway_geom,
       a.id as otherway_osmid,
       b.id as motorway_osmid

from otherways a, motorways b
where ST_Intersects(a.geometry,b.geometry)),

Second, create a buffer around intersection point.

buffers AS (SELECT intersections.intersection, 
    ST_ExteriorRing (ST_Buffer(intersections.intersection, 1))AS extring,
    ST_Buffer(intersections.intersection, 1) as buffer,
    intersections.otherway_geom,
    intersections.motorway_geom,
    intersections.otherway_osmid,
    intersections.motorway_osmid
FROM 
    intersections),

Third, get the end point of each line, which preserve the direction of linestring.

points AS(

SELECT 
    ST_Intersection(buffers.buffer, buffers.otherway_geom) as otherway_intersected_lane, 
    ST_Intersection(buffers.buffer, buffers.motorway_geom) as motorway_intersected_lane,
    st_endpoint(ST_Intersection(buffers.buffer, buffers.otherway_geom)) as point1,
    st_endpoint(ST_Intersection(buffers.buffer, buffers.motorway_geom)) as point2,
    buffers.intersection,
    buffers.extring,
    buffers.buffer,
    buffers.otherway_geom,
    buffers.motorway_geom,
    buffers.otherway_osmid,
    buffers.motorway_osmid
FROM 
    buffers)

Finally, calculate the angle!

SELECT 

st_astext(points.point1) as point1,
st_astext(points.point2) as point2,
st_astext(points.extring) as extring,
st_astext(points.otherway_geom) as otherway_geom,
st_astext(points.motorway_geom) as motorway_geom,
abs
(
    round
    (
        degrees
        (
            ST_Azimuth
            (
                points.point2,
                points.intersection
            )

            -

            ST_Azimuth
            (
                points.point1,
                points.intersection
            )           
        )::decimal % 180.0
        ,2
    )
)AS angle_Azimuth,   
points.otherway_osmid AS otherway_id, 
points.motorway_osmid AS motorway_id 
FROM points