[OSRS-PROJ] Transverse Mercator revisited

[Gerald I. Evenden]

Maybe a month or so ago there was a claim about ESRI developing a Transverse
Mercator projection from an "obscure" NOAA publication.  Part of the claim was
that the projection was good to 90 degrees from the CM.  I more or less said
"bovine excrement" and off list got into a discussion with someone who
supported the claim.

Subsequently I managed to get a reference to the NOAA article and locate
the author who was kind enough to send me a copy.  After messing around with
the article I wrote the following critque and, as you will note, have not got
a projection that will extend to 90 degrees longitude---nor 90 NS for that
matter.

The reason I post it here is that if anyone has a copy of ESRI's software and
could puch in a few bench mark conversions I would be most grateful.
Does ESRI's software go to 90 degree longitude??

A copy of this critique has also been sent to Jeff Dozier.  Sorry, but some
comments relate to material of the article and are not meaningful without
having a copy.

-----------------------------------------------------------------------------------------

A programmers review of "Improved Algorithm for
Calculation of UTM and Geodetic Coordinates" by
J. Dozier.

FIRST TRY.

First attempt to encode the C code was by typing and
simultaneous conversion of 70's C to C99 level and
thus using the current C language complex arithmetic,
subroutines as well as prototyping capabilities.
This allowed removal of the "8. complex" section
as well as several complex transcendental functions.
Similarly, ".im" and ".re" arithmetic was change to
straight forward complex arithmetic notation and
operations.

Original code sections were placed in individual
files and each contains an '#include "dozier.h"' that
gives full prototyping for all procedures as well as
some constant definitions.

Basic projection entry was changed to the 'gkfor'
procedure modified so that the input arguments are
only eccentricity, geodetic latitude and longitude
and with x and y output.  Note that x-y in the current
version are reversed when compared to most western
notation.  The ellipse major axis, a, was omitted
(as per usage in the PROJ.4 system) and assumed to
be 1.  Conversion of geographic latitude to isometic
latitude was made inside 'gkfor.'

In addition, the failure tests and calls to "abend"
were changed to a long-jump return to the main entry
into the projection so as to allow a graceful return
with error code rather than aborting the program.
'gkfor' now has a zero functional value if the projection
was a success or non-zero value for failure along with
the returned x-y values set to HUGE_VAL.

An error in coding is believed to have occur in the original
article code in file tmw.c, procedure tmfd, pg. 16.
Original code reads:

ms = snw;
ms.re *= m;
ms.im *= m;

where m was assigned k_tran*k_tran a few lines earlier.
k_tran should be eccentricity at that point, thus m now
becomes eccentricity squared.  In the statement following
the above the operation (1-ms)*(1+ms) is made which is
effectively 1-ms^2 or, in this case (1-(e^4*snw^2)).
I believe k_tran should have been in the aforementioned
operations rather than m.  The interesting thing about
this error is that it was involved with the computation
of the derivitive in the Newton-Raphson iteration.
My experience is that an error in determining this
deriviative may only slightly impair the covergence
of the algorithm as long as it is "nearly" correct,
but accuracy of the result is not affected.

Currently, the software fails when longitude exceeds
about 82 degrees as the jseries function fails to
converge.  This is apparent from looking at the
defining equations which involve sin(z) where z
is complex.  Unlike real sine functions whose
results stay within +-1, complex sine values may
exceed 1 and thus the increasing size of the
multiple complex sine may negate the decreasing
size of the series coefficient.  My question is:
is the value of the sn w function finite as the
imaginary part of the argument approaches pi/2?
The jseries function is modified to have an iteration
check and abort functioning if excessive iterations
have been made.

It should be noted that when eccentricity is set to
zero (spherical earth) the software performs flawlessly
except for failure to compute near lat,lon=90,0 and 0,90.

Note that in both the spherical equations and the
elliptical series equations a latitude of 90
degrees is computable.  but in this system the latitude
is converted to and from isometric latitude which is
indeterminate at 90 degrees.  Note that this causes
the procedure geodat to immediately fail because it
sets a range to the Brent algorithm of 0 to 90
degrees which causes function 'zerobr' to fail almost
immediately.  At the moment, I have been concentrating
on the forward projection and have not made patches
to the inverse problem---in this case, limiting the
polar range to something less than pi/2.

Comparisons with lproj's tmerc show agreement to
within about 1e-10 over most of the +-3.5 lon
and latitude (<90) range. Because the certified
accuracy of tmerc is unknown this is not the
greatest test, but along the central meridian
(lon=0) the northing is known to about 1e-14
and checkable against extended precision numerical
integration values BUT this version of the Dozier
software only gets 1e-10 along the central meridian.
Roughly 1mm on the ground.

TRY AGAIN:

For the second attempt, the article was scanned and
OCR was used to make a transcription of the text to
computer readable text.  OCR is certainly not fool
proof and additional careful editing of the text was
required as well as exhaustive cross checking of the
new code with the source.  In this operation care
was made to ensure that the computer readable text
matched as closely as possible to the original article.

Again, each section was put into an individual file
and it was necessary to add a header inclusion to
"lcomplex.h" that contains:

typedef struct {
    double re, im;
} complex;

Attempts were made to remain as close to 70's K&R
C as possible.

Some other code additions were required. The function
"cdiv" was not declare in gk.c and tmwi.c and "csub"
was not declared in tmw.c.  Compilation error occurred
in two cases when not corrected and segmentation fault
occurred in the other.

The first, "recoded" version works better at the moment
so I suspended work on the scanned code.

CONCLUSIONS

I am having considerable problem in assembling a
functional program from the material in this article.
I have been able to verify some components of the
code by use of the GNU Scientific Library and the
Numerical Recipes library but many parts do not have
available counterparts---at least none that I can
locate.

A minor comment is that the parts of the paper
related to UTM are irrelevant.  In reality, UTM is
"defined" by the series solution to the transverse
mercator projection.  Even if the series solution is
not precise, it still is the UTM "solution."
This method cannot replace the current UTM solution
nor other established cadastral coordinate systems
but it is of academic interest in that it could give a
better idea of the accuracy of the series solutions.

Some people have expressed an in interest in using
this method when the longitude range extends well
beyond the 3.5 degree limits of current transverse
mercator series solutions.  But as extent expands
the scale shrinks and resolutions needs are lessened.

A side note here: part of my personal interest was
the result of a discussion where someone claimed that
the transverse meridian projection  produced finite
results with longitude values 90 degrees from the
central meridian.  The Dozier paper also makes this
claim on page 1: "The method can be used for all
longitudes up to and including 90 deg from the central
meridian."  This claim was counter intuitive to me and
thus made me look into this assertion.  So far, indications
that I observe are that as I approach 90 degrees
longitude the value of the easting starts to rapidly
increase much in the same manner as in the case where
polar latitudes create ever larger northings as one
approaches the Pole in normal mercator projection.

Lastly, this paper speaks volumes for including
numeric examples when publishing articles related to
practical solutions and where there is no machine
readable source code.  Even when machine readable
code is available, test results should be available
to verify the installation.
----------------------------------------
PROJ.4 Discussion List
See http://www.remotesensing.org/proj for subscription, unsubscription
and other information.