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

Subject: Re: [geometry] Support for geographic coordinate system
From: Barend Gehrels (barend_at_[hidden])
Date: 2014-11-26 14:39:17

Hi Adam,

Adam Wulkiewicz wrote On 26-11-2014 0:16:
> Hi,
>
> Barend Gehrels wrote:
>> Adam Wulkiewicz wrote On 23-11-2014 13:54:
>>> Barend Gehrels wrote:
>>>> Adam Wulkiewicz wrote On 7-11-2014 0:00:
>>>>> Barend Gehrels wrote:
>>>>>> Adam Wulkiewicz wrote On 6-11-2014 17:06:
>>> <snip>
>>>>>>>
>>>>>>>>
>>>>>>>> Basically, to support this CS we need a distance strategy and 
>>>>>>>> side strategy. They're already implemented (more or less) 
>>>>>>>> though in the extensions so not officially released yet. To 
>>>>>>>> release them we need to think about the interface first. The 
>>>>>>>> geographical strategies should have an interface similar to the 
>>>>>>>> one used in spherical strategies.
>>>>>>>>
>>>>>>>> GEOGRAPHICAL:
>>>>>>>>
>>>>>>>> - a spherical side strategy is just used by default which I 
>>>>>>>> guess should be changed. I think the most precise would be a 
>>>>>>>> strategy finding the side using the geographic courses found 
>>>>>>>> using the inverse Vincenty's formula but this requires some 
>>>>>>>> testing
>>>>>>>
>>>>>>> As far as I know, the spherical side strategy works for the 
>>>>>>> geographic Earth too, please indicate if that is not the case.
>>>>>>>
>>>>>>> See also
>>>>>>> http://barendgehrels.blogspot.nl/2011/06/spherical-side-formula.html 
>>>>>>>
>>>>>>>
>>>>>>> But I might be wrong, it is good to research this.
>>>>>>
>>>>>> I feel that Vincenty's formula should give different results in 
>>>>>> some edge cases because the shape of geodesic on ellipsoid is 
>>>>>> different than on a sphere but I must perform some tests to be sure.
>>>>
>>>> I'm curious about this. 
>>>
>>> I've made some quick test where I compared the result of spherical 
>>> side formula and a side found by comparison of azimuths calculated 
>>> using Vincenty's formula. The method comparing azimuths is very 
>>> simple and probably not good enough to be released nevertheless the 
>>> error is too small to be seen in this case. The difference between 
>>> spherical and geographical geodesic seems to be a lot greater.
>>
>> 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.

> 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...

>
>> I cannot completely read that from your story, but the 
>> color-descriptions indicate the Vincenty azimuth comparisons are OK? 
>> It looks good indeed.
>>
>
> Most importantly they give different results. The difference is marked 
> with bright colors. Some users may be ok with spherical calculations 
> some may not.
>
>>> It's also possible that I made something wrong. In that case don't 
>>> hesitate to point it out :)
>>
>> bgd::vincenty_inverse<double> vi2(lon_s1 * bg::math::d2r, lon_s1 * 
>> bg::math::d2r,
>> the second lon should be lat. Makes no difference for this test 
>> because both are 51, but for other tests it should be changed.
>>
>> you skip collinear cases?
>>
>
> 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.

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

Thanks, Barend