How to get intersection point from two lines (two coordinated points and two bearings)

Question

I have tried for ages to find the algorithm for getting the coordinate of two intersecting lines from two given coordinated (latitude and longitude) points.

As the picture below suggests, I have two points (A and B) with known coordinates. With these two coordinates, come two bearings, or directions. I am looking for a way to find coordinate of C, where the lines-of-sight intersect.

proof

Note: I am going to build it into a Java / C# Application. If there exists a library that already contains the math, that would help too.

As stated, this isn't a GIS question, just basic math (convert lines into Ax + By = C form, then solve with matrix algebra). Most GIS packages have a way to locate the intersection of a pair of lines, and they also will handle non-intersection, and intersection along a line. Then you get into non-Cartesian solutions, when the coordinates are angular (either spherical or spheroidal)... — Vince, Sep 05 '14 at 18:42
Please edit this question to clarify that the coordinates are latitudes and longitudes: that completely changes the correct answer! — whuber, Sep 05 '14 at 19:44
@Vince: That's a strange argument. Any technical discipline can be broken down to "just basic math" (or statistics). — Martin F, May 15 '15 at 18:54
@whuber: Given that OP accepted a planar coord answer, is it fair to edit the Q to reflect that, and let this be the canonical Q for line-intersection on the plane? — Martin F, May 15 '15 at 22:56
@Martin That's one approach. Another is to offer a better answer that includes planar and spherical (or even ellipsoidal) solutions. — whuber, May 16 '15 at 19:38
https://math.stackexchange.com/questions/1934683/how-to-calculate-an-intersection-point-between-two-coordinates-and-their-bearing — neogeomat, Aug 02 '19 at 07:59

Martin F · Answer 1 · 2015-05-16T23:16:43.927

Assuming that we're dealing with the planar coordinate case (that is not actually what the OP suggested, but I offer this as a better answer to the one given so far – and so far, accepted, by the OP – for the planar case), it helps to first determine the direction cosines from the two clockwise bearings, β_AC and β_BC, from known points A and B to unknown intersection point C:

f_AC = sin β_AC

g_AC = cos β_AC

f_BC = sin β_BC

g_BC = cos β_BC

then calculate a determinant

d = f_BC g_AC - f_AC g_BC

(If d is very close to zero, the two directions are near parallel and, essentially, there is no intersection.)

Next, calculate a couple of coordinate differences between known points A and B,

Dy = y_B - y_A

Dx = x_B - x_A

(where x is Easting and y is Northing) and use them to obtain two very useful distances

s_AC = (f_BC Dy - g_BC Dx) / d

s_BC = (f_AC Dy - g_AC Dx) / d

where s_AC is the distance from A to the intersection point, C, along its line-of-sight,
and s_BC is the distance from B to C, along its line-of-sight.

If either one of those distances is negative, the intersection point is behind the line-of-sight. That is, there is no real intersection in the direction of the bearing.

Finally, calculate the intersection coordinates via those of known point A and its distance and direction cosines

x_C = x_A + f_AC s_AC

y_C = y_A + g_AC s_AC

And while you're at it, do a check from the other point

x_C = x_B + f_BC s_BC

y_C = y_B + g_BC s_BC

score 1 · Answer 2 · edited Aug 02 '19 at 03:11

Or if you are happy with just plane trig solution in just 19 lines of code. Here it is in U-Basic.

This differs from the algorithms for geodesics solution by about 6mm at a distance of 175 meters.

    5   'interx.ub adapted from code at
    6   'https://rosettacode.org/wiki/Find_the_intersection_of_two_linesSinclair_ZX81_BASIC
    7   'In U Basic by yuji kida https://en.wikipedia.org/wiki/UBASIC
   10   XA=48.7815144526:'669595.708
   20   YA=-117.2847245001:'2495736.332
   30   XB=48.7815093807:'669533.412
   40   YB=-117.2901673467:'2494425.458
   50   XC=48.7824947147:'669595.708
   60   YC=-117.28751374:'2495736.332
   70   XD=48.77996737:'669331.214
   80   YD=-117.2922957:'2494260.804
   90   print "THE TWO LINES ARE:"
  100   print "YAB=";YA-XA*((YB-YA)/(XB-XA));"+X*";((YB-YA)/(XB-XA))
  110   print "YCD=";YC-XC*((YD-YC)/(XD-XC));"+X*";((YD-YC)/(XD-XC))
  120   X=((YC-XC*((YD-YC)/(XD-XC)))-(YA-XA*((YB-YA)/(XB-XA))))/(((YB-YA)/(XB-XA))-((YD-YC)/(XD-XC)))
  130   print "Lat =  ";X
  140   Y=YA-XA*((YB-YA)/(XB-XA))+X*((YB-YA)/(XB-XA))
  150   print "Lon = ";Y
  160   'print "YCD=";YC-XC*((YD-YC)/(XD-XC))+X*((YD-YC)/(XD-XC))

score 0 · Answer 3 · edited Apr 13 '17 at 12:33

Assuming that you are dealing with geographic coordinates and geographic (true north) bearings, you have a couple of options:

Project the situation onto the plane, solve using the simple planar coordinates method, and reverse project back to the spheroid.
Perform geodetic or spherical geometry calculations.

If your distances are not too great and your accuracy requirements are not too tight, the first is probably the easier solution. Broken down, this amounts to

a. Project the two geographic points to plane grid points via a suitable map projection. Ideally, this would be Gnomonic projection because it projects great circles as straight lines. (For a treatise on that, see Algorithms for geodesics by CFF Karney, aka @cffk.) However, a UTM projection method is probably the more common and suitable alternative available. Yet another compromise, suggested by @whuber in a comment to calculating-the-distance-between-a-point-and-a-virtual-line-of-two-lat-lngs, is to treat the geographic coordinates as though they were planar, pre-multiply the difference in longitudes (between A and B) by the cosine of their average latitudes.

b. Convert the geographic azimuths to grid bearings via the meridian convergence. In theory, I believe all well-defined projections have computable grid convergences, as well as scale factors. (See calculate grid convergence: True North to Grid North for an example for the UTM case.)

c. Calculate the point of intersection via the planar method. (See my other answer here for this same question.)

d. Reverse project the intersection point back from grid to geographic coordinates. (Projections, in either direction, are the subject of too many questions on GIS SE to mention here.)

I believe you can achieve this effect in PostGIS via the ST_Project() and ST_Intersection() functions:

geogC ← ST_Intersection (ST_Project (geogA, bigdist, azimA), ST_Project (geogB, bigdist, azimB))

where

geogA and azimA are geographic point A and azimuth at A, respectively,
geogB and azimB are geographic point B and azimuth at B, respectively, and
bigdist is an arbitrary but appropriately large enough distance to ensure an intersection.

I conclude that such use of the ST_Project() and ST_Intersection() is equivalent to the algorithm I describe because of the note in the user manual for ST_Intersection():

For geography this is really a thin wrapper around the geometry implementation. It first determines the best SRID that fits the bounding box of the 2 geography objects (if geography objects are within one half zone UTM but not same UTM will pick one of those) (favoring UTM or Lambert Azimuthal Equal Area (LAEA) north/south pole, and falling back on mercator in worst case scenario) and then intersection in that best fit planar spatial ref and retransforms back to WGS84 geography.

As for the second, more rigorous solution, I do not know of even a spherical trigonometry solution, let alone an ellipsoidal trigonometry one. I suspect the solution is very difficult to explain here in simple terms.

score 0 · Answer 4 · answered May 19 '15 at 21:45

The solution in terms of geodesics is given in Section 8 of Algorithms for geodesics. This uses the ellipsoidal gnomonic projection to convert the problem to an equivalent planar problem. An implementation of the solution in C++ is given here. GeographicLib includes the necessary underlying geodesic routines in both C# and Java.

score -1 · Accepted Answer · answered Sep 05 '14 at 18:25

-1

This is basicmath: With

(xa -xc)/(ya-yc) = d_a
(xb -xc)/(yb-yc) = d_b

or

xa -xc =d_a(ya-yc)
xb - xc = d_b(yb-yc)

substracting both from eachother:

xa -xb = d_a(ya-yc) -  d_b(yb-yc)
d_a*yc -d_b*yc = d_a*ya - d_b*yb -xa +xb

so basically

yc  = (d_a*ya - d_b*yb -xa +xb)/(d_a -d_b)

and

xc = xa-d_a(ya-yc)

Do a check that d_a<>d_b, because you will not have an intersection if your bearing is equal.

answered Sep 05 '14 at 18:25

johanvdw

6,207
27
42

1

You also need to trap for an undefined slope (y1 == y2) – Vince Sep 05 '14 at 18:29
in my case xa and ya refer to the point coordinates, not slope. They can be zero. – johanvdw Sep 05 '14 at 18:30
In your first pair of equations, ya==yc or yb==yc – Vince Sep 05 '14 at 18:31
2

I see that it is basic algebra. I was just sure, that the sphere of earth had to be taken into account, when doing calculations on lat long, hence the reason I wrote it here. Thank you. – Wildcard Sep 05 '14 at 19:28
4

The question needs clarification then, because this answer is correct only when the coordinates are projected and the bearing has been computed in that projection. It does not apply to (lat, lon) ("geographic") coordinates with geodetic bearings. – whuber Sep 05 '14 at 19:46
Unstated in the answer is that d_a and d_b are tangents of the bearings from A and B. And those will break down for bearings close to +/- 90. Also, it doesn't state how to detect diverging bearings. – Martin F May 15 '15 at 23:03

How to get intersection point from two lines (two coordinated points and two bearings)

5 Answers5