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
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
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
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
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
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.
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?
The algorithm is described in section 8 of "Algorithms for geodesics". I hope this provides you with the information you need.
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)
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
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