Intersection between two geodesic lines

[Gonzalo Pesquero]

Would it be possible to calculate the intersection point between two geodesic
lines using GeographicLib?

I've been reading the documentation but I haven't been able to find anything
regarding this functionality.

Thanks in advance,
Gonzalo

[Charles Karney]

The code below is a simple solution to the intersection problem
using the gnomonic projection. In this projection geodesics are
nearly straight; and they are exactly straight if they go
through the center of projection. So the intersection of two
geodesics can be found as follows: Guess an intersection point.
Project the two line segments into gnomonic, compute their
intersection in this projection, use this intersection point as
the new center, and repeat.

The attached code implements this (without any attempt to figure
out a termination condition). Here's an example of its
operation. Note that the 4 points must be in the same
hemisphere centered at the intersection point for the gnomonic
projection to be defined.

Line A: Istanbul - Washington
Line B: Reyjkavik - Accra

g++ -I/usr/local -L/usr/local -Wl,"-rpath=/usr/local/lib" \
-lGeographic -o intersect intersect.cpp

echo 42 29 39 -77 64 -22 6 0 | ./intersect

-->

Initial guess 37.75 -17.5
Increment 16.97176336894064 2.933476949591608
Increment -0.004733759685677796 0.002667305917045226
Increment -3.072813115068129e-10 1.835989138498917e-10
Increment 1.4210854715202e-14 -5.329070518200751e-15
Increment 0 0
Increment 0 0
Increment 0 0
Increment 0 0
Increment 0 0
Increment 0 0
Final result 54.7170296089477 -14.56385574430775
Azimuths on line A -84.14506440848534 -84.14506440848535
Azimuths on line B 160.9174944697586 160.9174944697586

intersect.cpp is attached

[Gonzalo Pesquero]

Hello,

Thank you very much for your fast response. I’ve compiled your sample program
and it works OK.

The problem is that my knowledge regarding geographical calculations and
projections is quite limited, so it’s going to take me some time to understand
the code.

My goal is to check if I can implement by myself other calculations like for
example: distance between a location and a geodesic line or check if three
locations are aligned. I’ll get to you with my conclusions.

Best Regards,
Gonzalo

[Charles Karney]

For more information on how to do these problems see Sections 11
and 13 of

C. F. F. Karney,
Geodesics on an ellipsoid of revolution,
Feb. 2011, http://arxiv.org/abs/1102.1215

The code below is a simple solution to the interception problem
using the gnomonic projection.

The problem is as follows:

A plane is traveling from A1 to A2. A missile battery is
located at B1. Find the point B2 on the path of the plane, such
that the distance from B1 to B2 is minimum.

The solution is similar to the intersection problem:

guess a position B2; make this the center of a gnomonic
projection; map A1, A2, B1 to this projection; solve the planar
interception problem to obtained an updated position B2; repeat.
The final azimuths printed allow you to verify that the
interception angle is 90 deg.

Example:

A1 = Istanbul, A2 = Washington, B1 = Reyjkavik

g++ -I/usr/local -L/usr/local -Wl,"-rpath=/usr/local/lib" \
-lGeographic -o intercept intercept.cpp

echo 42 29 39 -77 64 -22 | ./intercept

-->

Initial guess 40.5 -24
Increment 14.43198468492627 2.186112681554786
Increment -0.003453185661449254 -0.1234026147103826
Increment -2.147899635929207e-09 -1.132893178379391e-06
Increment -1.4210854715202e-14 3.552713678800501e-15
Increment 7.105427357601002e-15 -3.552713678800501e-15
Increment -7.105427357601002e-15 0
Increment 1.4210854715202e-14 7.105427357601002e-15
Increment -1.4210854715202e-14 0
Increment 7.105427357601002e-15 -7.105427357601002e-15
Increment -7.105427357601002e-15 0
Final result 54.92853149711691 -21.93729106604878
Azimuth to A1 89.82530404926784
Azimuth to A2 -90.17469595073216
Azimuth to B1 -0.1746959507321521

intercept.cpp is attached

[Anonymous]

It a is very useful tool, but which is the way to perform intersection of two geodesics in a similar manner to FAAs Compsys?

[Charles Karney]

I'm not sure what you're asking for here. If you want a slick
GUI such as Compsys offers, then I can't help.

[Anonymous]

Ok, Compsys is an achievement of many people, but your achievement is also remarkable by alone person.

I do not understand the problem, that is why i can not do the code or follow the instructions you gave.

My question is to intersect two geodesic path, which is the first step? I suppose to calculate first the two path then transform into gnomonic projection and finally calculate the intersection point. It is that correct?

[Charles Karney]

The algorithm is described in section 8 of "Algorithms for geodesics". I hope this provides you with the information you need.

[Anonymous]

The solution for calculating the intersection always returns the nearest intersection point without taking directionality into account. But what if you need the antipodal intersection point?
Example: plane A travels in a straight line from N52 E15 with an initial bearing of 80 degrees; plane B travels in a straight line from N51 E15 with an initial bearing of 110 degrees. Where do they collide?
The solution above returns N51.41228, E10.46628, but that would only be true when the planes traveled in the opposite directions. Instead, it should be something like S51.4117, W169.9028 (approximately; determined visually)

[Charles Karney]

For each plane, pick two points (1000km apart, let's say) on its path
which straddle the intended intersection point (visually determined as
you say). (These points can be computed by solving the direct geodesic
problem with distances of around 20000km.) Then run the algorithm with
these points.

Here's what I get

$ echo N52 E15 80 19e6 | GeodSolve 
-49.63686773 -178.95752685 110.55708754

$ echo N52 E15 80 20e6 | GeodSolve 
-52.01314639 -165.24507698 99.90432510

$ echo N51 E15 110 19e6 | GeodSolve 
-53.23887237 -179.40552276 81.11117316

$ echo N51 E15 110 20e6 | GeodSolve 
-50.97579152 -165.25227608 69.91829965

$ echo -49.63686773 -178.95752685 \
       -52.01314639 -165.24507698 \
       -53.23887237 -179.40552276 \
       -50.97579152 -165.25227608 | ./intersect
Initial guess -51.4661695025 -172.2151006675
Increment -0.1705621047736869 3.87998172225744
Increment 4.762136285307861e-08 -2.013739219819399e-07
Increment 0 0
...
Final result -51.63673155965233 -168.3351191466165
Azimuths on line A 102.3337261784672 102.3337261784673
Azimuths on line B 72.32472232694377 72.32472232694394

So the answer is S51.63673 W168.33512

[Anonymous]

Thank you very much!
I've incorporated it on this page: http://www.javawa.nl/coordcalc_en.html
It's a tool for projecting points, calculate distances and intersections etc.
(BTW: the coordinates in my example weren't correct, because I used a bearing of 100 degrees instead of 110)

[Anonymous]

On the 2nd question (about the distance from a point to a geodesic), using the 2nd provided source file intercept.cpp, i downloaded the last version of the library and installed it in std place (/usr/local, w.o. errors, tests succeeded too)

/usr/local/lib: ldconfig -p | grep -i geographic
libGeographic.so.19 (libc6,x86-64) => /usr/local/lib/libGeographic.so.19
libGeographic.so (libc6,x86-64) => /usr/local/lib/libGeographic.so
/usr/local/lib:

but when i try to compile the older code in intercept.cpp from
https://sourceforge.net/p/geographiclib/discussion/1026621/thread/21aaff9f/8a93/attachment/intercept.cpp
i get some linking errors (i guess):

$ g++ intercept.cpp
/tmp/ccQGg0bv.o: In function main': intercept.cpp:(.text+0xbe): undefined reference toGeographicLib::Geodesic::Geodesic(double, double)'
intercept.cpp:(.text+0xd7): undefined reference to GeographicLib::Gnomonic::Gnomonic(GeographicLib::Geodesic const&)' /tmp/ccQGg0bv.o: In functionGeographicLib::Geodesic::Inverse(double, double, double, double, double&, double&) const':
intercept.cpp:(.text.ZNK13GeographicLib8Geodesic7InverseEddddRdS1[ZNK13GeographicLib8Geodesic7InverseEddddRdS1]+0x8f): undefined reference to GeographicLib::Geodesic::GenInverse(double, double, double, double, unsigned int, double&, double&, double&, double&, double&, double&, double&) const' /tmp/ccQGg0bv.o: In functionGeographicLib::Gnomonic::Forward(double, double, double, double, double&, double&) const':
intercept.cpp:(.text.ZNK13GeographicLib8Gnomonic7ForwardEddddRdS1[ZNK13GeographicLib8Gnomonic7ForwardEddddRdS1]+0x77): undefined reference to GeographicLib::Gnomonic::Forward(double, double, double, double, double&, double&, double&, double&) const' /tmp/ccQGg0bv.o: In functionGeographicLib::Gnomonic::Reverse(double, double, double, double, double&, double&) const':
intercept.cpp:(.text.ZNK13GeographicLib8Gnomonic7ReverseEddddRdS1[ZNK13GeographicLib8Gnomonic7ReverseEddddRdS1]+0x77): undefined reference to GeographicLib::Gnomonic::Reverse(double, double, double, double, double&, double&, double&, double&) const' collect2: error: ld returned 1 exit status

I am trying to see if i can get the c++ solution to work, if it's not too much trouble.

[Charles Karney]

The comments at the top of intercept.cpp say to compile with

g++ -o intercept -I/usr/local intercept.cpp -lGeographic -Wl,-rpath=/usr/local/lib

Please report what you get when you do this.

[Anonymous]

Thank you very much, it solved the problem.

[ the compile command i tried 1st was from same 2nd question, "g++ -I/usr/local -L/usr/local -Wl,"-rpath=/usr/local/lib" -lGeographic -o intercept intercept.cpp". Then i tried without specifying include and linkage paths (since install was std and system found it) and got the same error, save for temp object files' names. The problems was the order of the options; it compiled now with what you just gave me, as well as with "g++ -o intercept2 intercept.cpp -lGeographic" and both produced the same final result as years ago on your example input:

$ echo 42 29 39 -77 64 -22 | ./intercept2
Initial guess 40.5 -24
Increment 14.43198468492627 2.186112681554786
Increment -0.003453185661449254 -0.1234026147103826
Increment -2.147906741356564e-09 -1.132893178379391e-06
Increment 0 0
Increment 0 0
Increment 0 0
Increment 0 0
Increment 0 0
Increment 0 0
Increment 0 0
Final result 54.92853149711691 -21.93729106604878
Azimuth to A1 89.82530404926784
Azimuth to A2 -90.17469595073216
Azimuth to B1 -0.1746959507321525
]

Thank you very, very much!!

[Anonymous]

I've written a Java implementation of the point to a line equation based on the provided cpp file above. However, after testing, it appears to be treating the line as a great circle and not as a geodesic. The closest points calculated are correct, but do not lay within the bounds of the geodesic. Is there a way to calculate the distance from a point to a geodesic? Or is that what is written above, and my implementation is flawed.

//Points A and B make up the geodesic (GeodesicData geo1), point C is the point part of point to a line

GnomonicData geo1A =
            gnom.Forward(estimatedCenter.latitude(), estimatedCenter.longitude(), geo1.lat1, geo1.lon1);
        GnomonicData geo1B =
            gnom.Forward(estimatedCenter.latitude(), estimatedCenter.longitude(), geo1.lat2, geo1.lon2);

        GnomonicData geo2C =
            gnom.Forward(estimatedCenter.latitude(), estimatedCenter.longitude(), distancedPoint.latitude(),
                distancedPoint.longitude());

        Vector3d vA1 = new Vector3d(geo1A.x, geo1A.y, 1);
        Vector3d vB1 = new Vector3d(geo1B.x, geo1B.y, 1);

        Vector3d la = new Vector3d();
        la.cross(vA1, vB1);

        Vector3d va1 = new Vector3d(geo2C.x, geo2C.y, 1);

        Vector3d lb = new Vector3d(la.y, (la.x *-1), (la.x * geo2C.y) - (la.y * geo2C.x));

        Vector3d p0 = new Vector3d();
        p0.cross(la, lb);
        norm(p0);

        GnomonicData intersectionData =
            gnom.Reverse(estimatedCenter.latitude(), estimatedCenter.longitude(), p0.x, p0.y);

[Charles Karney]

Please try with the example given above. If you get different answers then your Java code is behaving differently from the cpp code. Narrow down where this difference is occurring. Eventually you will find your problem!

[Anonymous]

Hi, thanks for the response. I've compared results of the provided point to a line cpp file and my java implementation, and I am getting the same results.

for the following input

 38.960844, -77.091494 Point A on geodesic
38.956638, -77.090743 Point B on geodesic
38.95617176919712, -77.09123611450195 Point C, distanced

The result is 38.95623330904161, -77.09067074538673

When plugging this into google maps or the mapping software of choice, the returned point doesn't lay on the geodesic., but on the great circle that the provided geodesic is a part of.

Is there a way to calculate it for solely the geodesic itself?

Thanks for the help!

[Charles Karney]

If you declare that the intercept point must lie within AB, and if the algorithm returns a point outside AB, then just return whichever of A and B is closer to C.

(By the way, your explanation is confusing because you incorrectly use the term "great circle". This lead me to believe that you were getting results which applied to a spherical, rather than an ellipsoidal, earth.)

[Anonymous]

Ok, sounds good. Thanks for the help! Sorry about my explanation being confusing, my quick googling led me astray :D, what would be the correct term instead of great circle?

[Charles Karney]

In contexts where you need to distinguish a bounded portion of a geodesic AB from the geodesic continued beyond AB, I recommend using "geodesic segment" for the former and "geodesic" for the latter. "Great circle" should be used only when talking about a spherical earth.

[Anonymous]

Ok, thanks for the help and the information! I really appreciate the quick responses, and the library is pretty great as well!

[Anonymous]

What about rhumb line . Are there any solutions to the same kind of problems for rhumb line in Geographiclib?

[Anonymous]

the link for intercept.cpp seems to be broken. Is there another way to get c copy of this file ?

[Charles Karney]

I posted the code below (you might have to go to the next page).

[Charles Karney]

The intersection of 2 rhumb lines reduces to a simple planar 2d problem (rhumb lines are straight lines in the Mercator projection). For the interception problem, you would need to specify how the closest point is determined (via a geodesic, rhumb line, etc.) and then take it from there. Even if the path is a rhumb line, you have the problem that in general the shortest path won't intersect the original rhumb line at right angles.