2

I am trying to find the intersection point of two lines of which I know 2 points.

For example, I have these 2 line segments:

select st_geometryfromtext('LINESTRING (4.50554487780521 51.922119943633575, 4.504656820078167 51.92231795217855)', 4326),
st_geometryfromtext('LINESTRING (4.505476754241158 51.92221504789901, 4.505379267847784 51.92221833721103)', 4326)

As it can be seen on the map,they do not have an intersection point:

enter image description here

What I am interested in finding is the intersection point of the lines which these segments belong to, as shown below:

enter image description here

1
  • 1
    The 2 line segments are the ones in the photo, and the ones you would see if you run the query. I want to find the intersection of the lines passing from these segments (because the segments do not intersect; so st_intersection does not return anything)
    – GStav
    Feb 24, 2022 at 14:57

2 Answers 2

4

Just extend (extrapolate) one of your lines and get the intersection point as:

select st_intersection((select st_geometryfromtext('LINESTRING (4.50554487780521 51.922119943633575, 4.504656820078167 51.92231795217855)', 4326)),
(select ST_MakeLine(ST_TRANSLATE(a, sin(az1) * len, cos(az1) * 
len),ST_TRANSLATE(b,sin(az2) * len, cos(az2) * len))
  FROM (
    SELECT a, b, ST_Azimuth(a,b) AS az1, ST_Azimuth(b, a) AS az2, ST_Distance(a,b) + 1 AS len
      FROM (
        SELECT ST_StartPoint(st_geometryfromtext('LINESTRING (4.505476754241158 51.92221504789901, 4.505379267847784 51.92221833721103)', 4326)) AS a, 
        ST_EndPoint(st_geometryfromtext('LINESTRING (4.505476754241158 51.92221504789901, 4.505379267847784 51.92221833721103)', 4326)) AS b
          FROM ST_MakeLine(ST_MakePoint(1,2), ST_MakePoint(3,4)) AS the_geom
    ) AS sub
) AS sub2));

extrapolation example from How to extend a straight line in postgis?

1
  • I have already accepted your answer but it is good to notice that it assumes we know which line has to be extended. I did manage to solve my problem though - so thanks a looot!
    – GStav
    Feb 25, 2022 at 9:41
1

If someone is looking for a version where you don't know which line has to be extended, here is a variation of Pepe N O's answer:

with segments as ( 
    SELECT ST_StartPoint(st_geometryfromtext('LINESTRING (4.505476754241158 51.92221504789901, 4.505379267847784 51.92221833721103)', 4326)) AS s1_a, 
           ST_EndPoint(st_geometryfromtext('LINESTRING (4.505476754241158 51.92221504789901, 4.505379267847784 51.92221833721103)', 4326)) AS s1_b,
           ST_StartPoint(st_geometryfromtext('LINESTRING (4.50554487780521 51.922119943633575, 4.504656820078167 51.92231795217855)', 4326)) AS s2_a, 
           ST_EndPoint(st_geometryfromtext('LINESTRING (4.50554487780521 51.922119943633575, 4.504656820078167 51.92231795217855)', 4326)) AS s2_b
)
, azimuths as (
select 
    *,
    ST_Azimuth(s1_a, s1_b) AS s1_az1,
    ST_Azimuth(s1_b, s1_a) AS s1_az2,
    1 AS s1_len,
    ST_Azimuth(s2_a, s2_b) AS s2_az1,
    ST_Azimuth(s2_b, s2_a) AS s2_az2,
    1 AS s2_len 
from segments
)
select
st_intersection(
ST_MakeLine(ST_TRANSLATE(s1_b::geometry, sin(s1_az1) * s1_len, cos(s1_az1) * s1_len), ST_TRANSLATE(s1_a::geometry,sin(s1_az2) * s1_len, cos(s1_az2) * s1_len)),
ST_MakeLine(ST_TRANSLATE(s2_b::geometry, sin(s2_az1) * s2_len, cos(s2_az1) * s2_len), ST_TRANSLATE(s2_a::geometry,sin(s2_az2) * s2_len, cos(s2_az2) * s2_len))
)
from azimuths;

This is the intersection of the extended lines the above query would give you:

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.