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 17:22:46

Hi Adam,

Adam Wulkiewicz wrote On 26-11-2014 21:25:
> 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 

Ah, so you did *not *read my blog, and apparently also not my statements 
that SSF can be used on a spheroid.

So I repeat it again, there it goes:

http://barendgehrels.blogspot.nl/2011/06/spherical-side-formula.html

Quoting, a.o.:
"Summary: Using some high-school mathematics I presented an algorithm 
and a formula to calculate at which side a point is with resepect to a 
segment, on a sphere or on the Earth"

Beside the summary, it has a whole paragraph about the Ellipsoid.

Because this is exact, and Vincenty is an approximation (close to 
reality, but still, an approximation), I now just assume that these 
results are correct and the Vincenty approximation is off within ranges 
of 0.6.

I assume that, until it is proven that it is the other way round...

If SSF was only applicable for 100% spheres, the error would be much 
larger. Compare Haversine/Vincenty, it can be off many kilometers.

Note: it is different from distance, where an exact mathematical formula 
for ellips or spheroid just does not exist.

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

Nope.

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

Yes we do. Side is crucial in all intersection operations. Left or right 
is a very large difference.

> 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

Thanks.

Regards, Barend