[Proj] Troubles with Newton-Raphson inverse projections

[Gerald I. Evenden]

On Sunday 19 October 2008 3:48:28 am OvV_HN wrote:
> One shouldn't laugh too soon.
> With Winkel-Tripel on a sphere with a radius, equal to the major axis of
> the WGS84 ellipsoid, I got the following results.
>
> Starting with: lat = 0.1 deg;  lon = 14 deg (east); lat_0 = 89.9 deg;
> Full round trip forward - inverse needed 286 iterations to get lat/lon back
> within an error of 1e-8 in the iteration loop. The loop needed 579
> iterations for an error of 1e-14.
> So it really needed that lot of iterations.
>
> Two things are rather stupid though.
> 1) How silly is this type of projection with a lat_0 near the pole, being
> interested in locations near the equator? Locations near the equator are
> not extremely silly, because the Winkel-Tripel will be probably used for
> thematic maps of the whole world. But a lat_0 near the pole is a bit
> strange.

I agree but someone always comes along and gets bent out of shape because 
their precious data does not process well in certain areas and explaining the 
problem can be difficult.

Take the case where data contains lat=90, lon=50 and they expect that when 
retrieving the information from a Cartesian system like Hammer that will get 
lon=50 back.  The idea that longitude at the pole is a meaningless item can 
sometime be a difficult concept for some to grasp.  ;-)
It is kinda like trying to retrieve the essence of something that has passed 
the event horizon of a black hole.

> 2) The Newton-Raphson method fails miserably in the known mathematical
> areas. We've already seen this happening in the Transverse Mercator
> projection, performed with complex elliptic functions, and near the
> mathematical boundaries.
> One should use another method in the difficult areas or restrict or warn
> the user for these areas.

Those who develop procedures to evaluate functions like cosine. sine etc. run 
into the same problems.  Arccosine behaves miserably as the argument nears 0 
and that is why Snyder developed the alternative inverse spherical geodesic 
formula to replace the arccosine version with an arcsine version in order to 
retain significant figures for closely spaced points.

This is why I was not surprised when these troubles arose because the shape of 
the projection near the poles indicated potential problems.  If the poles 
were cusp shape the we would probably not have a problem.

The results are not likely to be 8 digit precision near the poles but if we 
can keep the lid on runaway looping and get at least 3 or 4 digits I will 
consider it a success.

I fear those in this group that will expect type double machine 
precision.  :-(
 
> Oscar van Vlijmen
>
>
> ----- Original Message -----
> From: "Gerald I. Evenden" <geraldi.evenden at gmail.com>
> To: "PROJ.4 and general Projections Discussions" <proj at lists.maptools.org>
> Sent: Saturday, October 18, 2008 9:17 PM
> Subject: Re: [Proj] Troubles with Newton-Raphson inverse projections
>
> > On Saturday 18 October 2008 2:10:18 pm OvV_HN wrote:
> >> Some time ago I programmed the inverse of the Winkel-Tripel according to
> >> the Turkish paper.
> >> I found no serious difficulties apart from the following.
> >> The lat_0 must be smaller than for instance 80d absolute and larger than
> >> say 1e-6 deg.
> >> In these circumstances the maximum number of iterations in the inverse
> >> should be smaller than 150 or so.
> >> If a larger value of lat_0 must be used: in one test I came up with 579
> >> iterations at a
> >> lat_0 of 89.9d.
> >
> > LOL, wow!!  I cutoff at 10 iterations and running at a moderately loose
> > tolerance of 10^-8 radians.
>
> ....
>
>
>
>
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj



-- 
The whole religious complexion of the modern world is due
to the absence from Jerusalem of a lunatic asylum.
-- Havelock Ellis (1859-1939) British psychologist