[Proj] Exact transverse Mercator code

[Karney, Charles]

> From: strebe <strebe at aol.com>
> Date: Sunday, February 01, 2009 04:23
>
> I found your work very interesting. Thank you for the excellent depth.

Why, thank you.

> I've posted here and image of yet another ellipsoidal transverse
> Mercator variation:
>      http://www.mapthematics.com/Wallis1TransverseMercator.tif

By the way, why plot TM with grid north to the left instead of up?  It
seems to be yet another possible source of confusion or the unwary.

> This one's properties are:
>      Conformal (of course)
>      Straight central meridian (of course)
>      Parallels are equally spaced.
>
> ... The projection is cropped at a distance of 4 arc seconds around
> the singularities; in point of fact the projection is infinite in
> extent.

What point(s) project to infinity with this projection?  (And to answer
Gerald's question, if the parallels are equally spaced then the scale on
the central meridian must vary.)

The beauty of Lee's exposition (following Thompson, I suppose) is the
economical notation (when written in complex form).  Thus:

Define
    Mercator coordinates, complex zeta = psi + i lambda
    Thompson TM,          complex w
    Gauss-Krueger TM,     complex s
    Gauss-Schreiber TM,   complex z

Then Lee gives (with k = e):
    (9.4)   psi  = atanh(sin(phi)) - k atanh(k sin(phi))
    (54.18) zeta = atanh( sn(w)  ) - k atanh(k  sn(w)  )
    (12.3)  zeta = atanh(sin(z)  )
    (55.8)  s    = ed(w) = integral k'^2 nd(w)^2

Lee's equations (9.4), (54.18), (55.8) completely specify the
Gauss-Krueger TM.  And presumably something similar can be written to
specify Wallis' development.

The remaining tasks are:

(1) separating the equations into real and imaginary components,
(2) inverting these equations,
(3) computing derivatives to give the convergence and scale,
(4) recasting the equations for numerical accuracy.

These are non-trivial, but are merely crank-turning exercises.

(1) is not strictly necessary (but Lee provides the separated
equations).  You could just implement the complex arithmetic.  But you
are then left with the task of checking that right roots are taken with
the inverse trig and inverse hyperbolic functions.

For (2), I used Newton's method with suitable starting points, which are
typically found using low-order expansions about special points.  I
treat the singularity at lat = 0, lon = 90 (1 - e) by retaining the
starting point as the solution if the next term in the expansion is
negligible.

> I have not looked in detail at how our results compare for the
> Gauss-Krueger, though spot-checking several points showed we agree to
> sub-micrometer precision.

For systematic testing, I use the dataset available at

    http://charles.karney.info/geographic/TMcoords.dat.gz

This is a random collecting of points (geographic and TM) which I read
into testing code to compute the errors for the forward and backward
transformation.  There are extra points around potential singularities
(real or numerical).

The results are read into Matlab to analyse these with, e.g.,

    max(err(lat0 >= 0))
    scatter(x0, y0, 2, log10(err));
    and so on...

This approach was very helpful in identifying areas where there was a
loss of accuracy due to cancellation.

Incidentally, I also tested the accuracy of TransverseMercatorExact with
"double" changed to "long double" (64-bit fraction with g++).  The
errors went from 8 nm (double) to 4 pm (long double) consistent with the
extra 11 bits of precision.  (The use of numeric_limits ensured that all
the tolerances change appropriately.)

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