Intersection between two geodesic lines

  • 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,

  • Charles Karney
    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

    Last edit: Charles Karney 2014-11-09
  • 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,

  • Charles Karney
    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,

    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.


    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

    Last edit: Charles Karney 2014-11-09

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

    Last edit: Charles Karney 2014-02-14

  • Anonymous

    Thank you very much!
    I've incorporated it on this page:
    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)



Cancel   Add attachments