geographiclib-users

[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

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

Re: [Geographiclib-users] Find S12 using planimeter

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

On 2014-09-18 09:54, Charles Karney wrote:
> However the geodesic HB only follows the equator if the longitude
> difference is less than (1-f)*180 deg = 179.396494 deg.

Sorry, this should be:

However the geodesic FH only follows the equator if the longitude
difference is less than (1-f)*180 deg = 179.396494 deg.

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

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

Re: [Geographiclib-users] Find S12 using planimeter

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

Great, thanks, I think that clears up my questions with my assumptions
on the quadrilateral for S12. And yes, the absolute error of residual
with the technique described previously increases approaching the
(1-f)*180 deg edges.

The planimeter approach is what I'm interested in. Out of curiosity,
are there any test datasets that can be used to check the accuracy of
polygon area calculations? Also, how does the error of polygon area
calculation accumulate? For instance, I see the expected error of S12
is 0.1 m^2. Does the error relate to the number of vertices in the
polygon in combination to the length of each edge? I just looking for
a general idea, nothing exhaustive.

Thanks in advance,
-Mike

Re: [Geographiclib-users] Find S12 using planimeter

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

The errors you get will very much depend on the type of polygon.  0.1m^2
per edge is a worst case; it appears that the errors cancel to some
extent.

Version 1.37 allows you to compile GeographicLib with Boost or MPFR
allowing you to use either quad or arbitrary precision.  So perhaps you
can gauge the errors in your case yourself.  (I can possibly get you a
binary version of GeographicLib + MPFR for Windows.)  If you are on a
Linux machine you can also compile using long doubles and get an extra
11 bits of precision.

As a practical example...  I have a dataset of the boundaries of Poland
which has

     67801 edges
     mean edge 53m
     min edge 0.06m
     max edge 10km

Planimeter gives the area as 312679715911.98608 m^2.  The true area is
312679715911.9856346137 m^2.  So the error in this case is less than
0.001 m^2.

   --Charles

On 2014-09-18 15:39, Mike Toews wrote:
> Great, thanks, I think that clears up my questions with my assumptions
> on the quadrilateral for S12. And yes, the absolute error of residual
> with the technique described previously increases approaching the
> (1-f)*180 deg edges.
>
> The planimeter approach is what I'm interested in. Out of curiosity,
> are there any test datasets that can be used to check the accuracy of
> polygon area calculations? Also, how does the error of polygon area
> calculation accumulate? For instance, I see the expected error of S12
> is 0.1 m^2. Does the error relate to the number of vertices in the
> polygon in combination to the length of each edge? I just looking for
> a general idea, nothing exhaustive.
>
> Thanks in advance,
> -Mike
>

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

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

Re: [Geographiclib-users] Find S12 using planimeter

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

On 19 September 2014 01:54, Charles Karney <cha...@sr...> wrote:
> 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

But of course finding the intersection point (i.e. 0 0) to make the
valid geometry is something altogether different! And yes, I see the
basic approach using a Gnomic projection. What I was looking to see if
computing an area on a non-simple polygon provides a similar answer as
other algorithms, and that it doesn't catch an unexpected exception
(e.g. divide by zero).

On 19 September 2014 10:27, Charles Karney <cha...@sr...> wrote:
> The errors you get will very much depend on the type of polygon.  0.1m^2
> per edge is a worst case; it appears that the errors cancel to some
> extent.
>
> Version 1.37 allows you to compile GeographicLib with Boost or MPFR
> allowing you to use either quad or arbitrary precision.  So perhaps you
> can gauge the errors in your case yourself.  (I can possibly get you a
> binary version of GeographicLib + MPFR for Windows.)  If you are on a
> Linux machine you can also compile using long doubles and get an extra
> 11 bits of precision.
>
> As a practical example...  I have a dataset of the boundaries of Poland
> which has
>
>     67801 edges
>     mean edge 53m
>     min edge 0.06m
>     max edge 10km
>
> Planimeter gives the area as 312679715911.98608 m^2.  The true area is
> 312679715911.9856346137 m^2.  So the error in this case is less than
> 0.001 m^2.

These errors are small enough to me. I might expect the largest errors
to be for simplified polygons with few vertices and many sharp angles.

I'm unable to build a higher precision version. On Debian Wheezy, I
can't compile with either GEOGRAPHICLIB_PRECISION 4 or 5, since the
supported releases of Boost 1.49 or MPFR 3.1.0 don't meet the
requirements for v 1.37 of GeographicLib, which requires either Boost
1.54 or MPFR 3.5.9 (in CMakeLists.txt). So I'll test this out some
other time.

-Mike

Re: [Geographiclib-users] Find S12 using planimeter

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

On 2014-09-25 09:01, Mike Toews wrote:
> I'm unable to build a higher precision version. On Debian Wheezy, I
> can't compile with either GEOGRAPHICLIB_PRECISION 4 or 5, since the
> supported releases of Boost 1.49 or MPFR 3.1.0 don't meet the
> requirements for v 1.37 of GeographicLib, which requires either Boost
> 1.54 or MPFR 3.5.9 (in CMakeLists.txt). So I'll test this out some
> other time.
>

You've got the MPFR requirements wrong...

I'm using 3.1.2 but I would expect 3.1.0 to be OK too.  The 3.5.9 refers
to MPFR C++ which is a single header file (mpreal.h) you can stick into
the include directory.  The MPFR C++ download link is

   http://www.holoborodko.com/pavel/mpfr/#download

I will fix the documentation in CMakeLists.txt and in doxygen to clarify
this.

By the way, at some point as you increase the precision you need to
switch to using Planimeter with the -E option.  The geodesics are then
computed with elliptic integrals.  The area calculation uses a 30th
order series (I never figured out how to express the are integral in
terms of elliptic integrals); so eventually you will run out of accuracy
on areas (but it's probably good for at least 1000 digits for WGS84).

   --Charles

Re: [Geographiclib-users] Find S12 using planimeter

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

On 26 September 2014 03:56, Charles Karney <cha...@sr...> wrote:
> You've got the MPFR requirements wrong...
>
> I'm using 3.1.2 but I would expect 3.1.0 to be OK too.  The 3.5.9 refers
> to MPFR C++ which is a single header file (mpreal.h) you can stick into
> the include directory.
>

OK thanks, I see the distinction between MPFR and MPFR C++; I had only
installed the former. The later is not yet available in Debian stable,
but is in testing (Debian Jessie). Regardless, it's a simple header
file that can easily be downloaded and placed in /usr/local/include,
and the project compiles normally.

-Mike

Re: [Geographiclib-users] Find S12 using planimeter

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

On 09/25/2014 04:40 PM, Mike Toews wrote:
> On 26 September 2014 03:56, Charles Karney <cha...@sr...> wrote:
>> You've got the MPFR requirements wrong...
>>
>> I'm using 3.1.2 but I would expect 3.1.0 to be OK too.  The 3.5.9 refers
>> to MPFR C++ which is a single header file (mpreal.h) you can stick into
>> the include directory.
>>
>
> OK thanks, I see the distinction between MPFR and MPFR C++; I had only
> installed the former. The later is not yet available in Debian stable,
> but is in testing (Debian Jessie). Regardless, it's a simple header
> file that can easily be downloaded and placed in /usr/local/include,
> and the project compiles normally.
>
> -Mike
>

The way to set the MPFR precision in the utility programs is as follows:

   echo '0 0;0 90;90 0' |
   tr ';' '\n' |
   GEOGRAPHICLIB_DIGITS=1000 tools/Planimeter -p 100

which means calculate with 1000 bits of precision and print the results
with ~100 decimal digits of precision.

   --Charles