2

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.

  • 1
    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 5 '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 5 '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
  • 1
    @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
-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.

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

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

0

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.

0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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