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-25 18:16:20

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

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

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

Regards,
Adam