Boost logo

Geometry :

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:

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

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:
>>>> 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:
>> 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:
> also in a PR:
> So e.g. to run 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, Barend

Geometry list run by mateusz at