[Geographiclib-users] Find S12 using planimeter

[Mike T. <mw...@gm...>]

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