[Proj] Re: Relationship betweeen transverse Mercator algorithms

Tue Sep 9 15:31:43 EDT 2008

On Sep 9, 2008, at 6:52:27 AM, "Charles Karney" <ckarney at sarnoff.com> wrote:

>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

And then there's Dozier's, which is not a series expansion. I haven't looked much at his paper and I don't have it in front of me, but as I recall, it's Lee's development through Jacobian elliptic functions. Naturally it's much less efficient than series expansion, but it gives a few bits less than machine precision over the complete globe — with some caveats about a few numerically problematic regions, if I recall.

There is also Wallis (unpublished) that deviates a little from steps 1-3 below. Wallis projects directly to the plane first via the ellipsoidal polar stereographic. Then, treating the (x,y) results as the complex plane, numerically computes the inverse of of the stereographic. For the central meridian this merely gives the original coordinates; and in general it gives the complex amplitude needed for the elliptic integral that computes the meridian segment length. In practice Wallis takes advantage of some minor efficiencies that fall out of using the colatitude as the elliptic integral's amplitude, and also in practice, the numerical difficulties must be carefully  managed in some regions.

Because the ellipsoidal polar stereographic is so structurally similar to using the conformal latitude, the numerically troublesome areas are really all the same in the end. It's all about the conformal latitude function inverse...

Regards,
-- daan Strebe

On Sep 9, 2008, at 6:52:27 AM, "Charles Karney" <ckarney at sarnoff.com> wrote:
From:   "Charles Karney" <ckarney at sarnoff.com>
Subject:    [Proj] Relationship betweeen transverse Mercator algorithms
Date:   September 9, 2008 6:52:27 AM PDT
To: Proj at lists.maptools.org
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
_______________________________________________
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/20080909/68b895ec/attachment.html