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

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

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