https://svn.boost.org/trac10/ticket/11931 <p> The following code illustrates a problem with handling of meridian crossing and pole encircling when computing spherical areas. It computes the area of a 4 identical diamond shaped areas at different positions on the sphere. Only the area of the last of these is computed correctly. The others are badly wrong. </p> <div class="wiki-code"><div class="code"><pre><span class="cp">#include</span> <span class="cpf">&lt;iostream&gt;</span><span class="cp"></span> <span class="cp">#include</span> <span class="cpf">&lt;string&gt;</span><span class="cp"></span> <span class="cp">#include</span> <span class="cpf">&lt;boost/geometry.hpp&gt;</span><span class="cp"></span> <span class="k">namespace</span> <span class="n">bg</span> <span class="o">=</span> <span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="p">;</span> <span class="k">typedef</span> <span class="n">bg</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">bg</span><span class="o">::</span><span class="n">cs</span><span class="o">::</span><span class="n">spherical_equatorial</span><span class="o">&lt;</span><span class="n">bg</span><span class="o">::</span><span class="n">degree</span><span class="o">&gt;</span> <span class="o">&gt;</span> <span class="n">pt</span><span class="p">;</span> <span class="kt">int</span> <span class="nf">main</span><span class="p">()</span> <span class="p">{</span> <span class="n">std</span><span class="o">::</span><span class="n">string</span> <span class="n">polys</span><span class="p">[]</span> <span class="o">=</span> <span class="p">{</span> <span class="s">&quot; 0 89, 90 89, 180 89, 270 89&quot;</span><span class="p">,</span> <span class="s">&quot; 0 -89, -90 -89, -180 -89, -270 -89&quot;</span><span class="p">,</span> <span class="s">&quot;-1 0, 0 -1, 1 0, 0 1&quot;</span><span class="p">,</span> <span class="s">&quot;89 0, 90 -1, 91 0, 90 1&quot;</span> <span class="p">};</span> <span class="c1">// Correct area for all these polygons is</span> <span class="kt">double</span> <span class="n">area0</span> <span class="o">=</span> <span class="mf">0.00060926577032113655</span><span class="p">;</span> <span class="n">bg</span><span class="o">::</span><span class="n">model</span><span class="o">::</span><span class="n">polygon</span><span class="o">&lt;</span><span class="n">pt</span><span class="p">,</span> <span class="nb">false</span><span class="p">,</span> <span class="nb">false</span><span class="o">&gt;</span> <span class="n">poly</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; area relerr&quot;</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="k">for</span> <span class="p">(</span><span class="kt">int</span> <span class="n">i</span> <span class="o">=</span> <span class="mi">0</span><span class="p">;</span> <span class="n">i</span> <span class="o">&lt;</span> <span class="mi">4</span><span class="p">;</span> <span class="o">++</span><span class="n">i</span><span class="p">)</span> <span class="p">{</span> <span class="n">bg</span><span class="o">::</span><span class="n">read_wkt</span><span class="p">(</span><span class="s">&quot;POLYGON((&quot;</span> <span class="o">+</span> <span class="n">polys</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">+</span> <span class="s">&quot;))&quot;</span><span class="p">,</span> <span class="n">poly</span><span class="p">);</span> <span class="kt">double</span> <span class="n">area</span> <span class="o">=</span> <span class="n">bg</span><span class="o">::</span><span class="n">area</span><span class="p">(</span><span class="n">poly</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="n">area</span> <span class="o">&lt;&lt;</span> <span class="s">&quot; &quot;</span> <span class="o">&lt;&lt;</span> <span class="p">(</span><span class="n">area</span> <span class="o">-</span> <span class="n">area0</span><span class="p">)</span> <span class="o">/</span> <span class="n">area0</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="p">}</span> <span class="p">}</span> </pre></div></div><p> Running this code gives </p> <pre class="wiki"> area relerr 0.000304633 -0.5 -6.28288 -10313.2 0 -1 0.000609266 2.31338e-15 </pre><p> Discussion: </p> <ul><li>I can't make out the handling of this issue in the current code at all. It seems to be irredeemably broken. </li><li>You need to do two things to handle these issues: <ol><li>Make sure that the sign of lambda_2 - lambda_1 correctly corresponds to east and west going edges (with some consistent convention of edges which pass over the pole). </li><li>Count the number of times the polygon encircles a pole, e.g., by keeping track of when an edge crosses the prime meridian (again with some consistent treatment of cases when one or both ends of an edge lie on the prime meridian). </li></ol></li><li>If this count is odd, then the total area needs to be adjusted by 2*pi, half the area of the sphere. </li><li>I believe I've handled these issues correctly in the PolygonAreaT class in <a class="ext-link" href="http://geographiclib.sourceforge.net/html/"><span class="icon">​</span>GeographicLib</a>. </li></ul> en-us /htdocs/site/boost.png https://svn.boost.org/trac10/ticket/11931 Trac 1.4.3 awulkiew Wed, 29 Mar 2017 23:00:57 GMT https://svn.boost.org/trac10/ticket/11931#comment:1 https://svn.boost.org/trac10/ticket/11931#comment:1 <ul> <li><strong>keywords</strong> area added </li> <li><strong>status</strong> <span class="trac-field-old">new</span> → <span class="trac-field-new">closed</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.64.0</span> </li> </ul> <p> Thanks! </p> <p> Fixed by this PR: <a class="ext-link" href="https://github.com/boostorg/geometry/pull/359"><span class="icon">​</span>https://github.com/boostorg/geometry/pull/359</a> </p> <p> The new results are: </p> <pre class="wiki"> area relerr 0.000609266 -1.40992e-12 0.000609266 -1.40992e-12 0.000609266 1.77952e-16 0.000609266 1.42362e-15 </pre> Ticket