[The following is a highly technical note (i.e, odds are you don't want to read it at all), which I'm writing in the form of a blog post mainly because I'm too lazy to start another Web page—but it might get copied elsewhere eventually (if someone wants to dump it on Wikipedia, feel free to).]
I've long been intrigued by the oft-used (especially in conjunction with GPS units) UTM coordinates: explanations as to how UTM coordinates are defined exactly is very hard to find (the spherical case is simple enough, but the ellipsoidal one is a much tougher nut), and although there are many online tools (such as this one) to convert from latitude+longitude to UTM and back, they use black box formulæ, converging power series, and looking at the source will give you very little insight on what is going on. So here is an attempt at a mathematically precise definition (it took me a whole day of angry formula-crunching before I came up with something entirely correct, so I won't spare the details).
First of all, what is a transverse mercator projection? It can be defined by starting from a central meridian (UTM uses 60 possible central meridians, one for each UTM zone) and constructing the mathematically unique projection from a spheroid (= rotation of an ellipse about one of its axes) to a plane which is conformal (= preserves angles) and maps the central meridian to a straight line with constant scale. So here are the defining features of Universal Transverse Mercator:
- It divides the Earth in 60 longitude zones, each 6° wide: zone 1 ranges from 180° to 174°W, zone 2 from 174°W to 168°W and so on through zone 30 from 6°W to 0°, zone 31 from 0° to 6°E, up to zone 60 from 174°E to 180°. Each zone's central meridian is halfway between the limiting longitudes (e.g., zone 33's central meridian is 15°E). UTM coordinates are given relative to a longitude zone and a hemisphere (North or South).
- In each zone, UTM is strictly conformal (= preserves angles). In practice, this means that it does not distort shapes anywhere, or that the UTM grid is a square grid everywhere, and it appears square when displayed on any conformal map projection. (Note that this does not imply that the lines of the grid are aligned north-south and east-west! But they are always perpendicular.)
- Each zone's central meridian (which receives an
eastingcoordinate of 500000, see below) is mapped to a vertical line along which the vertical (northing) coordinate is simply proportional to distance (on the ellipsoid) along that meridian; however,- the scale along the central meridian is not 1:1, it is 9996:10000. In other words, two points lying at 10km's distance from one another on the central meridian will have vertical (
northing) coordinates differing by only 9996. The reason for this is that the map scale increases on either side of the meridian (as is unavoidable when mapping a curve surface) and the 0.04% decrease at the center is chosen to make the average scale on a 6°-wide zone roughly 1:1.- The two UTM coordinates, easting (given first) and northing (given second) are orthogonal and measured in meters (after the transverse mercator projection and the 9996:10000 scale are applied). Easting is measured horizontally, with the central meridian having value 500000 (in practice, it can take values from 166021 to 833979 at the equator, the range being more narrow at higher latitudes). Northing is measured vertically, with the equator having value 0 in the northern hemisphere and 10000000 in the southern hemisphere; in the northern hemisphere, it ranges from 0 at the equator to 9328094 at 84°N on the central meridian, or even 9997965 if used all the way to the pole (but this normally isn't supposed to happen: beyond 84°N, and beyond 80°S Universal Polar Stereographic coordinates should be used instead of UTM).
- UTM is (nowadays) almost exclusively used with the WGS84 ellipsoid: this has a semimajor axis of a = 6378137m (exactly) and a flattening of f = 1 − b/a = 1/298.257223563 (exactly).
This is a complete definition, however it isn't a very usable one
because it lacks a description of how to use the conformal
part
of the definition to actually convert geodetic coordinates to
Universal Transverse Mercator. If we were to use a spherical Earth
model, computations would be easy enough:
(Transverse Mercator computations, spherical Earth:)
- Start with latitude χ and longitude ϖ, and let ϖ0 be the central meridian's longitude;
- if we rotate the Earth to make the central meridian become the new (transverse) equator and the old equator become the new (transverse) reference meridian, then simple spheric trigonometry shows that the given point's transverse (i.e., new, i.e., rotated) longitude is ϖ† = arctan[tan(χ)/cos(ϖ−ϖ0)] and its transverse latitude is χ† = arcsin[sin(ϖ−ϖ0)·cos(χ)];
- now if we apply the formula for Mercator projection to these transverse coordinates (namely: longitude in radians gives one coordinate and the other is given by log(tan(π/4+χ†/2)) where χ† is longitude), we get the transverse Mercator x = log(tan(π/4+χ†/2)) = arctanh[sin(ϖ−ϖ0)·cos(χ)] and y = ϖ† = arctan[tan(χ)/cos(ϖ−ϖ0)] (in radians).
- It is then a simple matter of scaling (multiply x and y by the Earth's radius and by the 0.9996 scale factor) and translating (add 500000 meters to the easting coordinate and possibly 10000000 to the northing if in the southern hemisphere) to obtain UTM coordinates.
Here is a numerical example: consider the point with latitude χ=45° and longitude ϖ=0° and refer it to UTM zone 31's central meridian of ϖ0=3°. We find x = −0.037024017523 and y = 0.786083865778. If we convert these in meters by taking the Earth's radius to be such that the meridian's length is the same as for the WGS84 model (10001965.7293m), and taking into accound the 0.9996 scale factor, then we get an easting coordinate of 264345.75067 and a northing of 5003346.90008. As we shall see below, these coordinates differ quite sensibly (by about 16km) from the correct values obtained for an ellipsoidal Earth.
For future reference, the complex number z = y − i·x will be called the complex latitude (this is my terminology: but it makes sense since when the point is on the central meridian ϖ=ϖ0 then z is real and equals the latitude χ; so in the spherical Earth case, computing the transverse Mercator coordinates is basically just the same as computing this complex latitude).
Now in the case of a spheroidal Earth, things are much more
complicated. For one thing, whereas longitude remains a simple
concept (since we are assuming the ellipsoid in question has equal
equatorial axes), the simple spherical latitude gives rise to
half
a dozen different concepts on a spheroid. The usual
(geodetic) latitude, φ, is the angle which a
line perpendicular to the ellipsoid at a given point forms with the
equatorial plane. This is what the unqualified word latitude
almost always means. But there are two other latitudes which are of
importance:
- Conformal latitude χ is so defined that
it is an increasing function of (geodetic) latitude which coincides
with it at the equator (χ=φ=0°) and the
poles (χ=φ=±90°), and such that mapping
the ellipsoid to a sphere (of arbitrary radius) by taking a given
point on the ellipsoid to the point on the sphere having the same
longitude and latitude χ (the conformal latitude
we are trying to define) gives a conformal map from the
ellipsoid to the sphere.
(To restate this in a less convoluted way: if you parametrize the spheroid by conformal latitude and longitude, you can pretend that it's a sphere when defining conformal projections.)
A rather simple computation shows that χ is related to φ by dχ/cos(χ) = dφ/cos(φ) · (cos²(ε)/[cos²(φ) + cos²(ε)·sin²(φ)]) with ε the angular eccentricity (e=sin(ε)), and on integrating this gives the horrendous formula which is displayed in the Wikipedia article. - Rectifying latitude μ is so defined that it is an increasing function of (geodetic) latitude which coincides with it at the equator (χ=φ=0°) and the poles (χ=φ=±90°) and which varies linearly with distance measured along the meridian. This is simple to define but it has no closed-form expression except involving incomplete elliptic integrals.
Obviously conformal latitude is of importance to us because we are
trying to define a conformal map. If we simply apply the
spherical-Earth formulæ with conformal latitude instead of (geodetic)
latitude we get a conformal map but it is not the transverse Mercator
we are looking for because coordinates along the central meridian will
not be proportional to distance as they are supposed to be: for that,
rectifying latitude needs to play a role. In fact, rectifying
latitude is exactly what we need (up to proportionality and a trivial
additive constant in the southern hemisphere) to get the vertical
(northing
) coordinate along the central meridian. So somehow
we need to use both the conformal and the rectifying latitudes.
How?
The trick is to consider the function M which takes a conformal latitude χ to the corresponding rectifying latitude μ=M(χ). That function is very difficult to write down explicitly, because it involves taking the reciprocal of a horrendous formula and then yet applying an (incomplete) elliptic integral. Nevertheless, it exists conceptually, and it is an analytic real function: so it has a locally unique holomorphic extension to a complex neighborhood of its real domain of definition; the reason we care is that holomorphic complex functions are conformal mappings. In practice, M can be approached by various power series techniques.
So the procedure to compute UTM coordinates on an elliptic Earth is this:
(Transverse Mercator computations, ellipsoidal Earth:)
- Start with latitude φ and longitude ϖ: compute the conformal latitude χ corresponding to φ (there is a closed-form expression for this, so we are happy);
- apply the spherical-Earth computations, using χ as latitude: let z = y − i·x be what we have previously called the complex latitude (here the complex conformal latitude, of course): at this point, we have defined a conformal map, but the central meridian is parametrized by conformal latitude, not rectifying latitude as we would like it to be;
- compute z♠ = M(z), where M is the (complex holomorphic extension of the) function taking a conformal latitude to the corresponding rectifying latitude: then z♠ = y♠ − i·x♠ where x♠ and y♠ are the desired UTM coordinates
{proof: since M is holomorphic, we have defined a conformal map, and since M takes a real conformal latitude to the corresponding rectifying latitude, it is what it should be on the real axis, but since real values suffice to identify an analytic function we are done}.Here is a numerical example: consider the point with latitude φ=45° and longitude ϖ=0° and refer it to UTM zone 31's central meridian of ϖ0=3°. We find a conformal latitude of χ=44°48′27.662601″. Now applying the spherical-Earth transverse Mercator formulæ, we find x = −0.037148195538 and y = 0.782727307059. Then M(y−i·x) = y♠−i·x♠ where x♠ = −0.037148414843 and y♠ = 0.783567347070. If we convert these in meters by taking into account the meridian's length for the WGS84 model (10001965.7293m), and taking into accound the 0.9996 scale factor, then we get an easting coordinate of 263553.97390 and a northing of 4987329.50469. This is the correct conversion: the point (45°N,0°) has UTM coordinates (263553.974,4987329.505) in zone 31 North.
Note that if we had naïvely applied the spherical-Earth transverse Mercator formulæ to the conformal latitude we would have found the incorrect value (263555.370,4981982.732), off by 5km; if we had applied the same formulæ to the rectifying latitude immediately, we would have found (263752.383,4987314.788), still off by 200m: so even for GPS-precision computations we can't simply forget about the complexity.
If we are aiming at a precision of slightly better than 1 metre, then the following numerical series expansions for M(z) (with the constants of WGS84) can be used: the first is around z=0 (to be used for latitudes less than 40°, say) and the second is around z=π/2 (to be used for latitudes above 40°):
- M(z) ≈ 1.00167851426z − 0.00112513485z3 + 0.00022996604z5 − 0.00002381665z7 + 1.76580·10−6z9 + O(z11)
- M(z) ≈ π/2 + 0.99832757260(z-π/2) + 0.00110890291(z-π/2)3 − 0.00021697949(z-π/2)5 + 0.00001886787(z-π/2)7 − 6.6472·10−7(z-π/2)9 + O((z-π/2)11)
I'm not sure how well these power series expansions compare with other classical approaches to computing UTM coordinates (e.g., pages 60–64 of this book).