c++boostgisboost-geometry
Boost point_circle coming out in weird shapes
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;
Solution
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 ;)
- C typedef name conflict
- How to use 32-bit pointers in 64-bit application?
- problems with scanf input stream
- Is it guaranteed to be safe to perform memcpy(0,0,0)?
- Need help understanding pointer to a type of n elements in C (int (*p)[10])?
- Local variables performance vs global in example
- How to Print SOAP Message Request Body generated using GSOAP when calling service
- Big Number Subtraction in C
- Docker custom program stop
- GTK4 Windows Multiple File Choosing
- Cannot get lldb to read file input through redirect
- Modeling nested lists in C
- Linker and Memory map for shared memory in autosar embedded environment
- Why is only the numerator cast to float to obtain a float quotient while dividing two integers?
- $gcc problemMatcher doesn't show all errors or warnings
- What registers to save in the ARM C calling convention?
- What does ** do in C language?
- How to install pywin32 module in windows 7
- How to differentiate a Cypher clause from an SQL clause in C?
- How do I check if my executable is running for a service or for a normal program?
- How to fix a simple error in Atmel studio
- Given day, month, year calculate day in year
- Definition of int64_t
- Can we have multiple c files to define functions for one header file?
- Check if array is ASCII
- What does C or GCC do when implicit function declaration exists
- How can I programmatically find the CPU frequency with C
- Can I use a binary literal in C or C++?
- Want to Convert Integer to String without itoa function
- How to permanently configure VSCode to use cl.exe without launching it from the Developer Command Prompt?