Geometry Boost List: Re: [geometry] Support for geographic coordinate system

Subject: Re: [geometry] Support for geographic coordinate system
From: Adam Wulkiewicz (adam.wulkiewicz_at_[hidden])
Date: 2014-11-26 15:25:28

Hi Barend,

Barend Gehrels wrote:
> Adam Wulkiewicz wrote On 26-11-2014 0:16:
>> Barend Gehrels wrote:
>>> Thanks! But...
>>> What is the conclusion?
>>> - SSF can only be used for spherical and not for geographic (all 
>>> non-spheres)?
>>
>> It can be used but will give wrong results for spheroids.
>
> If I ask if it can be used, I mean of course with correct results. So 
> your answer is no.
>

Well, yes and no. It's more about choosing a default for geographic CS. 
Similar but not the same - distance Andoyer vs Vincenty vs haversine. 
But in the case of distance often a distance between far Points is 
measured. On the other hand, in the case of side calculation, I'm 
guessing that it isn't very common to define Polygons containing 
150km-long segments.

>
>> SSF gives the same results as Vincenty for sphere/flattening=0.
>>
>>> - the method comparing azimuths (you mentioned is probably not good 
>>> enough) is not sufficient?
>>>
>>
>> It's because the further the Point the calculation becomes less 
>> accurate. A "real" cross-track distance should probably be calculated 
>> and compared with EPS (if needed), similar to side_by_cross_track (or 
>> side_by_azimuth). However even there the distance nor radius isn't 
>> taken into account and without it we can't calculate a value of XTD. 
>> But maybe it's sufficient to do it this way... At least for doubles, 
>> I'm guessing that for floats it'd be different.
>
> So for points up to about 0.6 meter left or right of a segment of 153 
> km, there is a difference in ssf and Vincenty. And you state that 
> Vincenty is correct and ssf is wrong, and you mention the Vincenty 
> calculation is less accurate (but maybe not in this case). If doubles, 
> for floats it is unknown. Besides that, Vincenty is very good, but it 
> is an approximation.
>
> Just a question, I'm just curious - how do you know for sure then that 
> the Vincenty approximation is correct here, and the ssf results are 
> wrong? It is not that easy to verify...
>

Of course in all formulas some approximation of a globe is used. My 
assessment is based on a fact that SSF assume thet a globe is a sphere 
and Vincenty that it's a spheroid. So the approximation is closer to 
reality. And the greater the flattening, the greater the difference 
between SSF and Vincenty. I'm not considering precision, numerical 
errors, stability etc.
As for the Vincenty formula here: http://www.icsm.gov.au/gda/gdav2.3.pdf 
states:
"Vincenty's formulae (Vincenty, 1975) may be used for lines ranging from 
a few cm to nearly 20,000 km, with millimetre accuracy."
Calculation done with SSF is probably also precise but it uses a sphere 
model, that's all.

And as mentioned above, I feel that the above is more important for 
calculating distances. Often the distance between far Points is 
measured. I don't know if we need such great precision for the sides 
calculation. Who defines a Polygon containing 153km-long segments? And 
is an error of 0.6m a big deal in this case?

Well, maybe on Saturn which has flattening=0.1 and is bigger ;)

>
>> Do you mean, when a point is on a segment? The epsilon is so small 
>> that all test points are missing the segment (at least for double).
>>
>> Btw, in spherical_side_formula EPS isn't used in the final result 
>> calculation:
>>
>> returndist>zero?1
>>       :  dist  <  zero  ?  -1
>>       :  0;
>>
>> AFAIU dist should be compared with 0 using math::equals(). Do you 
>> agree? Or is there a reason why it's implemented like this?
>
> Yes there is a reason: in the past we used it using FP-input a lot and 
> had to fight with robustness issues. I remember that it has been using 
> math::equals() but that is gone.
> However, all these robustness problems are gone for intersections so 
> we might add the math::equals again. Thanks for pointing this out.
>
>>
>>>
>>>> Here is the code and the results:
>>>> https://github.com/awulkiew/test-geoside
>>>
>>> Can you give a coordinate-pair where the deviation is large 
>>> (probably easy to read from the graph)?
>>>
>>
>> For the Earth the width of the erroreous part around lon=11.3, 
>> lat=17.3 is more or less 1.7 km but this is for this very long 
>> segment (-51 -51, 51 51). For a 1 deg segment (11 17, 12 18) (153km 
>> long) the width of the erroreous fragment is around 0.6m. So the 
>> difference seems to be relativaly small. I didn't check it for 
>> coordinates closer to poles.
>
> I assume lon=11,lat=17 and similarly. I tried that using a third point 
> somewhere in the middle.

Yes, (lon1 lat1, lon2 lat2) in degrees. I'm guessing that it'd be hard 
to find a point of difference becuase a segment is quite "curvy" when 
raw coordinates are considered, as shown here: 
https://github.com/awulkiew/test-geoside.

> However, vincenty_inverse is not in the library - you must use some 
> future improvement... Can you point where I can find that?

It's in my fork awulkiew/geometry branch feature/geographic:

https://github.com/awulkiew/geometry/tree/feature/geographic

also in a PR:

https://github.com/boostorg/geometry/pull/183

So e.g. to run https://github.com/awulkiew/test-geoside besides GLUT 
installed in the system (e.g. freeglut-dev package AFAIR) you'd have to 
check out this branch from my fork, e.g.:

In BG local clone>
git remote add awulkiew git_at_[hidden]:awulkiew/geometry
git fetch awulkiew
git checkout feature/geographic
git pull

Regards,
Adam