Boost point_circle coming out in weird shapes

Question

I'm attempting to create a 10m radius polygon on Earth using Boost's Geometry library.

Here's the tutorial.

To compile this example, I used Wandbox with the latest Clang and Boost 1.73.0.

I first discovered the issue in my production environment, which is Clang 12 and Boost 1.71.0.

Using a 1000m radius circle with 32 points yields expected results:

Shrinking that down to 10m however has unexpected results:

I used a WKT playground to display the results, and have confirmed that results are the same in other visualisation tools.

It seems to be a floating point rounding error, but everything here should be using double-precision floats which are more than enough to represent GPS coordinates. Something seems to be going wrong with the calculation.

The same thing happens with boost::geometry::point_circle using a radius of 0.0001.

What's going on, and should I just calculate the circle manually instead?

Edit 1

It gets even weirder if you use bg::area to calculate the area. I tried on a '10m' radius circle drawn around POINT(4.9 52.1) and got an area of 25984.4m. I tried the same at POINT(4.9 52.1000001) and got -1122.14.

See the following playground: https://godbolt.org/z/sTGqKK

Edit 2

I discovered that the issue with displaying the polygon is separate to the issue of the calculated area being incorrect. In fact, the display issue is as a result of rounding when printing to stdout. By increasing the precision of decimals, or using std::fixed, the display issue is resolved.

std::cout << std::fixed << bg::wkt(result) << std::endl;

Further details on how boost::geometry::buffer works: https://www.boost.org/doc/libs/master/libs/geometry/doc/html/geometry/reference/algorithms/buffer/buffer_7_with_strategies.html — Chris Watts, Dec 29 '20 at 14:49
If I read [this](https://sites.google.com/site/trescopter/Home/concepts/required-precision-for-gps-calculations#TOC-Precision-of-Float-and-Double) I don't think I see actual numbers for double precision, but the single precision accuracies list worst case 1.64-2.37m accuracy. If you divide a circle of radius 10m to that kind of precision, it's like drawing a circle r=5 on an integer grid (10m/~2m = 5). That would explain a shape like you show. That the "resolution is enough for GPS" is hardly relevant because you couldn't use GPS to accurately outline a 10m radius circle either. — sehe, Dec 29 '20 at 15:16
@sehe Yes, what we see is almost as if it's displaying a float. The machine epsilon for double is 2.22e-16, which according to [this table](https://en.wikipedia.org/wiki/Decimal_degrees#Precision) ought to have a scale precision of much less than 1mm (2.471304e-11 metres to be precise). — Chris Watts, Dec 29 '20 at 15:26

score 2 · Answer 1 · answered Dec 29 '20 at 17:20

It does seem that there are accuracy issues. I tried to work around things but didn't get as far as I'd like.

BGL uses some hard-qualified std::abs and std::acos calls that make it hard to use multiprecision types. I tried patching some of these, but the rabit hole was too deep for an afternoon.

Here's a testbed that might help pinpoint/debug/trace things a bit further. Note that

for float the accuracy is such that the library is_valid will report invalid due to spikes.
long double seems to do reasonable

The over-arching problem (lack of control/predictability) however remains.

Live On Compiler Explorer¹

#include <boost/geometry.hpp>
#include <iostream>

#ifdef TRY_BOOST_MULTIPRECISION
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
    namespace bmp = boost::multiprecision;
    using OctFloat    = bmp::cpp_bin_float_oct;
    using Decimal     = bmp::number<bmp::cpp_dec_float<50>,  bmp::et_off>;
    using LongDecimal = bmp::number<bmp::cpp_dec_float<100>, bmp::et_off>;

    namespace boost::multiprecision {
        inline auto mod(OctFloat    const& a, OctFloat    const& b) { return fmod(a, b); }
        inline auto mod(Decimal     const& a, Decimal     const& b) { return fmod(a, b); }
        inline auto mod(LongDecimal const& a, LongDecimal const& b) { return fmod(a, b); }
        inline auto abs(OctFloat    const& a) { return fabs(a); }
        inline auto abs(Decimal     const& a) { return fabs(a); }
        inline auto abs(LongDecimal const& a) { return fabs(a); }
    }

    namespace std { // sadly BG overqualifies std::abs in places
        inline auto abs(OctFloat    const& a) { return fabs(a); }
    }
#endif

template <typename F, typename DegreeOrRadian>
void do_test(int n, F offset = {}) {
    namespace bg = boost::geometry;
    std::cout << "----- " << __PRETTY_FUNCTION__ << " n:" << n << " offset: " << offset << " ----\n";
    bg::model::point<F, 2, bg::cs::geographic<bg::degree> > Amsterdam { 4.9, 52.1 + offset };
    typedef bg::model::point<F, 2, bg::cs::geographic<DegreeOrRadian> > point;

    // Declare the geographic_point_circle strategy (with n points)
    // Default template arguments (taking Andoyer strategy)
    bg::strategy::buffer::geographic_point_circle<> point_strategy(n);

    // Declare the distance strategy (one kilometer, around the point, on Earth)
    bg::strategy::buffer::distance_symmetric<F> distance_strategy(10.0);

    // Declare other necessary strategies, unused for point
    bg::strategy::buffer::join_round    join_strategy;
    bg::strategy::buffer::end_round     end_strategy;
    bg::strategy::buffer::side_straight side_strategy;

    // Declare/fill a point on Earth, near Amsterdam
    point p;
    bg::convert(Amsterdam, p);

    // Create the buffer of a point on the Earth
    bg::model::multi_polygon<bg::model::polygon<point> > result;
    bg::buffer(p, result,
                distance_strategy, side_strategy,
                join_strategy, end_strategy, point_strategy);

    std::string reason;
    is_valid(result, reason);
    //std::cout << "result: " << wkt(result) << "\n";
    std::cout << reason << "\n";
    std::cout << "result: " << (bg::is_simple(result)?"simple":"compound") << "\n";

    auto area = bg::area(result);

    std::cout << "reference: " << bg::dsv(Amsterdam)  << std::endl;
    std::cout << "point: " << bg::dsv(p)  << std::endl;
    std::cout << "area: " <<  area << " m²" << std::endl;
}

int main() {
    for (long double offset : { 0.l/*, 1e-7l*/ }) {
        for (int n : { 36 }) {
            do_test<float,       boost::geometry::degree>(n, offset);
            do_test<double,      boost::geometry::degree>(n, offset);
            do_test<long double, boost::geometry::degree>(n, offset);

            do_test<float,       boost::geometry::radian>(n, offset);
            do_test<double,      boost::geometry::radian>(n, offset);
            do_test<long double, boost::geometry::radian>(n, offset);

            // not working yet
            //do_test<OctFloat,    boost::geometry::radian>(n, offset);
            //do_test<Decimal,     boost::geometry::degree>();
            //do_test<LongDecimal, boost::geometry::degree>();
        }
    }
}

Prints

----- void do_test(int, F) [F = float, DegreeOrRadian = boost::geometry::degree] n:36 offset: 0 ----
Geometry has spikes. A spike point was found with apex at (4.9, 52.0975)
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: -1.37916e+07 m²
----- void do_test(int, F) [F = double, DegreeOrRadian = boost::geometry::degree] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 25984.4 m²
----- void do_test(int, F) [F = long double, DegreeOrRadian = boost::geometry::degree] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 301.264 m²
----- void do_test(int, F) [F = float, DegreeOrRadian = boost::geometry::radian] n:36 offset: 0 ----
Geometry has spikes. A spike point was found with apex at (-1.38318, -1.30708)
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 1.85308e+08 m²
----- void do_test(int, F) [F = double, DegreeOrRadian = boost::geometry::radian] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 6399.41 m²
----- void do_test(int, F) [F = long double, DegreeOrRadian = boost::geometry::radian] n:36 offset: 0 ----
Geometry is valid
result: simple
reference: (4.9, 52.1)
point: (4.9, 52.1)
area: 302.318 m²

On my machine

¹ exceeds processing time

A valiant effort so far, thanks for looking into this. As you say, it times out on Compiler Explorer, so I need to explore your test results further to verify. However, I have tried using `long double` in my production environment and I'm still getting wildly varying areas even with that. One thing I have noticed however - a clue - is that the in-memory representation of the circle polygon appears to be correct, but the WKT output is not. I manually plotted the outer ring of the polygon result (with num points = 8) and it was indeed an octagon. Does not explain the area issue though. — Chris Watts, Dec 29 '20 at 18:17
You might want to post on the [boost-geometry](https://lists.boost.org/mailman/listinfo.cgi/geometry) mailing list as well. @adamwulkiewicz @ barendgehrels @ mloskot and others are usually quite responsive there and have the GIS chops I lack to understand how these projections/coordinate systems are implemented. — sehe, Dec 29 '20 at 22:39
I've validated your findings, and written an error function to show how far out the areas are. Just waiting to hear back from the mailing list now. https://godbolt.org/z/sTGqKK — Chris Watts, Dec 30 '20 at 13:50
A bug report has been submitted for this, as I believe it's an avoidable underflow: https://github.com/boostorg/geometry/issues/799 The issue is such that if you compile on AArch64 (such as a smartphone, or even the new M1 Apple Macs) it's impossible to get a good result: `long double` doesn't appear to be supported. — Chris Watts, Jan 21 '21 at 23:17
Thanks for the update. I agree this seemed to be suboptimal. The ticket will be helpful for future visitors/reference! — sehe, Jan 21 '21 at 23:40

score 2 · Accepted Answer · answered Feb 08 '21 at 14:49

To my understanding there are two sources of inaccuracy, the area algorithm and the geographic buffer over a point algorithm.

Regarding the former https://github.com/boostorg/geometry/pull/801 proposes a fix. Using this fix the above error function (godbolt.org/z/sTGqKK) returns less than 1% relative error. The following piece of code extends this by using strategies.

#include <boost/geometry.hpp>
#include <cmath>
#include <iostream>
 
template <typename CT>
void error_function(CT area, CT theoreticalArea)
{
    std::cout << "area: " <<  area << " m², ";
    std::cout << "error: " <<  area - theoreticalArea << " m²,\t";
    std::cout << "normalised error: " <<  fabs(100 * (area - theoreticalArea)
        / theoreticalArea) << "%" << std::endl;
}

template <typename F, typename DegreeOrRadian>
void do_test(int n, F radius, F offset = {}) {
    namespace bg = boost::geometry;

    std::cout
        << "----- " << __PRETTY_FUNCTION__
        << " n:" << n << " radius:" << radius << " offset:" << offset
        << " ----\n";

    bg::model::point<F, 2, bg::cs::geographic<bg::degree> > Amsterdam { 4.9, 52.1 + offset };
    typedef bg::model::point<F, 2, bg::cs::geographic<DegreeOrRadian> > point;

    // Declare the geographic_point_circle strategy (with n points)
    // Default template arguments (taking Andoyer strategy)
    bg::strategy::buffer::geographic_point_circle<> point_strategy(n);

    // Declare the distance strategy (ten metres, around the point, on Earth)
    bg::strategy::buffer::distance_symmetric<F> distance_strategy(radius);

    // Declare other necessary strategies, unused for point
    bg::strategy::buffer::join_round    join_strategy;
    bg::strategy::buffer::end_round     end_strategy;
    bg::strategy::buffer::side_straight side_strategy;

    // Declare/fill a point on Earth, near Amsterdam
    point p;
    bg::convert(Amsterdam, p);

    // Create the buffer of a point on the Earth
    bg::model::multi_polygon<bg::model::polygon<point> > result;
    bg::buffer(p, result,
                distance_strategy, side_strategy,
                join_strategy, end_strategy, point_strategy);

    auto area = bg::area(result);
    auto areat = bg::area(result,bg::strategy::area::geographic<bg::strategy::thomas>());
    auto areav = bg::area(result,bg::strategy::area::geographic<bg::strategy::vincenty>());
    auto areak = bg::area(result,bg::strategy::area::geographic<bg::strategy::karney>());

    // Assumes that the Earth is flat, which it clearly is.
    // A = n/2 * R^2 * sin(2*pi/n) where R is the circumradius
    auto theoreticalArea = n * radius * radius * std::sin(2.0 * 3.142 / n) / 2.0;

    std::cout << "reference: " << bg::dsv(Amsterdam)  << std::endl;
    std::cout << "point: " << bg::dsv(p)  << std::endl;
    std::cout << "radius: " <<  radius << " m" << std::endl;
    error_function(area, theoreticalArea);
    error_function(areat, theoreticalArea);
    error_function(areav, theoreticalArea);
    error_function(areak, theoreticalArea);
}

int main() {
    double offset = 1e-7;
    int n = 8;

    do_test<double,      boost::geometry::degree>(n, 10.);
    do_test<long double, boost::geometry::degree>(n, 10.);

    do_test<double,      boost::geometry::radian>(n, 10.);
    do_test<long double, boost::geometry::radian>(n, 10.);

    do_test<double,      boost::geometry::degree>(n, 10., offset);
    do_test<long double, boost::geometry::degree>(n, 10., offset);

    do_test<double,      boost::geometry::degree>(n, 1000.);
    do_test<double,      boost::geometry::degree>(n, 1000., offset);

    do_test<double,      boost::geometry::degree>(n, 1.);
    do_test<long double, boost::geometry::degree>(n, 1.);
}

which returns:

----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 281.272 m²,   error: -1.59991 m²,     normalised error: 0.565596%
area: 282.843 m²,   error: -0.0284134 m²,   normalised error: 0.0100446%
area: 281.749 m²,   error: -1.12206 m²,     normalised error: 0.396666%
area: 282.843 m²,   error: -0.028415 m²,    normalised error: 0.0100452%
----- void do_test(int, F, F) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 283.57 m²,    error: 0.698736 m²,     normalised error: 0.247015%
area: 282.843 m²,   error: -0.0284201 m²,   normalised error: 0.010047%
area: 283.568 m²,   error: 0.696594 m²,     normalised error: 0.246258%
area: 282.843 m²,   error: -0.0284255 m²,   normalised error: 0.0100489%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::radian] n:8 radius:10 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 282.715 m²,   error: -0.156633 m²,    normalised error: 0.0553726%
area: 282.843 m²,   error: -0.0286857 m²,   normalised error: 0.0101409%
area: 280.578 m²,   error: -2.29311 m²,     normalised error: 0.810656%
area: 282.843 m²,   error: -0.0286896 m²,   normalised error: 0.0101423%
----- void do_test(int, F, F) [with F = long double; DegreeOrRadian = boost::geometry::radian] n:8 radius:10 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 283.135 m²,   error: 0.263058 m²,     normalised error: 0.0929955%
area: 282.843 m²,   error: -0.0287086 m²,   normalised error: 0.010149%
area: 283.164 m²,   error: 0.292786 m²,     normalised error: 0.103505%
area: 282.843 m²,   error: -0.0287018 m²,   normalised error: 0.0101466%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:1e-07 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 281.749 m²,   error: -1.12206 m²,     normalised error: 0.396666%
area: 282.843 m²,   error: -0.0283973 m²,   normalised error: 0.010039%
area: 281.749 m²,   error: -1.12206 m²,     normalised error: 0.396666%
area: 282.843 m²,   error: -0.0284534 m²,   normalised error: 0.0100588%
----- void do_test(int, F, F) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:10 offset:1e-07 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 10 m
area: 283.569 m²,   error: 0.697826 m²,     normalised error: 0.246694%
area: 282.843 m²,   error: -0.0284078 m²,   normalised error: 0.0100427%
area: 283.568 m²,   error: 0.696529 m²,     normalised error: 0.246235%
area: 282.843 m²,   error: -0.0283946 m²,   normalised error: 0.010038%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1000 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 1000 m
area: 2.82843e+06 m²,   error: -284.28 m²,  normalised error: 0.0100498%
area: 2.82843e+06 m²,   error: -284.27 m²,  normalised error: 0.0100494%
area: 2.82843e+06 m²,   error: -284.259 m², normalised error: 0.0100491%
area: 2.82843e+06 m²,   error: -284.27 m²,  normalised error: 0.0100494%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1000 offset:1e-07 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 1000 m
area: 2.82843e+06 m²,   error: -284.372 m², normalised error: 0.010053%
area: 2.82843e+06 m²,   error: -284.27 m²,  normalised error: 0.0100494%
area: 2.82843e+06 m²,   error: -284.282 m², normalised error: 0.0100499%
area: 2.82843e+06 m²,   error: -284.27 m²,  normalised error: 0.0100494%
----- void do_test(int, F, F) [with F = double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 1 m
area: 2.81749 m²,   error: -0.0112205 m²,   normalised error: 0.396663%
area: 2.8285 m²,    error: -0.000219998 m², normalised error: 0.0077773%
area: 2.83391 m²,   error: 0.0051987 m²,    normalised error: 0.183783%
area: 2.82848 m²,   error: -0.000234082 m², normalised error: 0.00827521%
----- void do_test(int, F, F) [with F = long double; DegreeOrRadian = boost::geometry::degree] n:8 radius:1 offset:0 ----
reference: (4.9, 52.1)
point: (4.9, 52.1)
radius: 1 m
area: 2.83535 m²,   error: 0.00663946 m²,   normalised error: 0.234717%
area: 2.82844 m²,   error: -0.000278463 m², normalised error: 0.00984417%
area: 2.83392 m²,   error: 0.005205 m²,     normalised error: 0.184006%
area: 2.82842 m²,   error: -0.000294424 m², normalised error: 0.0104084%

Some comments:

the use of different strategies (i.e. algorithms that perform the geographic computations in boost geometry) controls the accuracy but also the performance of algorithms.
there is still an issue in geographic buffer, feel free to file an issue on github to keep it in track
the "theoreticalArea" is accurate only for small areas, as the area grows the boost geometry algorithms are expected to be more accurate than that area.
the Earth is not flat ;)

@sehe just note that `bg::convert` does not convert between units so the example you have with radians is not working correctly. Instead you may use `bg::transform`. See https://github.com/boostorg/geometry/issues/818#issuecomment-785090121 — Vissarion, Feb 24 '21 at 16:16

Boost point_circle coming out in weird shapes

Edit 1

Edit 2

2 Answers2