
Geometry : 
Subject: Re: [geometry] Support for geographic coordinate system
From: Barend Gehrels (barend_at_[hidden])
Date: 20141126 17:22:46
Hi Adam,
Adam Wulkiewicz wrote On 26112014 21:25:
> Hi Barend,
>
> Barend Gehrels wrote:
>> Adam Wulkiewicz wrote On 26112014 0:16:
>>> Barend Gehrels wrote:
>>>> Thanks! But...
>>>> What is the conclusion?
>>>>  SSF can only be used for spherical and not for geographic (all
>>>> nonspheres)?
>>>
>>> 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 150kmlong 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" crosstrack 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/sphericalsideformula.html
Quoting, a.o.:
"Summary: Using some highschool 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 153kmlong 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 FPinput 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/testgeoside
>>>>
>>>> Can you give a coordinatepair 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/testgeoside.
>
>> 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/testgeoside besides GLUT
> installed in the system (e.g. freeglutdev 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
Geometry list run by mateusz at loskot.net