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

GStav
  • 309
  • 3
  • 12
  • 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 '22 at 14:57

2 Answers2

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?

Pepe N O
  • 676
  • 4
  • 10
  • 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 '22 at 09: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

GStav
  • 309
  • 3
  • 12