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

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.

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)... Commented Sep 5, 2014 at 18:42
• Please edit this question to clarify that the coordinates are latitudes and longitudes: that completely changes the correct answer! Commented Sep 5, 2014 at 19:44
• @Vince: That's a strange argument. Any technical discipline can be broken down to "just basic math" (or statistics). Commented May 15, 2015 at 18:54
• Commented May 15, 2015 at 19:32
• @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? Commented May 15, 2015 at 22:56

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:

fAC = sin βAC

gAC = cos βAC

fBC = sin βBC

gBC = cos βBC

then calculate a determinant

d = fBC gAC - fAC gBC

(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 = yB - yA

Dx = xB - xA

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

sAC = (fBC Dy - gBC Dx) / d

sBC = (fAC Dy - gAC Dx) / d

where sAC is the distance from A to the intersection point, C, along its line-of-sight,
and sBC 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

xC = xA + fAC sAC

yC = yA + gAC sAC

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

xC = xB + fBC sBC

yC = yB + gBC sBC

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

``````

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

1. Project the situation onto the plane, solve using the simple planar coordinates method, and reverse project back to the spheroid.

2. 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.

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.

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.

• You also need to trap for an undefined slope (y1 == y2) Commented Sep 5, 2014 at 18:29
• in my case xa and ya refer to the point coordinates, not slope. They can be zero. Commented Sep 5, 2014 at 18:30
• In your first pair of equations, ya==yc or yb==yc Commented Sep 5, 2014 at 18:31
• 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. Commented Sep 5, 2014 at 19:28
• 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. Commented Sep 5, 2014 at 19:46