Boost C++ Libraries: Ticket #11917: geometry default geographic (andoyer) antipodal distance is zero https://svn.boost.org/trac10/ticket/11917 <p> When calling the default (andoyer) <code>boost::geometry::distance()</code> function, with <code>boost::geometry::cs::geographic</code> points that are close to being antipodal (e.g. the North and South poles or points that are diametrically opposed around the equator), it returns reasonable results. However, when called with antipodal points it returns zero, e.g.: </p> <pre class="wiki">Andoyer distance near poles: 20003917.164970 near opposite equatorial points: 20037508.151445 Andoyer distance at poles: 0.000000 opposite equatorial points: 0.000000 </pre><p> Here the code I used to obtain the results above: </p> <div class="wiki-code"><div class="code"><pre> <span class="n">std</span><span class="o">::</span><span class="n">cout</span><span class="p">.</span><span class="n">setf</span><span class="p">(</span><span class="n">std</span><span class="o">::</span><span class="n">ios</span><span class="o">::</span><span class="n">fixed</span><span class="p">);</span> <span class="k">typedef</span> <span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="o">::</span><span class="n">cs</span><span class="o">::</span><span class="n">geographic</span><span class="o">&lt;</span><span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="o">::</span><span class="n">radian</span><span class="o">&gt;</span> <span class="n">Wgs84Coords</span><span class="p">;</span> <span class="k">typedef</span> <span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="o">::</span><span class="n">model</span><span class="o">::</span><span class="n">point</span><span class="o">&lt;</span><span class="kt">double</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">Wgs84Coords</span><span class="o">&gt;</span> <span class="n">GeographicPoint</span><span class="p">;</span> <span class="c1">// Min value for delta. 0.000000014 causes Andoyer to fail.</span> <span class="k">const</span> <span class="kt">double</span> <span class="nf">DELTA</span><span class="p">(</span><span class="mf">0.000000015</span><span class="p">);</span> <span class="c1">// Note boost points are Long &amp; Lat NOT Lat &amp; Long</span> <span class="n">GeographicPoint</span> <span class="nf">near_north_pole</span> <span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="n">M_PI_2</span> <span class="o">-</span> <span class="n">DELTA</span><span class="p">);</span> <span class="n">GeographicPoint</span> <span class="nf">near_south_pole</span> <span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="o">-</span><span class="n">M_PI_2</span> <span class="o">+</span> <span class="n">DELTA</span><span class="p">);</span> <span class="kt">double</span> <span class="nf">andoyer_minor</span><span class="p">(</span><span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="o">::</span><span class="n">distance</span><span class="p">(</span><span class="n">north_pole</span><span class="p">,</span> <span class="n">south_pole</span><span class="p">));</span> <span class="n">std</span><span class="o">::</span><span class="n">cout</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;Andoyer distance near poles: &quot;</span> <span class="o">&lt;&lt;</span> <span class="n">andoyer_minor</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;</span><span class="se">\n</span><span class="s">&quot;</span><span class="p">;</span> <span class="n">GeographicPoint</span> <span class="nf">near_equator_east</span> <span class="p">(</span> <span class="n">M_PI_2</span> <span class="o">-</span> <span class="n">DELTA</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">);</span> <span class="n">GeographicPoint</span> <span class="nf">near_equator_west</span> <span class="p">(</span><span class="o">-</span><span class="n">M_PI_2</span> <span class="o">+</span> <span class="n">DELTA</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">);</span> <span class="kt">double</span> <span class="nf">andoyer_major</span><span class="p">(</span><span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="o">::</span><span class="n">distance</span><span class="p">(</span><span class="n">near_equator_east</span><span class="p">,</span> <span class="n">near_equator_west</span><span class="p">));</span> <span class="n">std</span><span class="o">::</span><span class="n">cout</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;near opposite equatorial points: &quot;</span> <span class="o">&lt;&lt;</span> <span class="n">andoyer_major</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;</span><span class="se">\n\n</span><span class="s">&quot;</span><span class="p">;</span> <span class="n">GeographicPoint</span> <span class="nf">north_pole</span> <span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="n">M_PI_2</span><span class="p">);</span> <span class="n">GeographicPoint</span> <span class="nf">south_pole</span> <span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="o">-</span><span class="n">M_PI_2</span><span class="p">);</span> <span class="n">andoyer_minor</span> <span class="o">=</span> <span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="o">::</span><span class="n">distance</span><span class="p">(</span><span class="n">north_pole</span><span class="p">,</span> <span class="n">south_pole</span><span class="p">);</span> <span class="n">std</span><span class="o">::</span><span class="n">cout</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;Andoyer distance at poles: &quot;</span> <span class="o">&lt;&lt;</span> <span class="n">andoyer_minor</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;</span><span class="se">\n</span><span class="s">&quot;</span><span class="p">;</span> <span class="n">GeographicPoint</span> <span class="nf">equator_east</span> <span class="p">(</span> <span class="n">M_PI_2</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">);</span> <span class="n">GeographicPoint</span> <span class="nf">equator_west</span> <span class="p">(</span><span class="o">-</span><span class="n">M_PI_2</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">);</span> <span class="n">andoyer_major</span> <span class="o">=</span> <span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="o">::</span><span class="n">distance</span><span class="p">(</span><span class="n">equator_east</span><span class="p">,</span> <span class="n">equator_west</span><span class="p">);</span> <span class="n">std</span><span class="o">::</span><span class="n">cout</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;opposite equatorial points: &quot;</span> <span class="o">&lt;&lt;</span> <span class="n">andoyer_major</span> <span class="o">&lt;&lt;</span> <span class="n">std</span><span class="o">::</span><span class="n">endl</span><span class="p">;</span> </pre></div></div><p> This is clearly a bug. The distance between the antipodal points should be slightly larger than the distance between the other points, and definitely NOT zero! </p> en-us Boost C++ Libraries /htdocs/site/boost.png https://svn.boost.org/trac10/ticket/11917 Trac 1.4.3 Charles Karney <charles@…> Sun, 17 Jan 2016 16:01:11 GMT <link>https://svn.boost.org/trac10/ticket/11917#comment:1 </link> <guid isPermaLink="false">https://svn.boost.org/trac10/ticket/11917#comment:1</guid> <description> <p> I recommend that boost::geometry include the algorithm for geodesic distance given in my paper, <a class="ext-link" href="https://dx.doi.org/10.1007/s00190-012-0578-z"><span class="icon">​</span>Algorithms for geodesics</a>. This returns an accurate result for the distance for arbitrary points on terrestrially relevant ellipsoids. <a class="ext-link" href="http://geographiclib.sourceforce.net"><span class="icon">​</span>GeographicLib</a> (version 1.25 and later) also includes the solution in terms of elliptic integrals which is accurate for ellipsoids of revolution with 1/100 &lt; b/a &lt; 100. </p> </description> <category>Ticket</category> </item> <item> <author>Ken Barker <ken.barker@…></author> <pubDate>Sun, 17 Jan 2016 19:20:16 GMT</pubDate> <title/> <link>https://svn.boost.org/trac10/ticket/11917#comment:2 </link> <guid isPermaLink="false">https://svn.boost.org/trac10/ticket/11917#comment:2</guid> <description> <p> I agree with Charles Karney's comment above. </p> <p> I originally raised this issue after finding an issue with the <code>Vincenty</code> function, see: <a class="ext-link" href="http://stackoverflow.com/questions/34767736/why-is-boostgeometry-geographic-vincenty-distance-inaccurate-around-the-equato"><span class="icon">​</span>http://stackoverflow.com/questions/34767736/why-is-boostgeometry-geographic-vincenty-distance-inaccurate-around-the-equato</a> </p> <p> On further investigation (and with a bit of help from Charles) I came to realise that both the <code>Vincenty</code> and <code>Andoyer</code> algorithms calculate incorrect distances for antipodal points on the equator. Whilst [<a class="missing wiki">GeographicLib</a>](<a class="ext-link" href="http://geographiclib.sourceforge.net/"><span class="icon">​</span>http://geographiclib.sourceforge.net/</a>) calculates the correct distance. </p> <p> It should be included in boost geometry and replace <code>Andoyer</code> as the default algorithm. </p> <p> BTW, please note that Charles's link to <code>GeographicLib</code> is incorrect, please use the one in this comment. </p> </description> <category>Ticket</category> </item> <item> <author>Charles Karney <charles@…></author> <pubDate>Sun, 17 Jan 2016 20:15:11 GMT</pubDate> <title/> <link>https://svn.boost.org/trac10/ticket/11917#comment:3 </link> <guid isPermaLink="false">https://svn.boost.org/trac10/ticket/11917#comment:3</guid> <description> <p> It might also be worth pointing out that "Algorithms for Geodesics" also provides a method for computing the area of a geodesic polygon (a polygon whose sides are geodesics) on an ellipsoid. </p> </description> <category>Ticket</category> </item> <item> <dc:creator>Barend Gehrels</dc:creator> <pubDate>Mon, 18 Jan 2016 07:52:43 GMT</pubDate> <title/> <link>https://svn.boost.org/trac10/ticket/11917#comment:4 </link> <guid isPermaLink="false">https://svn.boost.org/trac10/ticket/11917#comment:4</guid> <description> <p> Thanks, both, and we welcome the papers and involvement of Charles Karney! We were already aware of your library and your papers, also with respect to area, this is very helpful. </p> </description> <category>Ticket</category> </item> <item> <dc:creator>awulkiew</dc:creator> <pubDate>Mon, 25 Jan 2016 18:55:11 GMT</pubDate> <title>status, version, milestone changed; resolution set https://svn.boost.org/trac10/ticket/11917#comment:5 https://svn.boost.org/trac10/ticket/11917#comment:5 <ul> <li><strong>status</strong> <span class="trac-field-old">new</span> → <span class="trac-field-new">closed</span> </li> <li><strong>version</strong> <span class="trac-field-old">Boost 1.61.0</span> → <span class="trac-field-new">Boost 1.60.0</span> </li> <li><strong>resolution</strong> → <span class="trac-field-new">fixed</span> </li> <li><strong>milestone</strong> <span class="trac-field-old">To Be Determined</span> → <span class="trac-field-new">Boost 1.61.0</span> </li> </ul> Ticket