[Proj] Relationship betweeen transverse Mercator algorithms

[Charles Karney]

In the discussion of transverse Mercator projections, there are several
"competing" algorithms.  However, it is probably a good idea to
understand just how similar they are.  In

    http://lists.maptools.org/pipermail/proj/2006-June/002296.html

(dated 2006-06-13), Oscar van Vlijmen lists 3 "accurate" approximations
to the transverse Mercator projection

    [fi] http://docs.jhs-suositukset.fi/jhs-suositukset/JHS154/JHS154.pdf
    [fr] http://www.ign.fr/telechargement/MPro/geodesie/CIRCE/NTG_76.pdf
    [se] http://www.lantmateriet.se/upload/filer/kartor/geodesi_gps_och_detaljmatning/geodesi/Formelsamling/Gauss_Conformal_Projection.pdf

and he states

    "Each follow a slightly different route, but the differences in the
     results are small."

This comment is exactly right.  I might strengthen the statement as
follows:

    "The algorithms are all the same differing only in minor details."

To van Vlijmen's list can be added

    [dk] Engsager-Poder algorithm implemented in proj_etmerc

and my restatement applies to this algorithm too.

What is the basic forward algorithm (phi,lam) -> (x,y)?

(1) convert geodetic latitude, phi, to conformal latitude, beta
(2) convert (beta,lam) to spherical TM (x',y')
(3) apply a conformal mapping to (x',y') to get (x,y)

Note that (apart from scale factors)

  y' = conformal latitude, beta
  y = rectifying latitude

so the mapping in (3) is serves just to rescale y and the x coordinate
follows by the requirements of conformality.

The reverse algorithm (x,y) -> (phi,lam) just reverses these steps
(3') apply a conformal mapping to (x,y) to get (x',y')
(2') convert spherical TM (x',y') to (beta,lam)
(1') convert conformal latitude, beta, to geodetic latitude, phi.

So how do the algorithms differ:

(1) and (1'):

    [fi] and [fr] use exact formulas.  The formulas are different but
    mathematically equivalent.  This approach necessitates an iterative
    solution method in (1')

    [se] and [dk] use series for (1) and its reversion for (1').  [se]
    uses includes terms up to e^8 while [dk] includes terms up to e^10.

    The other difference is in the choice of small parameter used for
    the expansion.  [se] uses e^2 = f*(2-f), [dk] uses n = f/(2-f).  The
    use of different expansion parameters changes the numerical values
    at a given order.

    [se] also writes the series as a power series in sin(phi) but this
    is easily converted to a series in sin(2*k*phi).

    The use of a series expansion here breaks the conformality of the
    algorithm and for this reason I prefer the [fi] and [fr] approaches.
    However the e^8 (and certainly the e^10) series suffice to give
    accuracy close to round-off for doubles.

(2) and (2')

    All approaches use the exact transformations for thess steps.  Again
    the formulas are different (and this may affect the numerical
    accuracy) but mathematically equivalent.

(3) and (3')

    All approaches use a series and its reversion for these steps.

    [fi] and [se] use same series in n accurate to n^4 (e^8).

    [dk] extends these series out to n^5 (e^10).

    [fr] uses a series in e^2 accurate to e^8.  In addition it folds the
    overall scale normalization into the series.

    It is easy to convert from the series used by [fi], [se], and [dk]
    (truncated to n^4) to the ones used by [fr] by substituting for n in
    terms of e^2 multiplying by the normalizing factor (called ahat in
    [se] and A1 in [fr]) and expanding in e^2.

    The use of a series and its reversion means that the forward and
    backward transformation are not exact inverses of one another.  The
    discrepancy here is of the same magnitude as the error in the
    forward transformation.

So there you have it.  The methods are all essentially equivalent.

[dk] has the accuracy edge because its series are truncated at n^5
instead of n^4.  However any of the other methods can leapfrog over [dk]
by using more terms in the series for (3) and (3') (and (1) and (1') in
the case of [se]).  I posted series accurate to n^8 in

   http://charles.karney.info/geographic/UTM-fi.txt

These can be directly used in [fi] and with minor modifications are
applicable to the other methods too.  In the same file, I also give a
table of errors for the different order approximations (n^4, n^5, n^6,
n^7, and n^8).  [The maxima code that I used to compute these series is
relatively short (40 lines or so), and I will post this at some point.]

We are then just left with relatively pedestrian numerical issues which
are important in the control of round-off:

  * use of Horner representation of polynomials
  * stable summation of trig series
  * recognize that asin(tanh(Q)) is inaccurate if Q is large and use
    the mathematically equivalent atan(sinh(Q)) instead
  * and so on.

The use of these series expansions in (3) and (3') limits its
applicability.  The transformation from conformal to rectifying latitude
when extended to the complex domain has a mild singularity a phi = 0,
lam = pi/2*(1-e) = 82.62d which is enough to make the series divergent
there.  Indeed the zone of 1mm accuracy for the n^8 series is limited to
72 degrees from the central meridian.

Finally, all the methods which I have labeled by country codes derive
from the work of 2 Germans, Gauss and Krueger.  The series (accurate to
n^4) needed for (3) and (3') are given by equations 26* (p. 18) and 41
(p. 21) in Krueger's 1912 paper

    http://www.gfz-potsdam.de/bib/pub/digi/krueger2.pdf

So perhaps they all should be called [de].

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