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:

return dist > 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@github.com:awulkiew/geometry
git fetch awulkiew
git checkout feature/geographic
git pull

Regards,
Adam