Re: [Geographiclib-users] Find S12 using planimeter

[Charles K. <cha...@sr...>]

You raise four issues:

(A) The cause of the discrepancy between GeodSolve's area and that
computed by Planimeter is as follows:

GeodSolve -f reports the area of the quadrilateral

     AF along a meridian
     FH along the equator
     HB along a meridian
     BA along a geodesic

In contrast, Planimeter reports the area of the quadrilateral

     AF along the shortest geodesic
     FH along the shortest geodesic
     HB along the shortest geodesic
     BA along the shortest geodesic

The geodesics AF and HB will be along meridians (for oblate ellipsoids).
However the geodesic HB only follows the equator if the longitude
difference is less than (1-f)*180 deg = 179.396494 deg.  In your
example, it picks a shorter geodesic through the northern hemisphere.
In fact, there are two equally short geodesics and it might be picked
the other one.

So, here's what the man page for Planimeter has to say about this:

        The edges of the polygon are given by the shortest geodesic
        between consecutive vertices.  In certain cases, there may be two
        or many such shortest geodesics, and in that case, the polygon is
        not uniquely specified by its vertices.  This only happens with
        very long edges (for the WGS84 ellipsoid, any edge shorter than
        19970 km is uniquely specified by its end points).  In such
        cases, insert an additional vertex near the middle of the long
        edge to define the boundary of the polygon.

(B) As far as getting the same results with the C library as with
GeodSolve,

Call either geod_geninverse or geod_gendirect.  As you can see from

   http://geographiclib.sourceforge.net/html/C/geodesic_8h.html

These take an optional pointer argument pS12 which returns the same area
as reported by GeodSolve.

(C) Do you need to worry about simple/non-simple polygons?  Definitely.
Areas for non-simple polygons can be defined in multiple ways and none
of them is the "right" way.

For simple polygons (not self intersecting), the area reported is
positive for counter-clockwise traversal.  If this area exceeds half the
area of the ellipsoid, the negative of the "rest" of the area of the
ellipsoid is returned.  (The -s and -r flags to Planimeter govern this
behavior.)  The poles of the ellipsoid are taken care of for you.

For non-simple polygons, the simple calculus definition is applied.  For
for the bow tie example, the triangle which is included
counter-clockwise counts as a positive area and the other one as a
negative area and the sum will be close to zero for a symmetric bow tie.
Similarly the area a polygon which encircles a point twice will be two
times the area of the same sized polygon which encircle it one.

If you want different rules applied for non-simple polygons, then you
must decompose such polygons into the union of simple polygons.  E.g.,
change the bow tie quadrilateral

     1 -1
     -1 -1
     1 1
     -1 1

to the hexagon

     1 -1
     -1 -1
     0 0
     -1 1
     1 1
     0 0

(D) Do you needs to worry about the poles?  No, not if you're using
Planimeter.  You can supply the boundary of Antarctica to Planimeter and
get the correct area.

I hope this answers your questions.

   --Charles

On 2014-09-18 07:03, Mike Toews wrote:
> Hi,
>
> I'm trying to use the C library (v1.32) to find polygon areas, but I'm
> having difficulty getting consistent results with the GeodTest.dat
> data set. From what I understand, to find the area under the geodesic,
> S12, the AFHB quadrilateral can be built with corners (lat1,lon1),
> (0,lon1), (0,lon2), and (lat2,lon2).
>
> Example 1 is line 1 of GeodTest.dat
>
> $ sed -n 1p GeodTest.dat
> 36.530042355041 0 176.125875162171 -48.164270779097768864
> 5.762344694676510456 175.334308316285410561 9398502.0434687
> 84.663858149358862201 6333544.7732452481809 -559418252332.321555
>
> or input for planimeter:
>
> 36.530042355041 0
> 0 0
> 0 5.762344694676510456
> -48.164270779097768864 5.762344694676510456
>
> With output:
> 4 19421014.749204 -559418252332.3
>
> This is great, since the areas are within the expected precision.
> However, do note that this quadrilateral is not simple, since it has a
> bow-tie, which is understood not be used by this function. Similar
> success with absolute residuals <= 0.1 m^2 is found from 55.3498% of
> the 500k lines in the dataset, even with simple quadrilaterals like
> line 1270.
>
> Example 2 is line 381 of GeodTest.dat, which forms a simple quadrilateral
>
> $ sed -n 381p GeodTest.dat
> 42.496576764546 0 .165548738367 61.200095175201199603
> 179.666196014678657282 179.746893753529690115 8509305.8221784
> 76.480111135775561028 6211557.6845057436865 127219347433402.487683
>
> or input for planimeter:
>
> 42.496576764546 0
> 0 0
> 0 179.666196014678657282
> 61.200095175201199603 179.666196014678657282
>
> With output:
> 4 39997539.52627776 47480061825561.484
>
> This isn't too great, since the absolute residual of area calculations
> is 79739285607841. or a ratio of about 0.373
>
> Now I know that GeographicLib (v1.37) can get the correct S12 for this
> example to GeodTest.dat, e.g.:
>
> $ ./GeodSolve -i -f
> 42.496576764546 0 61.200095175201199603 179.666196014678657282
> 42.49657676 0.00000000 0.16554874 61.20009518 179.66619601
> 179.74689375 8509305.822 76.48011114 6211557.685 0.2368660383
> 0.2380528208 127219347433402
>
> However, with the C++ tool Planimeter from this same version (v1.37),
> I get similar results to the C library. The only difference between
> the two is the precision of output, which is really just a formatting
> difference. So I suspect I'm doing my area under the geodesic
> calculation with the planimeter utilities incorrectly.
>
> My question is how can I achieve the same answer to GeodSolve with the
> C library? Also, do I need to worry about simple/non-simple polygons?
> Is there any further post-analysis if the polygon encircles a pole (in
> regards to last paragraph of section 6 of "Algorithms for geodesics"
> paper), or does the library handle these situations seamlessly?
>
> Many thanks,
>
> -Mike


-- 
Charles Karney <cha...@sr...>
SRI International, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662