geometry
geometry copied to clipboard
boostorg
Geographic Area returns wrong results
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)))
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.
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.
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.
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
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
@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?
@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.
@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?
@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
@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.
Indeed, I can probably get such a threshold, in meters, for you. I try to get it tonight.
@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.
@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).
Right, I see. I did not use it exactly like this. I started with bg::area(polygon) which obviously goes wrong (for my case). I will post a complete minimal scenario, one of these days. If necessary - I will substitute what you did first to verify my findings.
Hi @vissarion - below is my new table, followed by the test program. In the table: first column is radius, then all kinds of areas with ratio to projected area (for now considered as the truth). First for float (f), then double (d) and then long double (e).
Thomas performs the best. I didn't know that because I didn't try it before. For long double, it is perfect even for one meter! Float is usually bad - most probably we should promote that internally and it will be fixed. Based on the table below, if we promote all to long double and use thomas - it will work correctly...
So (based on your hint, and I also changed my test program a bit) my findings are different now and mainly depend on used inverse_strategy.
testing for f
1 geographic -2616650.022 (-55248.998) andoyer -2616650.015 (-55248.997) thomas -2616884.666 (-55253.952) vincenty -2616880.404 (-55253.862) spherical 47.162 (0.996) projected 47.361
2 geographic 2044950.046 (18246.779) andoyer 2044950.041 (18246.779) thomas -8326487.667 (-74295.984) vincenty -8326396.453 (-74295.171) spherical 111.666 (0.996) projected 112.072
4 geographic 1741861.440 (7072.469) andoyer 1741861.429 (7072.469) thomas 6661512.407 (27047.696) vincenty 6596567.215 (26784.000) spherical 245.457 (0.997) projected 246.288
8 geographic -8468850.767 (-14662.330) andoyer -8468850.757 (-14662.330) thomas 578.399 (1.001) vincenty -3908.772 (-6.767) spherical 575.644 (0.997) projected 577.592
16 geographic 20164620.681 (13140.102) andoyer 20164620.604 (13140.102) thomas 1533.261 (0.999) vincenty 18778.747 (12.237) spherical 1529.411 (0.997) projected 1534.586
32 geographic 150640065.989 (32514.002) andoyer 150640065.466 (32514.002) thomas 4638.348 (1.001) vincenty 13597.633 (2.935) spherical 4617.457 (0.997) projected 4633.083
64 geographic -306584731.272 (-19904.868) andoyer -306584730.490 (-19904.868) thomas 15417.078 (1.001) vincenty 15913.975 (1.033) spherical 15350.553 (0.997) projected 15402.500
128 geographic -694773483.095 (-12455.078) andoyer -694773480.986 (-12455.078) thomas 55817.069 (1.001) vincenty 56655.112 (1.016) spherical 55594.213 (0.997) projected 55782.347
256 geographic 190560276.129 (906.418) andoyer 190560275.143 (906.418) thomas 210317.872 (1.000) vincenty 210709.415 (1.002) spherical 209525.477 (0.997) projected 210234.526
512 geographic 1178331674.125 (1444.692) andoyer 1178331670.301 (1444.692) thomas 815989.978 (1.000) vincenty 816433.473 (1.001) spherical 812877.278 (0.997) projected 815628.126
1024 geographic -527570208.535 (-157.729) andoyer -527570208.727 (-157.729) thomas 3346138.595 (1.000) vincenty 3346159.365 (1.000) spherical 3333501.613 (0.997) projected 3344782.350
2048 geographic 8959153318.996 (674.797) andoyer 8959153288.916 (674.797) thomas 13282138.011 (1.000) vincenty 13282195.756 (1.000) spherical 13232038.375 (0.997) projected 13276816.420
4096 geographic -15045873900.053 (-287.173) andoyer -15045873851.174 (-287.173) thomas 52414219.343 (1.000) vincenty 52414204.389 (1.000) spherical 52216378.481 (0.997) projected 52393088.902
8192 geographic -4675164620.233 (-22.239) andoyer -4675164605.840 (-22.239) thomas 210303212.277 (1.000) vincenty 210303227.115 (1.000) spherical 209509936.500 (0.997) projected 210219019.846
16384 geographic 841593524.763 (1.002) andoyer 841593523.057 (1.002) thomas 840475852.740 (1.000) vincenty 840475842.658 (1.000) spherical 837302908.768 (0.997) projected 840137824.184
testing for d
1 geographic -590757.726 (-189078.175) andoyer -590757.725 (-189078.174) thomas 0.000 (0.000) vincenty -0.927 (-0.297) spherical 3.114 (0.997) projected 3.124
2 geographic -75092.034 (-6008.497) andoyer -75092.034 (-6008.497) thomas 12.476 (0.998) vincenty -2.411 (-0.193) spherical 12.455 (0.997) projected 12.498
4 geographic -44830.186 (-896.773) andoyer -44830.186 (-896.773) thomas 50.041 (1.001) vincenty -9.739 (-0.195) spherical 49.822 (0.997) projected 49.991
8 geographic -45412.033 (-227.103) andoyer -45412.033 (-227.103) thomas 199.987 (1.000) vincenty -38.992 (-0.195) spherical 199.288 (0.997) projected 199.962
16 geographic -2415.338 (-3.020) andoyer -2415.338 (-3.020) thomas 800.173 (1.000) vincenty -155.643 (-0.195) spherical 797.151 (0.997) projected 799.849
32 geographic 12972.208 (4.055) andoyer 12972.208 (4.055) thomas 3200.669 (1.000) vincenty 13070.195 (4.085) spherical 3188.605 (0.997) projected 3199.396
64 geographic 10564.452 (0.826) andoyer 10564.452 (0.826) thomas 12802.682 (1.000) vincenty 21033.070 (1.644) spherical 12754.422 (0.997) projected 12797.584
128 geographic 48182.060 (0.941) andoyer 48182.060 (0.941) thomas 51210.939 (1.000) vincenty 50955.079 (0.995) spherical 51017.687 (0.997) projected 51190.336
256 geographic 204964.881 (1.001) andoyer 204964.881 (1.001) thomas 204843.977 (1.000) vincenty 203901.637 (0.996) spherical 204070.754 (0.997) projected 204761.346
512 geographic 819471.930 (1.001) andoyer 819471.930 (1.001) thomas 819375.949 (1.000) vincenty 819152.330 (1.000) spherical 816283.046 (0.997) projected 819045.416
1024 geographic 3277740.667 (1.000) andoyer 3277740.670 (1.000) thomas 3277504.330 (1.000) vincenty 3277381.984 (1.000) spherical 3265132.427 (0.997) projected 3276181.924
2048 geographic 13109957.429 (1.000) andoyer 13109957.440 (1.000) thomas 13110019.639 (1.000) vincenty 13110025.356 (1.000) spherical 13060531.615 (0.997) projected 13104729.856
4096 geographic 52440457.107 (1.000) andoyer 52440457.149 (1.000) thomas 52440093.556 (1.000) vincenty 52440113.002 (1.000) spherical 52242141.106 (0.997) projected 52418938.059
8192 geographic 209761676.873 (1.000) andoyer 209761677.040 (1.000) thomas 209760481.410 (1.000) vincenty 209760493.006 (1.000) spherical 208968671.714 (0.997) projected 209675922.902
16384 geographic 839047232.375 (1.000) andoyer 839047233.041 (1.000) thomas 839042616.289 (1.000) vincenty 839042622.334 (1.000) spherical 835875386.941 (0.997) projected 838705402.291
testing for e
1 geographic 189.969 (60.801) andoyer 189.969 (60.801) thomas 3.126 (1.000) vincenty -0.608 (-0.195) spherical 3.114 (0.997) projected 3.124
2 geographic 137.277 (10.984) andoyer 137.277 (10.984) thomas 12.503 (1.000) vincenty -16.970 (-1.358) spherical 12.455 (0.997) projected 12.498
4 geographic 38.234 (0.765) andoyer 38.234 (0.765) thomas 50.011 (1.000) vincenty 67.352 (1.347) spherical 49.822 (0.997) projected 49.991
8 geographic 209.906 (1.050) andoyer 209.906 (1.050) thomas 200.043 (1.000) vincenty 208.040 (1.040) spherical 199.288 (0.997) projected 199.962
16 geographic 796.485 (0.996) andoyer 796.485 (0.996) thomas 800.172 (1.000) vincenty 790.812 (0.989) spherical 797.151 (0.997) projected 799.849
32 geographic 3199.594 (1.000) andoyer 3199.594 (1.000) thomas 3200.688 (1.000) vincenty 3201.491 (1.001) spherical 3188.605 (0.997) projected 3199.396
64 geographic 12797.747 (1.000) andoyer 12797.747 (1.000) thomas 12802.750 (1.000) vincenty 12804.544 (1.001) spherical 12754.422 (0.997) projected 12797.584
128 geographic 51211.510 (1.000) andoyer 51211.510 (1.000) thomas 51211.000 (1.000) vincenty 51209.779 (1.000) spherical 51017.687 (0.997) projected 51190.336
256 geographic 204845.098 (1.000) andoyer 204845.098 (1.000) thomas 204844.006 (1.000) vincenty 204843.414 (1.000) spherical 204070.754 (0.997) projected 204761.346
512 geographic 819380.579 (1.000) andoyer 819380.580 (1.000) thomas 819376.053 (1.000) vincenty 819376.109 (1.000) spherical 816283.046 (0.997) projected 819045.416
1024 geographic 3277522.367 (1.000) andoyer 3277522.370 (1.000) thomas 3277504.458 (1.000) vincenty 3277504.375 (1.000) spherical 3265132.427 (0.997) projected 3276181.924
2048 geographic 13110091.522 (1.000) andoyer 13110091.532 (1.000) thomas 13110019.743 (1.000) vincenty 13110019.993 (1.000) spherical 13060531.615 (0.997) projected 13104729.856
4096 geographic 52440380.626 (1.000) andoyer 52440380.668 (1.000) thomas 52440093.629 (1.000) vincenty 52440094.598 (1.000) spherical 52242141.107 (0.997) projected 52418938.060
8192 geographic 209761629.539 (1.000) andoyer 209761629.705 (1.000) thomas 209760481.451 (1.000) vincenty 209760485.160 (1.000) spherical 208968671.714 (0.997) projected 209675922.902
16384 geographic 839047208.844 (1.000) andoyer 839047209.511 (1.000) thomas 839042616.438 (1.000) vincenty 839042619.881 (1.000) spherical 835875386.940 (0.997) projected 838705402.291
And this is the standalone test program:
#include <boost/geometry.hpp>
#include <boost/geometry/srs/projection.hpp>
#include <boost/geometry/srs/transformation.hpp>
#include <boost/geometry/strategies/transform/srs_transformer.hpp>
#include <iostream>
namespace boost { namespace geometry
{
template
<
typename FormulaPolicy = strategy::andoyer,
typename Spheroid = srs::spheroid<double>,
typename CalculationType = void
>
class point_circle_geographic
{
public :
explicit point_circle_geographic(std::size_t count = 90)
: m_count((count < 3u) ? 3u : count)
{}
template
<
typename Point,
typename OutputRange,
typename Distance
>
inline void apply(Point const& point,
Distance const& buffer_distance,
OutputRange& output_range) const
{
typedef typename boost::range_value<OutputRange>::type output_point_type;
typedef typename select_most_precise
<
typename geometry::coordinate_type<Point>::type,
typename geometry::coordinate_type<output_point_type>::type,
CalculationType
>::type calculation_type;
typedef typename FormulaPolicy::template direct
<
calculation_type, true, false, false, false
> direct_t;
calculation_type const two_pi = geometry::math::two_pi<calculation_type>();
calculation_type const diff = two_pi / calculation_type(m_count);
// TODO: after calculation of some angles is corrected, we can start at 0.0
calculation_type angle = 0.001;
for (std::size_t i = 0; i < m_count; i++, angle += diff)
{
typename direct_t::result_type
dir_r = direct_t::apply(get_as_radian<0>(point), get_as_radian<1>(point),
buffer_distance, angle,
m_spheroid);
output_point_type p;
set_from_radian<0>(p, dir_r.lon2);
set_from_radian<1>(p, dir_r.lat2);
output_range.push_back(p);
}
{
// Close it:
const output_point_type p = output_range.front();
output_range.push_back(p);
}
}
private :
std::size_t m_count;
Spheroid m_spheroid;
};
}}
namespace bg = boost::geometry;
template <typename PointType>
void test(double radius, int point_count)
{
typedef typename bg::select_most_precise
<
typename bg::coordinate_type<PointType>::type,
double
>::type ct;
typedef bg::model::polygon<PointType> polygon_type;
PointType center(4.9, 52.0);
polygon_type buffer;
// Generate point buffer
bg::point_circle_geographic<> generator(point_count);
generator.apply(center, radius, bg::exterior_ring(buffer));
// Make a projection (~equal area)
typedef bg::projections::parameters<ct> parameters_type;
typedef bg::projections::sterea_ellipsoid<ct, parameters_type> projection_type;
bg::srs::detail::proj4_parameters params( "+lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m");
parameters_type par = bg::projections::detail::pj_init<ct>(params);
bg::strategy::transform::srs_forward_transformer<projection_type> transformer(params, par);
typedef bg::model::d2::point_xy<ct> cartesian_point_type;
typedef bg::model::polygon<cartesian_point_type> cartesian_polygon_type;
cartesian_polygon_type cart_buffer;
bg::transform(buffer, cart_buffer, transformer);
// Calculate various areas
bg::strategy::area::geographic
<
bg::strategy::vincenty, 5, bg::srs::spheroid<ct>, ct
> vincenty_strategy;
bg::strategy::area::geographic
<
bg::strategy::andoyer, 5, bg::srs::spheroid<ct>, ct
> andoyer_strategy;
bg::strategy::area::geographic
<
bg::strategy::thomas, 5, bg::srs::spheroid<ct>, ct
> thomas_strategy;
const ct default_area = bg::area(buffer);
const ct spherical_area = bg::area(buffer, bg::strategy::area::spherical<>(6371228.0));
const ct cart_area = bg::area(cart_buffer);
const ct andoyer_area = bg::area(buffer, andoyer_strategy);
const ct thomas_area = bg::area(buffer, thomas_strategy);
const ct vincenty_area = bg::area(buffer, vincenty_strategy);
// Report
std::cout << int(radius)
<< std::fixed << std::setprecision(3)
<< " geographic " << default_area << " (" << default_area / cart_area << ")"
<< " andoyer " << andoyer_area << " (" << andoyer_area / cart_area << ")"
<< " thomas " << thomas_area << " (" << thomas_area / cart_area << ")"
<< " vincenty " << vincenty_area << " (" << vincenty_area / cart_area << ")"
<< " spherical " << spherical_area << " (" << spherical_area / cart_area << ")"
<< " projected " << cart_area
<< std::endl;
}
template <typename CoordinateType>
void test_all(int point_count)
{
std::cout << "testing for " << typeid(CoordinateType).name() << std::endl;
double radius = 1.0;
for (int i = 0; i < 15; i++, radius *= 2.0)
{
test<bg::model::point<CoordinateType, 2, bg::cs::geographic<bg::degree> > >(radius, point_count);
}
}
int main(int, char* [])
{
test_all<float>(36);
test_all<double>(36);
test_all<long double>(36);
return 0;
}
@barendgehrels thanks for the results and code. Below I added one more formula, named series below, implemented here: https://github.com/boostorg/geometry/pull/500 (with some modifications to use radians instead of degrees). This is expected to be the most accurate method when used with long double and "verifies" the correctness of thomas formula.
Two comments here: a) the results are expected to be different for longer distances (see https://github.com/boostorg/geometry/pull/431#issuecomment-367357573), b) vincenty appears to be less accurate than (the more efficient) thomas in this dataset which is at least to my understanding unexpected.
Both comments are written with distance computation in mind but here we compute azimuths so it needs further investigation.
testing for f
1 geographic -2616650.022 (-55248.998) andoyer -2616650.015 (-55248.997) thomas -2616884.666 (-55253.952) vincenty -2616880.404 (-55253.862) series 884342.911 (18672.371) spherical 47.162 (0.996) projected 47.361
2 geographic 2044950.046 (18246.779) andoyer 2044950.041 (18246.779) thomas -8326487.667 (-74295.984) vincenty -8326396.453 (-74295.171) series -772899.346 (-6896.463) spherical 111.666 (0.996) projected 112.072
4 geographic 1741861.440 (7072.469) andoyer 1741861.429 (7072.469) thomas 6661512.407 (27047.696) vincenty 6596567.215 (26784.000) series 1509817.845 (6130.304) spherical 245.457 (0.997) projected 246.288
8 geographic -8468850.767 (-14662.330) andoyer -8468850.757 (-14662.330) thomas 578.399 (1.001) vincenty -3908.772 (-6.767) series 1972339.431 (3414.760) spherical 575.644 (0.997) projected 577.592
16 geographic 20164620.681 (13140.102) andoyer 20164620.604 (13140.102) thomas 1533.261 (0.999) vincenty 18778.747 (12.237) series -1513872.171 (-986.502) spherical 1529.411 (0.997) projected 1534.586
32 geographic 150640065.989 (32514.002) andoyer 150640065.466 (32514.002) thomas 4638.348 (1.001) vincenty 13597.633 (2.935) series -183191.011 (-39.540) spherical 4617.457 (0.997) projected 4633.083
64 geographic -306584731.272 (-19904.868) andoyer -306584730.490 (-19904.868) thomas 15417.078 (1.001) vincenty 15913.975 (1.033) series 625519.996 (40.612) spherical 15350.553 (0.997) projected 15402.500
128 geographic -694773483.095 (-12455.078) andoyer -694773480.986 (-12455.078) thomas 55817.066 (1.001) vincenty 56655.112 (1.016) series 294780.677 (5.284) spherical 55594.213 (0.997) projected 55782.347
256 geographic 190560276.129 (906.418) andoyer 190560275.143 (906.418) thomas 210317.872 (1.000) vincenty 210709.415 (1.002) series -165180.081 (-0.786) spherical 209525.477 (0.997) projected 210234.526
512 geographic 1178331674.125 (1444.692) andoyer 1178331670.301 (1444.692) thomas 815989.978 (1.000) vincenty 816433.473 (1.001) series -158161.915 (-0.194) spherical 812877.278 (0.997) projected 815628.126
1024 geographic -527570208.535 (-157.729) andoyer -527570208.727 (-157.729) thomas 3346138.595 (1.000) vincenty 3346159.365 (1.000) series 2977999.064 (0.890) spherical 3333501.613 (0.997) projected 3344782.350
2048 geographic 8959153318.996 (674.797) andoyer 8959153288.916 (674.797) thomas 13282138.011 (1.000) vincenty 13282195.756 (1.000) series 13503619.322 (1.017) spherical 13232038.375 (0.997) projected 13276816.420
4096 geographic -15045873900.053 (-287.173) andoyer -15045873851.174 (-287.173) thomas 52414219.343 (1.000) vincenty 52414204.389 (1.000) series 52647596.389 (1.005) spherical 52216378.481 (0.997) projected 52393088.902
8192 geographic -4675164620.233 (-22.239) andoyer -4675164605.840 (-22.239) thomas 210303212.277 (1.000) vincenty 210303227.115 (1.000) series 210596896.520 (1.002) spherical 209509936.500 (0.997) projected 210219019.846
16384 geographic 841593524.763 (1.002) andoyer 841593523.057 (1.002) thomas 840475852.736 (1.000) vincenty 840475842.661 (1.000) series 840286070.892 (1.000) spherical 837302908.768 (0.997) projected 840137824.184
testing for d
1 geographic -590757.714 (-189078.171) andoyer -590757.712 (-189078.170) thomas 0.000 (0.000) vincenty -0.928 (-0.297) series 3.134 (1.003) spherical 3.114 (0.997) projected 3.124
2 geographic -109060.586 (-8726.494) andoyer -109060.586 (-8726.494) thomas 12.474 (0.998) vincenty -2.411 (-0.193) series 12.659 (1.013) spherical 12.455 (0.997) projected 12.498
4 geographic -44830.186 (-896.773) andoyer -44830.186 (-896.773) thomas 50.041 (1.001) vincenty -9.739 (-0.195) series 50.082 (1.002) spherical 49.822 (0.997) projected 49.991
8 geographic -45412.033 (-227.103) andoyer -45412.033 (-227.103) thomas 199.987 (1.000) vincenty -38.992 (-0.195) series 199.974 (1.000) spherical 199.288 (0.997) projected 199.962
16 geographic -2415.338 (-3.020) andoyer -2415.338 (-3.020) thomas 800.173 (1.000) vincenty -155.643 (-0.195) series 800.143 (1.000) spherical 797.151 (0.997) projected 799.849
32 geographic 12972.208 (4.055) andoyer 12972.208 (4.055) thomas 3200.669 (1.000) vincenty 13070.195 (4.085) series 3200.655 (1.000) spherical 3188.605 (0.997) projected 3199.396
64 geographic 10564.452 (0.826) andoyer 10564.452 (0.826) thomas 12802.682 (1.000) vincenty 21033.070 (1.644) series 12802.642 (1.000) spherical 12754.422 (0.997) projected 12797.584
128 geographic 48182.060 (0.941) andoyer 48182.060 (0.941) thomas 51210.939 (1.000) vincenty 50955.078 (0.995) series 51211.032 (1.000) spherical 51017.687 (0.997) projected 51190.336
256 geographic 204964.881 (1.001) andoyer 204964.881 (1.001) thomas 204843.977 (1.000) vincenty 203901.637 (0.996) series 204843.926 (1.000) spherical 204070.754 (0.997) projected 204761.346
512 geographic 819471.930 (1.001) andoyer 819471.930 (1.001) thomas 819375.948 (1.000) vincenty 819152.330 (1.000) series 819375.861 (1.000) spherical 816283.046 (0.997) projected 819045.416
1024 geographic 3277740.667 (1.000) andoyer 3277740.670 (1.000) thomas 3277504.330 (1.000) vincenty 3277381.987 (1.000) series 3277504.419 (1.000) spherical 3265132.427 (0.997) projected 3276181.924
2048 geographic 13109957.429 (1.000) andoyer 13109957.440 (1.000) thomas 13110019.639 (1.000) vincenty 13110025.356 (1.000) series 13110019.757 (1.000) spherical 13060531.615 (0.997) projected 13104729.856
4096 geographic 52440457.103 (1.000) andoyer 52440457.145 (1.000) thomas 52440093.556 (1.000) vincenty 52440113.002 (1.000) series 52440093.394 (1.000) spherical 52242141.106 (0.997) projected 52418938.059
8192 geographic 209761676.873 (1.000) andoyer 209761677.039 (1.000) thomas 209760481.410 (1.000) vincenty 209760493.006 (1.000) series 209760481.321 (1.000) spherical 208968671.714 (0.997) projected 209675922.902
16384 geographic 839047232.375 (1.000) andoyer 839047233.041 (1.000) thomas 839042616.289 (1.000) vincenty 839042622.334 (1.000) series 839042616.172 (1.000) spherical 835875386.941 (0.997) projected 838705402.291
testing for e
1 geographic 189.969 (60.801) andoyer 189.969 (60.801) thomas 3.126 (1.000) vincenty -0.608 (-0.195) series 3.126 (1.000) spherical 3.114 (0.997) projected 3.124
2 geographic 137.277 (10.984) andoyer 137.277 (10.984) thomas 12.503 (1.000) vincenty -16.970 (-1.358) series 12.503 (1.000) spherical 12.455 (0.997) projected 12.498
4 geographic 38.234 (0.765) andoyer 38.234 (0.765) thomas 50.011 (1.000) vincenty 67.352 (1.347) series 50.011 (1.000) spherical 49.822 (0.997) projected 49.991
8 geographic 209.906 (1.050) andoyer 209.906 (1.050) thomas 200.043 (1.000) vincenty 208.040 (1.040) series 200.043 (1.000) spherical 199.288 (0.997) projected 199.962
16 geographic 796.485 (0.996) andoyer 796.485 (0.996) thomas 800.172 (1.000) vincenty 790.812 (0.989) series 800.172 (1.000) spherical 797.151 (0.997) projected 799.849
32 geographic 3199.594 (1.000) andoyer 3199.594 (1.000) thomas 3200.688 (1.000) vincenty 3201.491 (1.001) series 3200.688 (1.000) spherical 3188.605 (0.997) projected 3199.396
64 geographic 12797.747 (1.000) andoyer 12797.747 (1.000) thomas 12802.750 (1.000) vincenty 12804.544 (1.001) series 12802.750 (1.000) spherical 12754.422 (0.997) projected 12797.584
128 geographic 51211.510 (1.000) andoyer 51211.510 (1.000) thomas 51211.000 (1.000) vincenty 51209.779 (1.000) series 51211.000 (1.000) spherical 51017.687 (0.997) projected 51190.336
256 geographic 204845.098 (1.000) andoyer 204845.098 (1.000) thomas 204844.006 (1.000) vincenty 204843.414 (1.000) series 204844.006 (1.000) spherical 204070.754 (0.997) projected 204761.346
512 geographic 819380.579 (1.000) andoyer 819380.580 (1.000) thomas 819376.053 (1.000) vincenty 819376.109 (1.000) series 819376.053 (1.000) spherical 816283.046 (0.997) projected 819045.416
1024 geographic 3277522.367 (1.000) andoyer 3277522.370 (1.000) thomas 3277504.458 (1.000) vincenty 3277504.375 (1.000) series 3277504.457 (1.000) spherical 3265132.427 (0.997) projected 3276181.924
2048 geographic 13110091.522 (1.000) andoyer 13110091.532 (1.000) thomas 13110019.743 (1.000) vincenty 13110019.993 (1.000) series 13110019.739 (1.000) spherical 13060531.615 (0.997) projected 13104729.856
4096 geographic 52440380.626 (1.000) andoyer 52440380.668 (1.000) thomas 52440093.629 (1.000) vincenty 52440094.598 (1.000) series 52440093.610 (1.000) spherical 52242141.107 (0.997) projected 52418938.060
8192 geographic 209761629.539 (1.000) andoyer 209761629.705 (1.000) thomas 209760481.451 (1.000) vincenty 209760485.160 (1.000) series 209760481.384 (1.000) spherical 208968671.714 (0.997) projected 209675922.902
16384 geographic 839047208.844 (1.000) andoyer 839047209.511 (1.000) thomas 839042616.438 (1.000) vincenty 839042619.881 (1.000) series 839042616.309 (1.000) spherical 835875386.940 (0.997) projected 838705402.291
Hi @vissarion, cool, the series is indeed better and perfect for both double and long double. But for float it behaves worse. Maybe we can make series that the default strategy for area calculations? And internally please promote to at least double to avoid large deviations? How is it performance-wise (w.r.t. speed)? That PR is not merged yet I think? What is blocking it?
Maybe we can make series that the default strategy for area calculations? How is it performance-wise (w.r.t. speed)?
Roughly speaking for inverse problems (distance) andoyer is 2x faster than thomas which is 2x faster than vincenty which is faster than GeographicLib (that implements the series formulas we discuss). I should benchmark those formulas for this specific problem though.
That PR is not merged yet I think? What is blocking it?
Apart from the degrees/radians issue discussed above I am ok with merging. There are some comments by @awulkiew unanswered to the PR page though.
I posted some results on accuracy and performance for various strategies here: https://github.com/vissarion/geometry/wiki/Accuracy-and-performance-of-geographic-algorithms
I would like to add one (or more) projection methods for area but I think the one used in the example code above is not the right one, e.g. it uses the bessel ellipsoid instead of WGS84.
@awulkiew what projection method do you recommend for area computation?
@vissarion The choice of projection depends on place on the globe and the usage of the map.
For each region of the globe there are some standard ones, e.g. Barend above used this one: https://epsg.io/28992 for Netherlands. So this is a safe pick if you want to represent coordinates in Netherlands in a standard way.
If you want to calculate areas at various places of the globe and want to compare them then you could choose one of the equal-area projections, e.g.:
- aea - Albers Equal Area (https://proj4.org/operations/projections/aea.html)
- cea - Equal Area Cylindrical (https://proj4.org/operations/projections/cea.html)
- leac - Lambert Equal Area Conic (https://proj4.org/operations/projections/leac.html)
- laea - Lambert Azimuthal Equal Area (https://proj4.org/operations/projections/laea.html)
@awulkiew thanks! I updated the benchmarks. As expected the right projection could be very accurate for short distances but as the distance increases (typically >100m) then projection method becomes inaccurate. Interestingly, there is no sole formula that attains the highest accuracy in all the tested cases.
Returning to the current issue, there is a lot of freedom for the choice of default strategy for area depending on where we want to focus (e.g. do we want good accuracy for short or long distances?). Does it make sense to have more than one default depending on whether a segment is "long" or "short"?