<html><body name="Mail Message Editor">On Sep 9, 2008, at 6:52:27 AM, "Charles Karney" <ckarney@sarnoff.com> wrote:<div><br><div><span class="Apple-style-span" style="font-family: -webkit-monospace; font-size: 11px; ">>In the discussion of transverse Mercator projections, there are several<br>>"competing" algorithms. However, it is probably a good idea to<br>>understand just how similar they are. In<br><br></span></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>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...</div><div><br></div><div>Regards,</div><div>-- daan Strebe</div><div><br></div><div><br></div><br>On Sep 9, 2008, at 6:52:27 AM, "Charles Karney" <ckarney@sarnoff.com> wrote:<br><blockquote style="padding-left: 5px; margin-left: 5px; border-left-width: 2px; border-left-style: solid; border-left-color: blue; color: blue; "><span class="Apple-style-span" style="border-collapse: separate; color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: 2; text-align: auto; text-indent: 0px; text-transform: none; white-space: normal; widows: 2; word-spacing: 0px; -webkit-border-horizontal-spacing: 0px; -webkit-border-vertical-spacing: 0px; -webkit-text-decorations-in-effect: none; -webkit-text-size-adjust: auto; -webkit-text-stroke-width: 0; "><div style="width: 100%; "><div id="felix-mail-header-block" style="color: black; background-color: white; border-bottom-width: 1px; border-bottom-style: solid; border-bottom-color: silver; padding-bottom: 1em; margin-bottom: 1em; width: 100%; "><table border="0" cellpadding="1" cellspacing="1" width="100%"><tbody><tr><td width="70px" style="font-family: 'Lucida Grande'; font-size: 8pt; color: gray; text-align: right; vertical-align: top; font-weight: bold; "><span>From:</span></td><td style="font-family: 'Lucida Grande'; font-size: 8pt; color: black; text-align: left; vertical-align: top; padding-left: 5px; "><span title=""Charles Karney" <ckarney@sarnoff.com>">"Charles Karney" <ckarney@sarnoff.com></span></td></tr><tr><td width="70px" style="font-family: 'Lucida Grande'; font-size: 8pt; color: gray; text-align: right; vertical-align: top; font-weight: bold; "><span>Subject:</span></td><td style="font-family: 'Lucida Grande'; font-size: 8pt; color: black; text-align: left; vertical-align: top; padding-left: 5px; "><span style="font-weight: bold; ">[Proj] Relationship betweeen transverse Mercator algorithms</span></td></tr><tr><td width="70px" style="font-family: 'Lucida Grande'; font-size: 8pt; color: gray; text-align: right; vertical-align: top; font-weight: bold; "><span>Date:</span></td><td style="font-family: 'Lucida Grande'; font-size: 8pt; color: black; text-align: left; vertical-align: top; padding-left: 5px; "><span>September 9, 2008 6:52:27 AM PDT</span></td></tr><tr><td width="70px" style="font-family: 'Lucida Grande'; font-size: 8pt; color: gray; text-align: right; vertical-align: top; font-weight: bold; "><span>To:</span></td><td style="font-family: 'Lucida Grande'; font-size: 8pt; color: black; text-align: left; vertical-align: top; padding-left: 5px; "><span title="Proj@lists.maptools.org">Proj@lists.maptools.org</span></td></tr></tbody></table></div><div id="felix-mail-content-block" style="color: black; background-color: white; width: 100%; "><div style="font-family: monospace; color: black; background-color: white; font-size: 8pt; ">In the discussion of transverse Mercator projections, there are several<br>"competing" algorithms. However, it is probably a good idea to<br>understand just how similar they are. In<br><br>http://lists.maptools.org/pipermail/proj/2006-June/002296.html<br><br>(dated 2006-06-13), Oscar van Vlijmen lists 3 "accurate" approximations<br>to the transverse Mercator projection<br><br>[fi] http://docs.jhs-suositukset.fi/jhs-suositukset/JHS154/JHS154.pdf<br>[fr] http://www.ign.fr/telechargement/MPro/geodesie/CIRCE/NTG_76.pdf<br>[se] http://www.lantmateriet.se/upload/filer/kartor/geodesi_gps_och_detaljmatning/geodesi/Formelsamling/Gauss_Conformal_Projection.pdf<br><br>and he states<br><br>"Each follow a slightly different route, but the differences in the<br>results are small."<br><br>This comment is exactly right. I might strengthen the statement as<br>follows:<br><br>"The algorithms are all the same differing only in minor details."<br><br>To van Vlijmen's list can be added<br><br>[dk] Engsager-Poder algorithm implemented in proj_etmerc<br><br>and my restatement applies to this algorithm too.<br><br>What is the basic forward algorithm (phi,lam) -> (x,y)?<br><br>(1) convert geodetic latitude, phi, to conformal latitude, beta<br>(2) convert (beta,lam) to spherical TM (x',y')<br>(3) apply a conformal mapping to (x',y') to get (x,y)<br><br>Note that (apart from scale factors)<br><br>y' = conformal latitude, beta<br>y = rectifying latitude<br><br>so the mapping in (3) is serves just to rescale y and the x coordinate<br>follows by the requirements of conformality.<br><br>The reverse algorithm (x,y) -> (phi,lam) just reverses these steps<br>(3') apply a conformal mapping to (x,y) to get (x',y')<br>(2') convert spherical TM (x',y') to (beta,lam)<br>(1') convert conformal latitude, beta, to geodetic latitude, phi.<br><br>So how do the algorithms differ:<br><br>(1) and (1'):<br><br>[fi] and [fr] use exact formulas. The formulas are different but<br>mathematically equivalent. This approach necessitates an iterative<br>solution method in (1')<br><br>[se] and [dk] use series for (1) and its reversion for (1'). [se]<br>uses includes terms up to e^8 while [dk] includes terms up to e^10.<br><br>The other difference is in the choice of small parameter used for<br>the expansion. [se] uses e^2 = f*(2-f), [dk] uses n = f/(2-f). The<br>use of different expansion parameters changes the numerical values<br>at a given order.<br><br>[se] also writes the series as a power series in sin(phi) but this<br>is easily converted to a series in sin(2*k*phi).<br><br>The use of a series expansion here breaks the conformality of the<br>algorithm and for this reason I prefer the [fi] and [fr] approaches.<br>However the e^8 (and certainly the e^10) series suffice to give<br>accuracy close to round-off for doubles.<br><br>(2) and (2')<br><br>All approaches use the exact transformations for thess steps. Again<br>the formulas are different (and this may affect the numerical<br>accuracy) but mathematically equivalent.<br><br>(3) and (3')<br><br>All approaches use a series and its reversion for these steps.<br><br>[fi] and [se] use same series in n accurate to n^4 (e^8).<br><br>[dk] extends these series out to n^5 (e^10).<br><br>[fr] uses a series in e^2 accurate to e^8. In addition it folds the<br>overall scale normalization into the series.<br><br>It is easy to convert from the series used by [fi], [se], and [dk]<br>(truncated to n^4) to the ones used by [fr] by substituting for n in<br>terms of e^2 multiplying by the normalizing factor (called ahat in<br>[se] and A1 in [fr]) and expanding in e^2.<br><br>The use of a series and its reversion means that the forward and<br>backward transformation are not exact inverses of one another. The<br>discrepancy here is of the same magnitude as the error in the<br>forward transformation.<br><br>So there you have it. The methods are all essentially equivalent.<br><br>[dk] has the accuracy edge because its series are truncated at n^5<br>instead of n^4. However any of the other methods can leapfrog over [dk]<br>by using more terms in the series for (3) and (3') (and (1) and (1') in<br>the case of [se]). I posted series accurate to n^8 in<br><br>http://charles.karney.info/geographic/UTM-fi.txt<br><br>These can be directly used in [fi] and with minor modifications are<br>applicable to the other methods too. In the same file, I also give a<br>table of errors for the different order approximations (n^4, n^5, n^6,<br>n^7, and n^8). [The maxima code that I used to compute these series is<br>relatively short (40 lines or so), and I will post this at some point.]<br><br>We are then just left with relatively pedestrian numerical issues which<br>are important in the control of round-off:<br><br>* use of Horner representation of polynomials<br>* stable summation of trig series<br>* recognize that asin(tanh(Q)) is inaccurate if Q is large and use<br>the mathematically equivalent atan(sinh(Q)) instead<br>* and so on.<br><br>The use of these series expansions in (3) and (3') limits its<br>applicability. The transformation from conformal to rectifying latitude<br>when extended to the complex domain has a mild singularity a phi = 0,<br>lam = pi/2*(1-e) = 82.62d which is enough to make the series divergent<br>there. Indeed the zone of 1mm accuracy for the n^8 series is limited to<br>72 degrees from the central meridian.<br><br>Finally, all the methods which I have labeled by country codes derive<br>from the work of 2 Germans, Gauss and Krueger. The series (accurate to<br>n^4) needed for (3) and (3') are given by equations 26* (p. 18) and 41<br>(p. 21) in Krueger's 1912 paper<br><br>http://www.gfz-potsdam.de/bib/pub/digi/krueger2.pdf<br><br>So perhaps they all should be called [de].<br><br>--<span class="Apple-converted-space"> </span><br>Charles Karney <ckarney@sarnoff.com><br>Sarnoff Corporation, Princeton, NJ 08543-5300<br><br>URL: http://charles.karney.info<br>Tel: +1 609 734 2312<br>Fax: +1 609 734 2662<br>_______________________________________________<br>Proj mailing list<br>Proj@lists.maptools.org<br>http://lists.maptools.org/mailman/listinfo/proj<br><br></div></div></div></span></blockquote><br><div><br></div></div><div class="aol_ad_footer" id="uCF3455B1ED8347D98646985BDCF00227"><FONT style="color: black; font: normal 10pt ARIAL, SAN-SERIF;"><HR style="MARGIN-TOP: 10px">Looking for spoilers and reviews on the new TV season? <A title="http://television.aol.com/feature/fall_tv?ncid=aoletv00050000000037" href="http://television.aol.com/feature/fall_tv?ncid=aoletv00050000000037" target="_blank">Get AOL's ultimate guide to fall TV</A>.</FONT></div></body></html>