# 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?

• What was the original CRS? Angle calculations should be made with a conformal projection—not with unprojected lat/long (SRID=4326). May 9, 2012 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. May 9, 2012 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. May 9, 2012 at 21:11

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
``````
• 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. May 9, 2012 at 23:52
• Yep thanks, I use the PostgreSQL degrees() function for that. May 9, 2012 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(...`? May 10, 2012 at 0:24
• I was actually surprised that it is not a PostGIS function. Thanks for your feedback on this. May 10, 2012 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.

@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

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
``````