[Proj] C++ code for exact Transverse Mercator Projection

[Charles Karney]

First of all, I posted maxima code for the Transverse Mercator
Projection in

    http://lists.maptools.org/pipermail/proj/2008-September/003729.html

There was a 1-letter error in ellint.mac which will have resulted in
incorrect results.  The error is in the smallest term of one of the
expansions of the elliptic integrals and most like the resulting error
is tiny (1e-40 or so).  Anyway, the corrected version is now posted.

Secondly, I've posted a C++ implementation of the Transverse Mercator
projection in

    http://charles.karney.info/geographic/TransverseMercatorExact.tgz

Accuracy is about 12 nm over the whole globe.  Unpack the tar file
somewhere convenient, change to the GeographicLib directory, and read
the comment at the top of TransverseMercatorTest.cpp for how to compile
the test program.  The code compiles under Linux and Visual Studio.

This is just a straightforward implementation of transformations given
by Lee in his 1976 monograph.  The resulting code is quite fast.  I
haven't done detailing timing on it.  However I estimate that it's no
more than a "few" times slower than a standard series solution (a la
Krueger).

Here the blurb taken from TransverseMercatorExact.cpp

Implementation of the Transverse Mercator Projection given in

   L. P. Lee,
   Conformal Projections Based on Elliptic Functions,
   (B. V. Gutsell, Toronto, 1976), 128pp.,
   ISBN: 0919870163.

   [Also appeared as:
   Monograph 16, Suppl. No. 1 to Canadian Cartographer, Vol 13.]

The method entails using the Thompson Transverse Mercator as an
intermediate projection.  The projections from the intermediate
coordinates to [phi, lam] and [x, y] are given by elliptic functions.
The inverse of these projections are found by Newton's method with a
suitable starting guess.

The same basic technique was used by

  J. Dozier,
  Improved Algorithm for Calculation of UTM and Geodetic Coordinates
  NOAA Technical Report NESS 81, Sept. 1980

The differences here are:

* Use of real instead of complex arithmetic.  This uses the separation
  into real and imaginary parts given by Lee and allows better control
  over branch cuts.

* Different (and possibly better?) algorithms for certain elliptic
  functions.

* Better starting guesses for Newton's method.  Dozier's guesses lead to
  the wrong solution is some cases.  This disqualifies his method for
  routine use.

This method gives the correct results for forward and reverse
transformations with the proviso that the reverse transformation may not
yield a sensible result if the input [x,y] is too far outside the set of
[x,y] obtained from the forward transformation.  The maximum error is
about 12nm (ground distance) for the forward and reverse
transformations.  The error in the convergence is 3e-9", the relative
error in the scale is 1e-11.  The method is "exact" in the sense that
the errors are close to the round-off limit and that no changes are
needed in the algorithms for them to be used with reals of a higher
precision (e.g., long double).

This implementation and notation closely follows Lee, with the following
exceptions:

    Lee    here    Description
    x/a    xi      Northing (unit Earth)
    y/a    eta     Easting (unit Earth)
    s/a    sigma   xi + i * eta
    y      x       Easting
    x      y       Northing
    k      e       eccentricity
    k^2    mu      elliptic function parameter
    k'^2   mv      elliptic function complementary parameter
    m      k       scale

Minor alterations have been made in some of Lee's expressions in an
attempt to control round-off.  For example atanh(sin(phi)) is replaced
by asinh(tan(phi)) which maintains accuracy near phi = pi/2.  Such
changes are noted in the code.

Loose ends:

Testing only done for WGS84 eccentricity.  Things that might go wrong
are other eccentricities: (1) Failure to converge most likely for much
higher eccentricity.  This will require refining the starting guesses
for the Newton's iterations.  (2) Failure with a sphere.  The basic
formulation carries over to the sphere with no problem.  Some attention
might need to be paid to the treatment of the singularity for phi=0,
lam=90.

The singularity at phi=90 is handled safely.  There's another
singularity (with the intermediate projection) at phi=0, lam=90*(1-e).
This is handled by using a Taylor expansion about the singularity.  This
gives a good enough starting guess for Newton's method to converge.
However detailed testing in the immediate neighborhood of the
singularity has not been done.  If there is a problem it can be handled
easily by treating the singularity specially.

The initial guesses for Newton's method are a little ad hoc.  Probably
better guesses can be used and so one or more iterations of Newton's
method can be skipped.

An analysis of the reverse transformation of "out-of-bounds" points has
not been done.  In particular it would be nice to determine the limits
of the correct reverse transformation.  These limits will of course
depend sensitively on the initial guess for Newton's method for the
reverse transformation.  (The out-of-bounds points in the "opposite"
hemisphere.)

-- 
Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

URL: http://charles.karney.info
Tel: +1 609 734 2312
Fax: +1 609 734 2662