Geographic Area returns wrong results

[barendgehrels]

Hi @vissarion , @awulkiew ,
I'm busy producing buffers in geographical coordinate system around a point. That is work in progress now, and it starts to work. However, when I want to check the area (of the geographic buffer) I get unexpected results.

Below the debug output, containing:
number of input points, "geographic", number of output points (one more), area-in-geographic, "projected", number of output points (still one more as input), area-in-cartesian. Below that their WKT representations in both geographic and projected (using sterea for the Dutch coordinate system, but that is not relevant for this issue)

Because I specified one meter, one would expect an area of pi. Because the number of points is low, it is less, but in cartesian (I used Boost.Geometry projection/transformation), that all seems allright.

So I have to conclude that the area calculated by boost::geometry::area for a geographic coordinate system is not right... Or maybe I had to specify another strategy, but I debugged the code, and basically it seems to take the right strategies.

You see the geographic area fluctiating from -7106 via --51940 to +13597 etc. In only one case (not below, because from a previous debugging version) I got a value close to a right value...

5  geographic  6   -7106.98  projected  6   2.3767
     MULTIPOLYGON(((4.900000014692011 52.40000898675671,4.900013977465128 52.40000276851281,4.900008623854553 52.39999272427848,4.899991352373283 52.39999273484305,4.899986031615026 52.40000278560664,4.900000014692011 52.40000898675671)))
     MULTIPOLYGON(((121816.3935228445 490239.7834838452,121817.3390161667 490239.0853036283,121816.9671865939 490237.9703284562,121815.7918899588 490237.9794161207,121815.437346262 490239.1000077746,121816.3935228445 490239.7834838452)))
7  geographic  8   -51940.44953270301  projected  8   2.735323085114942
     MULTIPOLYGON(((4.900000014692011 52.40000898675671,4.900011495831415 52.40000559612448,4.90001432037284 52.39999799149619,4.900006361382522 52.39999189931258,4.899993612143401 52.39999190711099,4.899985673088598 52.39999800901907,4.899988522489213 52.40000561017674,4.900000014692011 52.40000898675671)))
     MULTIPOLYGON(((121816.3935228445 490239.7834838452,121817.1722615819 490239.40101414,121817.3587729585 490238.5537002923,121816.8126101027 490237.8795869154,121815.945044786 490237.8862951292,121815.4093713848 490238.5687735133,121815.6089628978 490239.4131019245,121816.3935228445 490239.7834838452)))
9  geographic  10   -28607.43220085142  projected  10   2.891395106998971
     MULTIPOLYGON(((4.900000014692011 52.40000898675671,4.900009455092996 52.4000068784781,4.900014471349158 52.4000015516828,4.900012716299604 52.39999549883819,4.900005011153675 52.39999155213724,4.89999496123439 52.39999155828457,4.899987269008388 52.39999551440371,4.899985533753316 52.40000156938328,4.899990567416462 52.40000689003125,4.900000014692011 52.40000898675671)))
     MULTIPOLYGON(((121816.3935228445 490239.7834838452,121817.0343521173 490239.5446116999,121817.3717130748 490238.9497042377,121817.2477507743 490238.2771252795,121816.720468555 490237.841581988,121816.0365876275 490237.8468699147,121815.5161034769 490238.2905147722,121815.4025564223 490238.9649302114,121815.7490763939 490239.554549747,121816.3935228445 490239.7834838452)))
11  geographic  12   144878.6498449351  projected  12   2.97234318849587
     MULTIPOLYGON(((4.900000014692011 52.40000898675671,4.90000795545699 52.40000755528196,4.900013370419436 52.40000372505826,4.900014540367223 52.39999871215498,4.900011093851504 52.39999410813363,4.900004125117584 52.39999137473812,4.899995846688667 52.39999137980185,4.89998888690606 52.3999941217171,4.899985455450992 52.39999872994556,4.899986641787119 52.40000374140755,4.899992069262416 52.40000756499918,4.900000014692011 52.40000898675671)))
     MULTIPOLYGON(((121816.3935228445 490239.7834838452,121816.9328105242 490239.6205935691,121817.2984238039 490239.1919980304,121817.3742830501 490238.6337732781,121817.1363034876 490238.1231517321,121816.6600419474 490237.8222521214,121816.0967081027 490237.8266079476,121815.6251564683 490238.1348362649,121815.3951013562 490238.6490767669,121815.4795836401 490239.2060617174,121815.8517807918 490239.6289523419,121816.3935228445 490239.7834838452)))
13  geographic  14   13597.64092986423  projected  14   3.019500571803292
     MULTIPOLYGON(((4.900000014692011 52.40000898675671,4.900006840724092 52.40000795320133,4.900012099627947 52.4000050976631,4.900014586651976 52.40000107431168,4.900013732049592 52.39999680484858,4.900009731600762 52.39999326735627,4.900003501760424 52.39999127223135,4.899996469709414 52.3999912765327,4.899990246404989 52.39999327927492,4.899986257530689 52.39999682165411,4.899985416889869 52.40000109215416,4.89998791706407 52.40000511245501,4.899993185294162 52.40000796155404,4.900000014692011 52.40000898675671)))
     MULTIPOLYGON(((121816.3935228445 490239.7834838452,121816.8572522716 490239.6653730333,121817.212975916 490239.3452833782,121817.3792017778 490238.8965435647,121817.3178495146 490238.4219544754,121817.042974191 490238.0302387534,121816.6175464315 490237.8111337412,121816.1390266085 490237.8148337687,121815.7170378461 490238.0404911955,121815.4482526841 490238.4364106264,121815.3942465646 490238.911891695,121815.5673916392 490239.3580074163,121815.9280224574 490239.6725580526,121816.3935228445 490239.7834838452)))
15  geographic  16   34802.96352132662  projected  16   3.049312928096697
     MULTIPOLYGON(((4.900000014692011 52.40000898675671,4.900005989198468 52.40000820615538,4.900010928117561 52.40000600663501,4.900013977465128 52.40000276851281,4.900014609981563 52.3999990516897,4.900012716299604 52.39999549883819,4.900008623854553 52.39999272427848,4.900003040267405 52.399991207757,4.899996930990699 52.39999121149386,4.899991352373283 52.39999273484303,4.899987269008388 52.39999551440371,4.899985386946967 52.39999906956475,4.899986031615026 52.40000278560664,4.899989091544179 52.40000601999194,4.899994037645167 52.40000821346587,4.900000014692011 52.40000898675671)))
     MULTIPOLYGON(((121816.3935228445 490239.7834838452,121816.7994963964 490239.693904431,121817.133936763 490239.4469437796,121817.3390161667 490239.0853036283,121817.3792745153 490238.6715148503,121817.2477507743 490238.2771252795,121816.9671865939 490237.9703284562,121816.5860940701 490237.8041723452,121816.1703675617 490237.80738683,121815.7918899588 490237.9794161183,121815.5161034769 490238.2905147722,121815.3906941045 490238.6868910086,121815.437346262 490239.1000077746,121815.6479933675 490239.4584334347,121815.9862126222 490239.7001929283,121816.3935228445 490239.7834838452)))

[vissarion]

Hi @barendgehrels I looked at the issue. It seems that the problem does not come from the area algorithm but from the initialization of the geometry. How are you initializing your geometries?

For your first example polygon in your list, I tried the following cases of coordinate types (used to parameterize everything, i.e. points, strategies, spheroid) and initialization methods:

coordinate type	initialization method	area result
double	read_wkt	-871.955
double	c++11 brace initialization	-871.955
long double	read_wkt	2.06332
long double	c++11 brace initialization	2.20349

In all cases, the resulting polygon was different and thus the area. Note that since the input coordinates are in degrees and the result in m^2 a small perturbation in the input (that may come from accuracy loss at initialization) will result a large difference in the output. This explains the large variance of your results.

[barendgehrels]

Hi @vissarion , these geometries are calculated by my concept buffer-around-point-in-geo implementation. Using the direct formula as used in the other ticket. So basically I think the initialization is perfect, above all because if I project it (as I did) and then calculate the area, I get the correct results. So the geometry itself should be fine...

In your table above it looks to me that just the coordinate type is the main difference. And the initialization results just in a 'small' difference.

I use doubles. And that should work (even float should work IMO...). The area results of around 2 that you get for long doubles look OK to me (though I don't understand the difference, and how a minor perturbation can lead to 10% of a difference in results - I don't get that yet).

I will do some comparisons with double/long double later, and also check the perimeter.

[vissarion]

You are right, I verify that the geometries are fine. It seems that the inaccuracy comes from summing up values that are very close to each other here : https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/strategies/geographic/area.hpp#L186
That is why computation with long double overcomes the inaccuracy. On the other hand, Kahan summation didn't seem to help. Needs further investigation.

[barendgehrels]

Hi @vissarion thanks. I now tried with long double and it is better - but also not perfect yet...

Below my new output. As you can see the area is more in the range. But it is still unusuable...

I now also added perimeter and that looks fine, for long double. And also for double it is reasonable (but not as good as their projected values).

Output for long double:

numpoints system numpoints+1 area perimeter system numpoints+1 area perimeter
5  geographic  6   -8.70934   5.87787  projected  6   2.3767   5.87668
6  geographic  7   -9.605667940099227   6.000014405269862  projected  7   2.597044063790236   5.998808060343662
7  geographic  8   -28.6840258191935   6.074389874202544  projected  8   2.735323084663833   6.073165633442998
8  geographic  9   11.04197348580413   6.122952157160927  projected  9   2.827303468897298   6.121718560789472
9  geographic  10   -7.609848928965851   6.156385383740395  projected  10   2.891395110484154   6.155139576603446
10  geographic  11   -17.09551017365113   6.180361705603756  projected  11   2.937758704316366   6.17911212344289
11  geographic  12   -4.66994831909932   6.198135973761448  projected  12   2.972343191977416   6.196884953050459
12  geographic  13   31.85406111249817   6.211662920149589  projected  13   2.998808179052503   6.210423096128881
13  geographic  14   -48.36471198818704   6.222224786814641  projected  14   3.019500574460835   6.220971190297066
14  geographic  15   54.35439501514531   6.230600659644504  projected  15   3.035979579948616   6.229348404247994
15  geographic  16   -20.65227759466741   6.237356143066346  projected  16   3.049312928953441   6.236111632974661
16  geographic  17   -13.85074591323568   6.242891438249192  projected  17   3.060251218981648   6.241650113948441
17  geographic  18   -73.97049289480156   6.24752317305389  projected  18   3.069334310953309   6.246242500893913
18  geographic  19   -28.63149905632986   6.25131874333177  projected  19   3.07695840705037   6.250092524928219
19  geographic  20   34.77809698644213   6.254614257579379  projected  20   3.083419511918692   6.253351917263175
20  geographic  21   56.66763616728842   6.257391110435357  projected  21   3.088942300268172   6.256135532030151
21  geographic  22   -142.6508088952368   6.259828080230406  projected  22   3.093699797324916   6.258531634024919
22  geographic  23   41.29627763257787   6.261880830127452  projected  23   3.097826950365743   6.260608925253983
23  geographic  24   44.21540304844772   6.263694919058841  projected  24   3.101430263120164   6.262421541250442
24  geographic  25   148.2382582753415   6.265246057958701  projected  25   3.104594677340629   6.264012592198701
25  geographic  26   -81.31339971108426   6.266732608215863  projected  26   3.107388616071603   6.26541676566701
26  geographic  27   -46.7369579968207   6.267942912965663  projected  27   3.109867678206228   6.266662214475424
27  geographic  28   -15.21189752858295   6.269056176908831  projected  28   3.112077415401018   6.26777197954567
28  geographic  29   -6.397964923270046   6.26998728525905  projected  29   3.114055449885655   6.26876508146113
29  geographic  30   -130.1152857882636   6.270888134743907  projected  30   3.115833052288053   6.269657312203098
30  geographic  31   -11.47674017764075   6.271697263318989  projected  31   3.117436394232755   6.270461879818648
31  geographic  32   6.880351493295484   6.272395970870092  projected  32   3.118887515793631   6.271189901466996
32  geographic  33   101.1286344815155   6.273091237410482  projected  33   3.12020508393789   6.271850788937381
33  geographic  34   50.09235745211329   6.273748100630176  projected  34   3.121404985066874   6.272452544650028
34  geographic  35   -124.6317943026502   6.274227133920037  projected  35   3.122500820942832   6.273002021307129
35  geographic  36   -43.07674325804758   6.274725652958071  projected  36   3.123504277623624   6.273505100582522
36  geographic  37   77.96110927570538   6.275227769598489  projected  37   3.124425451146408   6.273966863875912
37  geographic  38   93.88006571831122   6.275714758140596  projected  38   3.125273093441592   6.274391714174344
38  geographic  39   68.48816641905242   6.276038073775148  projected  39   3.126054823096638   6.274783481943314
39  geographic  40   -39.92998385169993   6.276431653077836  projected  40   3.126777296911314   6.275145513549985
40  geographic  41   113.7640605137189   6.276636867917465  projected  41   3.127446355162419   6.27548074536756
41  geographic  42   71.93103433272893   6.277014318975496  projected  42   3.128067133469131   6.275791757293996
42  geographic  43   110.9793271047548   6.277307831640422  projected  43   3.128644166934009   6.276080828097051
43  geographic  44   117.7244887709063   6.277602084472294  projected  44   3.129181469495279   6.276349973633893
44  geographic  45   293.2099410003796   6.277809590166146  projected  45   3.12968260750813   6.276600985915282
45  geographic  46   -146.8255608281281   6.278018187354725  projected  46   3.130150750192797   6.276835454094177
46  geographic  47   -55.85972673204054   6.278335461919328  projected  47   3.130588734272351   6.277054803584052
47  geographic  48   53.87766418303545   6.278610416958765  projected  48   3.13099909205414   6.277260304023793
48  geographic  49   144.2215191312397   6.278661087707991  projected  49   3.131384101293179   6.277453099946584
49  geographic  50   -345.4872121693842   6.27898300414949  projected  50   3.13174580628413   6.277634216836573
50  geographic  51   -43.76808189643426   6.279022799439124  projected  51   3.132086049988629   6.277804578903483

Output for double

5  geographic  6   -7106.98   5.87145  projected  6   2.3767   5.87668
6  geographic  7   9877.873026939991   5.990441149585183  projected  7   2.59704417915782   5.998808193597772
7  geographic  8   -51940.44953270301   6.076608045950909  projected  8   2.735323085114942   6.073165633977603
8  geographic  9   -15801.34382348935   6.123263339124103  projected  9   2.827303563128226   6.121718662777818
9  geographic  10   -28607.43220085142   6.153498692251656  projected  10   2.891395106998971   6.155139572854291
10  geographic  11   4168.449447946321   6.179248947943594  projected  11   2.937758780579315   6.179112203671658
11  geographic  12   144878.6498449351   6.157551484274717  projected  12   2.97234318849587   6.196884949366017
12  geographic  13   -28710.3576541802   6.243436350713854  projected  13   2.998808245720284   6.210423165152673
13  geographic  14   13597.64092986423   6.248798925754039  projected  14   3.019500571803292   6.220971187512639
14  geographic  15   -86283.87965229033   6.2561429224446  projected  15   3.035979636755656   6.229348462560575
15  geographic  16   34802.96352132662   6.206895958733149  projected  16   3.049312928096697   6.236111632198725
16  geographic  17   49095.97289276295   6.238712639300454  projected  17   3.060251271916059   6.241650167907099
17  geographic  18   73707.87018695462   6.223765623920101  projected  18   3.069334308874204   6.246242498768325
18  geographic  19   -17501.06691975416   6.288031507884794  projected  19   3.076958450678831   6.250092569270165
19  geographic  20   11355.60807726318   6.260196866396992  projected  20   3.08341950580143   6.25335191104325
20  geographic  21   54750.68440883196   6.174727768988833  projected  21   3.088942338168636   6.256135570454799
21  geographic  22   -91423.81211860616   6.28608924130507  projected  22   3.093699793505039   6.258531630092383
22  geographic  23   -40560.27572535454   6.208167131811362  projected  23   3.097826987302142   6.260608962512189
23  geographic  24   -149152.3787578737   6.134579246128111  projected  24   3.101430254951083   6.262421533036529
24  geographic  25   189158.8023008354   6.315632739554544  projected  25   3.104594711631762   6.264012626793592
25  geographic  26   253645.7206160811   6.17313195148979  projected  26   3.107388613682815   6.265416763258954
26  geographic  27   19394.72616995342   6.314308947007198  projected  27   3.109867707887588   6.266662244428333
27  geographic  28   25803.50541588549   6.226280367672984  projected  28   3.112077415569274   6.267771979632485
28  geographic  29   -289165.7023337696   6.42299700821913  projected  29   3.114055478587034   6.268765110350144
29  geographic  30   138687.6766455288   6.277061288152301  projected  30   3.115833049064349   6.269657308967967
30  geographic  31   133580.832379287   6.382206571971186  projected  31   3.117436420543072   6.270461906161152
31  geographic  32   -154021.4551534558   6.119565526811111  projected  32   3.11888751231163   6.271189898033787
32  geographic  33   65016.62865333557   6.2142592428263  projected  33   3.120205110878032   6.271850816016245
33  geographic  34   700411.0785136138   6.213222877295266  projected  34   3.121404981123305   6.272452540676222
34  geographic  35   -41061.30295670583   6.178680557974025  projected  35   3.122500844603564   6.273002045068829
35  geographic  36   -549351.9575878523   6.374040143397463  projected  36   3.123504272464515   6.273505095333792
36  geographic  37   -22154.422378672   6.17649265244964  projected  37   3.124425472392659   6.27396688534256
37  geographic  38   15854.78040639361   6.472128723775513  projected  38   3.125273089965503   6.27439171068878
38  geographic  39   171666.0336067806   6.283664699602097  projected  39   3.126054841625546   6.27478350054973
39  geographic  40   -471754.6958557257   6.073149205169562  projected  40   3.126777293243322   6.275145509904249
40  geographic  41   -98102.50219920019   5.799508218696191  projected  41   3.12744637304354   6.275480763371528
41  geographic  42   -296735.0928373209   5.934658108771179  projected  42   3.128067131410489   6.275791755199044
42  geographic  43   -409594.8254779964   6.140941248298602  projected  43   3.128644183363349   6.276080844577871
43  geographic  44   84753.42718479398   5.933002239498226  projected  44   3.129181464903468   6.276349969081677
44  geographic  45   284775.5204423611   5.902864570947456  projected  45   3.129682623452709   6.276601002008998
45  geographic  46   -82284.61634465648   5.970360164674502  projected  46   3.130150745918627   6.276835449649665
46  geographic  47   -702019.7854401803   6.186378752826677  projected  47   3.130588747509023   6.277054816761917
47  geographic  48   -933381.192933112   6.216717190584058  projected  48   3.130999090044384   6.277260301957213
48  geographic  49   372651.5107663177   6.202817750333533  projected  49   3.131384117016069   6.277453115733204
49  geographic  50   -642122.396644773   6.500462042484167  projected  50   3.131745803445483   6.277634213933462
50  geographic  51   -1615660.075003426   6.219329432941888  projected  51   3.132086063875477   6.277804592790042

[barendgehrels]

And indeed, if I use a buffer distance of 100 meter instead of 1 meter, then coordinates are further apart, and I get a neat output, and also much better for double

Output for long double

5  geographic  6   24406.7   587.787  projected  6   23767   587.668
10  geographic  11   29389.3917650142   618.0358232388136  projected  11   29377.58735880721   617.9112156599526
15  geographic  16   30505.42137202255   623.7369244367449  projected  16   30493.12967984087   623.6111672907244
20  geographic  21   30902.35272895053   625.7397181644264  projected  21   30889.42341354454   625.6135573623785
25  geographic  26   31088.36394180376   626.6680287578648  projected  26   31073.8865827361   626.5416808175752
30  geographic  31   31185.06358378872   627.1726429540689  projected  31   31174.36438939671   627.0461924742734
35  geographic  36   31248.2589853291   627.4770254838768  projected  36   31235.04324843551   627.3505147949721
40  geographic  41   31287.61311179867   627.6746299245902  projected  41   31274.46400490234   627.548079086611
45  geographic  46   31312.62648794206   627.8101281412048  projected  46   31301.50797059848   627.683550099199
50  geographic  51   31335.88638699717   627.9070597607381  projected  51   31320.86095820805   627.7804624773828

Output for double

5  geographic  6   24424.5   587.787  projected  6   23767   587.668
10  geographic  11   28954.61113643066   618.0358844756693  projected  11   29377.58735513408   617.9112156213394
15  geographic  16   30467.21497418519   623.7366198510647  projected  16   30493.12967984623   623.6111672908565
20  geographic  21   29999.62650359623   625.7399310790181  projected  21   30889.4234114656   625.6135573413259
25  geographic  26   32533.4898053848   626.6691447541973  projected  26   31073.88658259891   626.5416808161938
30  geographic  31   28521.01030769165   627.1729623150326  projected  31   31174.36438796759   627.0461924598916
35  geographic  36   26814.83778598347   627.4765069980294  projected  36   31235.0432483115   627.350514793699
40  geographic  41   26692.62545923042   627.6760558573422  projected  41   31274.46400359846   627.5480790734821
45  geographic  46   23815.63281618918   627.8084416608789  projected  46   31301.50797025645   627.683550095775
50  geographic  51   33650.51412680989   627.90694786202  projected  51   31320.86095703265   627.7804624655969

[awulkiew]

@barendgehrels Out of curiosity, you're using andoyer (so for direct this is thomas without the 2nd degree approximation), correct? Have you tried different FormulaPolicies too? In case we didn't know where is the issue this could help to localize it. If the accuracy was low for all geodesic formulas we'd know that the error is in the area.

@vissarion I'm also curious if maybe the SeriesOrder in area is to blame, i.e. are we sure that the order in area should be the same as in geodesic formulas? Area is squared length after all. I don't know what's the formula in area used for (I haven't looked so this is a wild guess) but if it takes some number of distances and returns the area is it possible that we need greater SeriesOrder to get comparable accuracy?

[barendgehrels]

@awulkiew yes, for buffer-generation I'm using the Andoyer indeed, by specifying typename FormulaPolicy = strategy::andoyer and I could change that. But the generation on itself seems OK now, as long as I don't touch the 180 degrees.

For area I also use the default, by calling: bg::area(buffer) - so maybe I should specify the strategy to get better results.

[vissarion]

@awulkiew SeriesOrder does not change the behaviour. Since the problem is on short segments (where the ellipsoidal term is small anyway) we can just avoid it. To be more clear the algorithm adds a spherical and an ellipsoidal term so I propose to avoid the latter for short segments. @barendgehrels how are your examples behaviour with spherical CS?

[barendgehrels]

@vissarion : they behave well! Here the new table, using double, with the spherical area right after the geographic area, calculated by bg::area(buffer, bg::strategy::area::spherical<>(6371228.0))

5  geographic  6   -7106.98   2.36845   5.87145  projected  6   2.3767   5.87668
10  geographic  11   4168.449447946321   2.927570082965481   6.179248947943594  projected  11   2.937758780579315   6.179112203671658
15  geographic  16   34802.96352132662   3.038737335821288   6.206895958733149  projected  16   3.049312928096697   6.236111632198725
20  geographic  21   54750.68440883196   3.078229307965614   6.174727768988833  projected  21   3.088942338168636   6.256135570454799
25  geographic  26   253645.7206160811   3.096611607288448   6.17313195148979  projected  26   3.107388613682815   6.265416763258954
30  geographic  31   133580.832379287   3.106624562744787   6.382206571971186  projected  31   3.117436420543072   6.270461906161152
35  geographic  36   -549351.9575878523   3.112671374623749   6.374040143397463  projected  36   3.123504272464515   6.273505095333792
40  geographic  41   -98102.50219920019   3.116599802490512   5.799508218696191  projected  41   3.12744637304354   6.275480763371528
45  geographic  46   -82284.61634465648   3.119294795856522   5.970360164674502  projected  46   3.130150745918627   6.276835449649665
50  geographic  51   -1615660.075003426   3.121223396694027   6.219329432941888  projected  51   3.132086063875477   6.277804592790042

[vissarion]

@barendgehrels can we find the threshold at which the geographic formula (as is) is inaccurate? As you tested for 1m is inaccurate but for 100m is fine.

[barendgehrels]

Indeed, I can probably get such a threshold, in meters, for you. I try to get it tonight.

[barendgehrels]

@vissarion : these are my found values for the radius where area calculations are stable:
Float: 19080
Double: 605
Long double: 45

So actually double was not yet fine at 100m, because the output areas were reasonable but not yet stable. My test increases the number of points in a circle, the area therefore should increase (until it is pi, for a radius of 1). If it decreases it is not stable.

So therefore the limit is higher. For float it is quite high (half the earth?)... But most probably people will use at least double for geographic cs (?).

I tested at 3 locations (one in Amsterdam, one at ~22 N, one at ~82N - but there is no much variation, I took the highest values of these three.

So the found radius is the radius of a buffer (around a point) with 5 up to 50 points (increasingly). So that might be an appropriate initial measurement for the threshold. It is not the distance from point to point. But for 5 or 6 points it has the same order, and that is where I started.

[vissarion]

@barendgehrels Thanks for the feedback. Still I cannot understand those numbers. Maybe there are several issues but I tried to focus a bit on the affect of inverse formulas (i.e. andoyer, thomas, vincenty). They are used to calculate the azimuth at the beginning of area computation. Surprisingly, for the data I tried (your data for 1m buffer) thomas resulted to the higher accuracy while vincenty to the lowest. Could you try to your whole dataset so we have a better idea, or tell me how to reproduce your data and I will do the test. To be clear I tried the following

bg::strategy::area::geographic
<
  inverse_strategy,
  5,
  bg::srs::spheroid<CT>,
  CT
> str;
bg::area(polygon, str);

Where inverse_strategy can be either bg::strategy::andoyer, bg::strategy::thomas, or bg::strategy::vincenty, CT is the number type, and 5 is the order of approximation (it does not affect the results in my case).