Coordinates Conversion

Note: This document is available in PDF.

Introduction

Coordinate conversion between Universal Transverse Mercator (UTM) and geographic coordinates is a very common problem when working with Geographic Information Systems.

This document is an English translation of a publication that is in Spanish [1] [2]. I use that method for coding my own program to convert coordinates.

There are some procedures to perform the conversion between geographic coordinates and UTM.

UTM Projection tables, from Servicio Geográfico del Ejército de España (SGE), Sección de Geodesia (1976): Proyección Universal Tranversa Mercator, SGE, Madrid. Vol. I: Sistemas conformes. Proyección U.T.M. Cuadrículas y Sistemas de referencia, (220 pp.); Vol. II: Tablas, (331 pp.).
US Army direct formulas for transformation, published in 1973 (USGS Bolletin No. 1532).
Coticchia-Surace formulas, from “Bolletino di Geodesia a Science Affini” No. 1.

This document is using the third method which is the easiest to code in a computer program. This procedure gets centimeter precision in the conversion if using enough decimal places, so floating point or double variables should be used.

Conversion from geographic coordinates to UTM

First we get a pair of geodesic geographic coordinates and the basic data of Hayford ellipsoid (as an example but you can use other ellipsoid), such as semi-major axis $(a)$ and semi-minor axis $(b)$ .

(1) $\begin{equation*}a = 6,378,388.0 \quad b = 6,356,911.94613\end{equation*}$

With this data we can start the operations.

Previous calculations

Ellipsoid geometry

We start to calculate the eccentricity, second eccentricity, polar radius of curvature and flattening. Eccentricity $(e)$ and second eccentricity $(e')$ are given by:

(2) $\begin{equation*}e = \frac{\sqrt{a^2 - b^2}}{a}\end{equation*}$

(3) $\begin{equation*}e' = \frac{\sqrt{a^2 - b^2}}{b}\end{equation*}$

The square of the second eccentricity is a very useful number in the operations so make sure to save it somewhere. The polar radius of curvature $(c)$ and flattening $(\alpha)$ are given by:

(4) $\begin{equation*}c = \frac{a^2}{b}\end{equation*}$

(5) $\begin{equation*}\alpha = \frac{a-b}{a}\end{equation*}$

Actually the flattening and eccentricity aren’t necessary for the Coticchia-Surace equations, but they are frequent parameters of the ellipsoid that sometimes are given as a pair of data $(a)$ and $(\alpha)$ ; and sometimes as $(a)$ and $(e)$ . Under those circumstances one can calculate $(b)$ if the corresponding equations are known.

Latitude and longitude

First, sexagesimal degrees (or degrees, minutes and seconds represented as $(DD)$ , $(MM)$ and $(SS)$ have to be converted to decimal notation.

(6) $\begin{equation*}dec = DD + \frac{MM}{60} + \frac{SS}{3600}\end{equation*}$

Once we have latitude $(\varphi)$ and longitude $(\lambda)$ in decimal notation, they have to be converted to radians. The operation is as follow

(7) $\begin{equation*}rad = \frac{dec \cdot \pi}{180}\end{equation*}$

Next, the longitude sign should be determined. This is as simple as

(8) $\begin{equation*}sign \left\{ \begin{array}{l}\text{If West, then negative (-)} \\ \text{If East, then positive (+)} \end{array} \right.\end{equation*}$

West and East are the longitude referenced from the Greenwich meridian.

Time zone

Once we have longitude and latitude ready, the time zone or UTM zone where the coordinates are located can be obtained.

(9) $\begin{equation*}zone = int \left[\frac{dec}{6} + 31\right]\end{equation*}$

Where $int$ represent only the integer part of the obtained value.

Next step is to get the central meridian of such time zone.

(10) $\begin{equation*}\lambda_0 = \left( zone \cdot 6 \right) - 183\end{equation*}$

Following we have to calculate the distance between the longitude of our original coordinates and the central meridian. Its important to notice that this operations have to be done in radians, so the central meridian just obtained have to be converted first. The original longitude is already in radians.

(11) $\begin{equation*}\Delta \lambda = \lambda - \lambda_0\end{equation*}$

Coticchia–Surace equations

In this section the Coticchia-Surace equations for the direct problem of converting geodesic geographic coordinates to UTM are presented.

Parameter calculation

A set of serial parameters should be calculated. There are a lot of equations but they are sequential and easy to code.

(12) $\begin{equation*}A = \cos \varphi \cdot \sin \Delta \lambda\end{equation*}$

(13) $\begin{equation*}\xi = \frac{1}{2} \cdot \ln \left[ \frac{1+A}{1-A} \right]\end{equation*}$

(14) $\begin{equation*}\eta = \arctan \left( \frac{\tan \varphi}{cos \Delta \lambda} \right) - \varphi\end{equation*}$

(15) $\begin{equation*}\nu = \frac{c}{(1 + e'^2 \cdot \cos^2 \varphi)^\frac{1}{2}} \cdot 0.9996\end{equation*}$

where $0.9996$ is the scale factor for the UTM projection.

(16) $\begin{equation*}\zeta = \frac{e'^2}{2} \cdot \xi^2 \cdot \cos^2 \varphi\end{equation*}$

(17) $\begin{equation*}A_1 = \sin (2 \cdot \varphi)\end{equation*}$

(18) $\begin{equation*}A_2 = A_1 \cdot \cos^2 \varphi\end{equation*}$

(19) $\begin{equation*}J_2 = \varphi + \frac{A_1}{2}\end{equation*}$

(20) $\begin{equation*}J_4 = \frac{3 \cdot J_2 + A_2}{4}\end{equation*}$

(21) $\begin{equation*}J_6 = \frac{5 \cdot J_4 + A_2 \cos^2 \varphi}{3}\end{equation*}$

(22) $\begin{equation*}\alpha = \frac{3}{4} \cdot e'^2\end{equation*}$

(23) $\begin{equation*}\beta = \frac{5}{3} \cdot \alpha^2\end{equation*}$

(24) $\begin{equation*}\gamma = \frac{35}{27} \cdot \alpha^3\end{equation*}$

(25) $\begin{equation*}B_\phi = 0.9996 \cdot c \cdot (\varphi - \alpha J_2 + \beta J_4 - \gamma J_6)\end{equation*}$

Final calculation of coordinates

Once we have all the parameters the final coordinates in UTM is obtained as follow

(26) $\begin{equation*}X = \xi \cdot \nu \cdot \left( 1 + \frac{\zeta}{3} \right) + 500,000\end{equation*}$

For the $Y$ solution is important to remember if the latitude of the geodesic geographic coordinates belongs to the southern hemisphere, we should add the value of $10,000,000$ to the result. In the case of northern hemisphere we have nothing to add.

(27) $\begin{equation*}Y = \eta \cdot \nu \cdot \left( 1 + \zeta \right) + B_\phi\end{equation*}$

Conversion from UTM to geographic coordinates

We are going to start by using a pair of UTM coordinates and their respective zone. Also we are going to need the basic data of Hayford ellipsoid (as an example but you can use another ellipsoid), such as semi-major axis $(a)$ and semi-minor axis $(b)$ .

(28) $\begin{equation*}a = 6,378,388.0 \quad b = 6,356,911.94613\end{equation*}$

With this data we can start the operations. The ellipsoid geometry calculations are exactly the same as the previous section, the parameter calculations are very similar.

Previous calculations

Ellipsoid geometry

We start to calculate the eccentricity, second eccentricity, polar radius of curvature and flattening. Eccentricity $(e)$ and second eccentricity $(e')$ are given by:

(29) $\begin{equation*}e = \frac{\sqrt{a^2 - b^2}}{a}\end{equation*}$

(30) $\begin{equation*}e' = \frac{\sqrt{a^2 - b^2}}{b}\end{equation*}$

The square of the second eccentricity is a very useful number in the operations so make sure to save it somewhere. The polar radius of curvature $(c)$ and flattening $(\alpha)$ are given by:

(31) $\begin{equation*}c = \frac{a^2}{b}\end{equation*}$

(32) $\begin{equation*}\alpha = \frac{a-b}{a}\end{equation*}$

Actually the flattening and eccentricity aren’t necessary for the Coticchia-Surace equations, but they are frequent parameters of the ellipsoid that sometimes are given as a pair of data $(a)$ and $(\alpha)$ ; and sometimes as $(a)) and \((e)$ . Under those circumstances one can calculate $(b)$ if the corresponding equations are known.

Previous calculations of (X, Y) coordinates

We start by modifying the X coordinate

(33) $\begin{equation*}X = X - 500,000\end{equation*}$

For the Y coordinate we should modify it only when it is on the southern hemisphere

(34) $\begin{equation*}sign \left\{ \begin{array}{l}\text{If Northern, then not modify Y} \\ \text{If Southern, then Y = Y - 10,000,000}\end{array} \right.\end{equation*}$

Central meridian of the time zone

The UTM zone of the X,Y coordinates to convert should be known as another parameter in the conversion. This operation is similar to the inverse conversion

(35) $\begin{equation*}\lambda_0 = \left( zone \cdot 6 \right) - 183\end{equation*}$

Coticchia-Surace equations

Parameter calculations

Most of the parameters are calculated in a similar way or even equal to the previous problem

(36) $\begin{equation*}\varphi' = \frac{Y}{6,366,197.724 \cdot 0.9996}\end{equation*}$

(37) $\begin{equation*}\nu = \frac{c}{(1 + e'^2 \cdot \cos^2 \varphi')^\frac{1}{2}} \cdot 0.9996\end{equation*}$

(38) $\begin{equation*}a = \frac{X}{\nu}\end{equation*}$

(39) $\begin{equation*}A_1 = \sin (2 \cdot \varphi')\end{equation*}$

(40) $\begin{equation*}A_2 = A_1 \cdot \cos^2 \varphi' \end{equation*}$

(41) $\begin{equation*}J_2 = \varphi' + \frac{A_1}{2}\end{equation*}$

(42) $\begin{equation*}J_4 = \frac{3 \cdot J_4 + A_2}{4}\end{equation*}$

(43) $\begin{equation*}J_6 = \frac{5 \cdot J_4 + A_2 \cdot \cos^2 \varphi'}{3}\end{equation*}$

(44) $\begin{equation*}\alpha = \frac{3}{4} \cdot e'^2\end{equation*}$

(45) $\begin{equation*}\beta = \frac{5}{4} \cdot \alpha^2\end{equation*}$

(46) $\begin{equation*}\gamma = \frac{35}{27} \cdot \alpha^3\end{equation*}$

(47) $\begin{equation*}B_\phi = 0.9996 \cdot c \cdot \left( \varphi' - \alpha \cdot J_2 + \beta \cdot J_4 - \gamma \cdot J_6 \right)\end{equation*}$

(48) $\begin{equation*}b = \frac{Y - B_\phi}{\nu}\end{equation*}$

(49) $\begin{equation*}\zeta = \frac{e'^2 \cdot a^2}{2} \cdot \cos^2 \varphi'\end{equation*}$

(50) $\begin{equation*}\xi = a \cdot \left[ 1 - \frac{\zeta}{3} \right]\end{equation*}$

(51) $\begin{equation*}\eta = b \cdot (1 - \zeta) + \varphi'\end{equation*}$

(52) $\begin{equation*}\sin h \xi = \frac{\exp^\xi - \exp^{-\xi}}{2}\end{equation*}$

(53) $\begin{equation*}\Delta \lambda = \arctan \frac{\sin h \xi}{\cos \eta}\end{equation*}$

(54) $\begin{equation*}\tau = \arctan (\cos \Delta \lambda \cdot \tan \eta)\end{equation*}$

Final calculation of coordinates

The longitude composition is very simple. We have to take care of performing this operation in degrees expressed in decimal notation, so $\Delta \lambda$ has to be divided by $\pi$ and multiplied by $180$ . $\lambda_0$ is already in decimal degrees

(55) $\begin{equation*}\lambda = \Delta \lambda + \lambda_0\end{equation*}$

The latitude calculation is a little more complicated:

(56) $\begin{equation*}\begin{array}\varphi = \varphi' + \left[ 1 + e'^2 \cdot \cos^2 \varphi' -\frac{3}{2} \cdot e'^2 \cdot \sin \varphi' \cdot \cos \varphi' \cdot (\tau - \varphi') \right] \cdot (\tau - \varphi')\end{array}\end{equation*}$

Now we have to convert the latitude from radians to degrees expressed in decimal notation

(57) $\begin{equation*}dec = \frac{\varphi (\text{in radians)}}{\pi} \cdot 180\end{equation*}$

Once we have latitude and longitude in degrees expressed in decimal notation, we have to convert the result to sexagesimal notation in degrees, minutes and seconds $(DD, MM, SS)$ . For the degrees we operate both latitude and longitude as follow

(58) $\begin{equation*}DD = int[dec]\end{equation*}$

where $int$ is the integer part of the value. Then minutes

(59) $\begin{equation*}MM = int[(dec - DD) \cdot 60]\end{equation*}$

again $int$ is the integer part of the value. Finally, seconds

(60) $\begin{equation*}SS = ((dec - DD) \cdot 60 - MM) \cdot 60\end{equation*}$

If longitude result with negative sign then the value its at west of Greenwich meridian.

References

Published: 5/1/2018, Minor fixes: 23/8/2024