[geotk] Question about behaviour of geotoolkit version 3.19

[Martin Desruisseaux]

Hello David

Welcome :-)

The two points used in the test file (38°57,9?N 121°32,3?E   and 
38°29,0?S 58°49,0?W) are close to antipodal points: they are 179.6° 
apart in longitude, and the difference in latitude (except for sign) is 
only 0.5°. Calculation of orthodromic distance is know to have 
convergence problems for points close to antipodal.

The reason why the loop stops after 12 iterations is for avoiding 
infinite loops. The value 12 is arbitrary. The real stop condition is 
the check for tolerance value (e.g. while (xy >= TOLERANCE_1)). This is 
a common practice in numerical computations. The idea is start with an 
approximation, then iteratively add a correction. The correction should 
be smaller and smaller after each iteration, and we stop when the 
correction is smaller than the desired accuracy. However in some 
particular cases, we fail to reach the accuracy target. This may be 
caused for example by a subtraction of two very close numbers, if the 
difference between those numbers is smaller than the IEEE 754 'double' 
accuracy (this problem is known as "/catastrophic cancellation/" [1]). 
So numerical algorithms typically put a limit in the number of 
iterations (12 in this particular case), meaning "if we didn't reached 
the accuracy target after 12 iterations, it is probably not going to work".

In the case of your test points, I saw that we get a convergence after 
59 iterations. So we could increase the iteration limit from 12 to 60 
for instance, but it would only delay the problem. If we get an answer 
for those two points (and I don't know what would be the accuracy of 
that answer if intermediate calculations are close to IEEE 754 precision 
limits), we would still have a failure for points closer to antipodes. I 
didn't tried, but my guess is that the number of iterations would grow 
exponentially as we get closer and closer to antipodal points.

There is some possible alternatives:

1) We could increase the iteration limit, but we would probably need to 
know how close to antipodal points you application may need to go.

2) The following code also compute the same kind of distance, but using 
a simpler algorithm. This method is more robust (it works for the points 
given in your test case) but I have not done a comparison of their 
accuracy. Note also that while more robust, this code will still fail 
for points too close to antipodal - it will just get closer:

    DefaultEllipdoid ellipsoid = DefaultEllipsoid.WGS84; // Choose your ellipsoid here.
    double distance = ellipspoid.orthodromicDistance(sourceLocLon, sourceLocLat, targetLocLon, targetLocLat);


3) If we replace WGS84 in the above code by SHERE, then the code will 
use yet simpler formulas. It will use spherical formulas instead than 
ellipsoidal ones. The spherical formulas should never fail, even for 
fully antipodal points. However I think that the difference between 
ellipsoidal and spherical formulas may be as much as a few hundred of 
kilometres in the worst case.

     Regards,

         Martin


[1] http://en.wikipedia.org/wiki/Loss_of_significance



Le 27/11/13 11:08, Cemernek, David a écrit :
> My name is David Cemernek and I am working for the company Infonova. 
> We are using the geotoolkit (version 3.19) in one of our projects for 
> calculating distances between locations.
>
> For a pair of locations, we get an exception while calculating the 
> distance for these two locations. I am sorry to bother you but I 
> downloaded the source code of geotoolkit and I do not know why the 
> attached program is producing an exception.
>
> Program to reproduce error: TestGeotoolkit.java (also with .txt ending 
> in case an antivirus software deletes it J
>
> Output of the program: Output_TestGeotoolkit.txt
>
> Due to that fact I don't know what the loop in line 1062 is actually 
> doing, I don't want to open an Jira issue, because I don't know why 
> this exception is thrown when the loop runs more than 12 times? 
> Because I think that the distance between the given source and target 
> location is below 20.000 km, and there for should be correct (or is 
> the max distance in the used ellipsoid below this limit?)
>
> I hope you can help me.
>
> Yours sincerely, David
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/geotoolkit/attachments/20131128/d7d2cb78/attachment.html>