[Proj] Exact transverse Mercator code

Tue Feb 3 05:54:30 EST 2009

Charles & the interested:

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

It is my convention to leave the orientation of the globe as it has been tilted by spherical coordinate transformation. If this were a sphere, it would become transverse by moving the pole of the original sphere to the equator and thence applying the equatorial Mercator formulæ. That puts the north pole along the left-right axis, and that is how I leave it, since I work with much more than just the Mercator and with any spherical coordinate transformation, not just transverse. In point of fact, this being an ellipsoid, there is no actual spherical coordinate transformation going on, but my software is completely general, and so I do not wish to abruptly change convention when the eccentricity edges from 0.0 to 0.000001. 

Presumably the presence of geography under the graticule suffices even for the unwary? — Not that I expect any unwary to be following this conversation!

>What point(s) project to infinity with this projection?

φ = 0.0, λ = ±90°(1±e)

Wallis happens to project to the ellipsoidal stereographic to arrive onto the plane from the ellipsoid. The specific projection to the plane is not important, other than the simpler the better for practical purposes. The point is to arrive at Wallis I, where parallels are evenly spaced on the complex plane.

Define
    Project to complex plane, ψ = X + iΥ
    Wallis I TM, complex ζ = π/2 - X⁻¹ (ψ)
    Wallis II TM, complex ξ = η⁻¹ (ζ)
    Gauß-Krüger, complex F (ξ)

where
    X (φ, λ), Υ (φ, λ) is any conformal mapping to the plane; e.g.,
        X = r cos λ 
        Υ = r sin λ
        r = [(1-sin φ)/cos φ][(1+e sin φ)/(1-e sin φ)]^(e/2).
    X⁻¹ inverts X, but acts on (ψ, 0.0), ψ being the complex "latitude"
    η⁻¹ is inverse parametric latitude, acting on a complex "latitude"
        η⁻¹ (p) = tan⁻¹ [tan (p)/√(1-e²)]
    F is the elliptic integral of the second kind with parameter e.
    e is eccentricity

As usual, the devil is in the details.

Regards,
— daan Strebe

On Feb 1, 2009, at 11:36:46 AM, "Karney, Charles" <ckarney at Sarnoff.com> wrote:
> 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
_______________________________________________
Proj mailing list
Proj at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20090203/d94d54ff/attachment.htm 