[Proj] Transverse Mercator algorithm

[Charles Karney]

Gerald I. Evenden wrote:
 > I appreciate your approach to error diagnosis in the comparison of the
 > ?tmercs however I must take a more brute force approach and merely
 > look at primitive results.  This is due to the fact that I am *not* a
 > theoretician but merely one looking for the bottom line.

It's probably good to have an idea of what the error might look like as
a function of position.  The error will approximately be the first
dropped term in the series expansions.  In the case of the 4-term JHS
154 expansion this is

     ~ e^10 sin(10 xi) cosh(10 eta)
     ~ e^10 sin(10 phi) cosh(10 x/a)

or some oscillating function of y multiplying an exponential function of
x.  We are, of course, primarily interested in the potentially large
variation with x.

 > I will get a 3D plot made sometime today.

There's a good chance that this will look messy because of the
oscillatory behavior of the error with latitude and the fact that
round-off errors dominate for small x.

My methodology is:

Generate a set of random latitudes (lat0), longitudes (lon0), exact
eastings (x0), and exact northings (y0).  I posted a set of such data
(1/4 million points) at

     http://charles.karney.info/geographic/TMcoords.dat.gz

For each geographic point (lat0,lon0), compute

     (x,y) = F(lat0,lon0)
     (lat,lon) = F^-1(x,y)

where F is the numerical transverse Mercator projection being
evaluated.  Define

     err = max( hypot(x-x0, y-y0),
                1e7/90 * hypot(lat-lat0, cos(lat0)*(lon-lon0)) )

     mu = asin( sin(lon0) * cos(lat0) )

The second term in err approximately converts the lat/lon discrepancy
into a distance.  mu is the angular distance away from the central
meridian and I use this as an ersatz x of the purpose of error analysis.

You can get a feel for the behavior of the error by plotting err against
mu.  E.g., in Matlab notation:

     plot(mu, log10(err), 'x');

The tabulated error data I sent out yesterday was computed with, e.g.,

     max(err(x0<5e5 & y0<96e5))
     max(err(mu<25))

and so on.  The second statement takes the maximum of the error for all
mu < 25 and this is useful because (from the plot) err is strongly
increasing with mu.

 > Of course, the above has one critical assumption we must not forget:
 > the maxima procedure is correct and god-like.

Quite so.  The good news is that when the procedure fails, it fails
badly.

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