User:Mark viking/Geographic coordinate conversion

In geodesy, conversion among different geographic coordinate systems is made necessary by the different geographic coordinate systems in use across the world and over time. Coordinate conversion comprises a number of different types of conversion: format change, conversion of coordinate systems, or transformation to different geodetic datums. Geographic coordinate conversion has applications in cartography, surveying, and geographic information systems.

In geodesy, geographic coordinate conversion is defined as translation among different coordinate formats or map projections all referenced to the same geodetic datum. A geographic coordinate transformation is a translation among different geodetic datums. Both conversion and translations will be considered in this article.

Coordinate format conversion

Informally, specifying a geographic location usually means giving the location's latitude and longitude. The numerical values for latitude and longitude can occur in a number of different formats:

• degrees minutes seconds: 40° 26′ 46″ N 79° 58′ 56″ W
• degrees decimal minutes: 40° 26.767′ N 79° 58.933′ W
• decimal degrees: 40.446° N 79.982° W

There are 60 minutes in a degree and 60 seconds in a minute, so to convert from degrees minutes seconds to decimal degrees,

${\rm {{decimal\ degrees}={\rm {{degrees}+{\rm {{minutes}/60+{\rm {{seconds}/3600}}}}}}}}$.

To convert back from decimal degrees to degrees minutes seconds,

${\begin{aligned}{\rm {degrees}}&=\lfloor {\rm {{decimal\ degrees}\rfloor }}\\{\rm {minutes}}&=\lfloor 60*({\rm {{decimal\ degrees}-{\rm {{degrees})\rfloor }}}}\\{\rm {seconds}}&=\lfloor 3600*({\rm {{decimal\ degrees}-{\rm {{degrees}-{\rm {{minutes}/60)\rfloor }}}}}}\\\end{aligned}}$,

where the notation $\lfloor x\rfloor$ means take the integer part of of $x$.

Specifying a geographic location

Main article: Geographic coordinate system

In geodesy to unambiguously specify a location in three dimensions requires in addition to latitude and longitude, a terrestrial reference system, and the height or altitude. To specify a location on a two-dimensional map requires in addition the specification of a map projection.

Terrestrial reference system

Main article: Geodetic datum

A terrestrial reference system, also called a geodetic datum, is a mathematical model of the earth's shape that defines the relation between a location and the coordinates that describe it. Modern datum models represent the earth as reference ellipsoids. Datums may be global, meaning that they represent the whole earth, or they may be local, meaning that they represent a best-fit ellipsoid to only a portion of the earth. Examples of global datums include World Geodetic System (WGS 84), the default datum used for Global Positioning System and the International Terrestrial Reference Frame (ITRF) used for estimating continental drift and crustal deformation. Given a location, the datum provides the latitude $\phi$ and longitude $\lambda$.

Height

The height is specified relative a baseline. Common baselines include

• The surface of the datum ellipsoid, resulting in an ellipsoidal height
• The mean sea level as described by the gravity geoid, yielding the orthometric height
• A vertical datum, yielding a dynamic height relative to a known reference height.

Along with the latitude $\phi$ and longitude $\lambda$, the height $h$ provides the three-dimensional geodetic coordinates for a location.

An alternative coordinate system to geodetic coordinates is the Earth-centered, Earth fixed (ECEF) coordinate system. This is a three dimensional Cartesian coordinate system with $X$, $Y$ and $Z$ coordinates. In ECEF, height is implicit in the coordinates. In general the particular ECEF coordinate axes chosen are related to geodetic datum specified; thus ECEF systems of different datums are not equivalent.

Map projection

Main article: Map projection

To establish the position of a geographic location on a map, a map projection is used to convert geodetic coordinates to two-dimensional coordinates on a map; it projects the ellipsoidal shape of a datum onto a flat surface of a map. The datum along with a map projection establishes a grid system for plotting locations. Common map projections in current use include the Universal Transverse Mercator (UTM), the Military grid reference system (MGRS), the United States National Grid (USNG), the Global Area Reference System (GARS) and the World Geographic Reference System (GEOREF). Coordinates on a map are usually in terms northing N and easting E offsets relative to a specified origin.

Coordinate system conversion

A coordinate system conversion is conversion from one coordinate system to another, with both corrdinate systems based on the same geodetic datum. Common conversion tasks include conversion between geodetic and ECEF coordinates and conversion from one type of map projection to another.

From geodetic to ECEF coordinates

The length PQ is called Normal ($\,N(\phi )$). The length IQ is equal to $\,e^{2}N(\phi )$. R = $\,(X,Y,Z)$.

Geodetic coordinates (latitude $\ \phi$, longitude $\ \lambda$, height $h$) can be converted into ECEF coordinates using the following formulae:

${\begin{aligned}X&=\left(N(\phi )+h\right)\cos {\phi }\cos {\lambda }\\Y&=\left(N(\phi )+h\right)\cos {\phi }\sin {\lambda }\\Z&=\left(N(\phi )(1-e^{2})+h\right)\sin {\phi }\end{aligned}}$

where

$N(\phi )={\frac {a}{\sqrt {1-e^{2}\sin ^{2}\phi }}},$

$a$ and $e$ are the semi-major axis and the first numerical eccentricity of the ellipsoid respectively.
$\,N(\phi )$ is called the Normal and is the distance from the surface to the Z-axis along the ellipsoid normal (see "Radius of curvature on the Earth"). The following equation holds:

${\frac {p}{\cos \phi }}-{\frac {Z}{\sin \phi }}-e^{2}N(\phi )=0,$

where $p={\sqrt {X^{2}+Y^{2}}}$.

The orthogonality of the coordinates is confirmed via differentiation:

${\begin{aligned}&{\big (}dX,\,dY,\,dZ{\big )}\\[6pt]=&{\big (}-\sin \phi \cos \lambda ,\,-\sin \phi \sin \lambda ,\,\cos \phi {\big )}\left(M(\phi )+h\right)\,d\phi \\[6pt]&{}+{\big (}-\sin \lambda ,\,\cos \lambda ,\,0{\big )}\left(N(\phi )+h\right)\cos \phi \,d\lambda \\[6pt]&{}+{\big (}\cos \lambda \cos \phi ,\,\cos \phi \sin \lambda ,\,\sin \phi {\big )}\,dh,\end{aligned}}$

where

$M(\phi )={\frac {a(1-e^{2})}{\left(1-e^{2}\sin ^{2}\phi \right)^{3/2}}}$

(see also "Meridian arc on the ellipsoid").

From ECEF to geodetic coordinates

The conversion of ECEF coordinates to geodetic coordinates (such WGS84) is a much harder problem, except for longitude, $\,\lambda$.

There exist two kinds of methods in order to solve the equation.

Newton–Raphson method

The following Bowring's irrational geodetic-latitude equation is efficient to be solved by Newton–Raphson iteration method:

$\kappa -1-{\frac {e^{2}a\kappa }{\sqrt {p^{2}+(1-e^{2})z^{2}\kappa ^{2}}}}=0,$

where $\kappa ={\frac {p}{z}}\tan \phi$. The height is calculated as follows:

$h=e^{-2}(\kappa ^{-1}-{\kappa _{0}}^{-1}){\sqrt {p^{2}+z^{2}\kappa ^{2}}},$

$\kappa _{0}=\left(1-e^{2}\right)^{-1}.$

The iteration can be transformed into the following calculation:

$\kappa _{i+1}={\frac {c_{i}+\left(1-e^{2}\right)z^{2}\kappa _{i}^{3}}{c_{i}-p^{2}}}=1+{\frac {p^{2}+\left(1-e^{2}\right)z^{2}\kappa _{i}^{3}}{c_{i}-p^{2}}},$

where $c_{i}={\frac {\left(p^{2}+\left(1-e^{2}\right)z^{2}\kappa _{i}^{2}\right)^{3/2}}{ae^{2}}}.$

$\,\kappa _{0}$ is a good starter for the iteration when $h\approx 0$. Bowring showed that the single iteration produces the sufficiently accurate solution. He used extra trigonometric functions in his original formulation.

Ferrari's solution

The quartic equation for this transformation can be solved by Ferrari's solution to yield:

${\begin{aligned}\zeta &=(1-e^{2})z^{2}/a^{2},\\[6pt]\rho &=(p^{2}/a^{2}+\zeta -e^{4})/6,\\[6pt]s&=e^{4}\zeta p^{2}/(4a^{2}),\\[6pt]t&={\sqrt[{3}]{\rho ^{3}+s+{\sqrt {s(s+2\rho ^{3})}}}},\\[6pt]u&=\rho +t+\rho ^{2}/t,\\[6pt]v&={\sqrt {u^{2}+e^{4}\zeta }},\\[6pt]w&=e^{2}(u+v-\zeta )/(2v),\\[6pt]\kappa &=1+e^{2}({\sqrt {u+v+w^{2}}}+w)/(u+v).\end{aligned}}$

Conversion across map projections

Conversion of coordinates and map positions among different map projections may be accomplished either through direct translation formulas from one projection to another, or by first converting from a projection $A$ to an intermediate coordinate system, such as ECEF, then converting from ECEF to projection $B$. The formulas involved can be complex and in some cases, such as in the ECEF to geodetic conversion above, the conversion has no closed-form solution and approximate methods must be used. References such as the DMA Technical Manual 8358.1 and the USGS paper Map Projections: A Working Manual contain formulas for conversion of map projections. It is more common to use computer programs to perform coordinate conversion tasks, such as the DoD and NGA supported GEOTRANS program.

Datum transformations

Transformations among datums can be accomplished in a number of ways. There are transformations that directly convert geodetic coordinates from one datum to another. There are more indirect transforms that convert from geodetic coordinates to ECEF coordinates, transform the ECEF coordinates from one datum to the another, then transform ECEF coordinated of the new datum back to geodetic coordinates. There are also grid-based transformations that directly transform from one (datum, map projection) pair to another (datum, map projection) pair.

Molodensky transformation

The Molodensky transformation converts directly between geodetic coordinate systems of different datums without the intermediate step of converting to geocentric ECEF coordinates. It requires the three shifts between the datum centers and the differences between the ellipsoid semi-major axes and flattening parameters.

The Molodensky transform is used by the National Geospatial-Intelligence Agency (NGA) in their standard TR8350.2 and the NGA supported GEOTRANS program. The Molodensky method used to popular before the advent of modern computers and the method is part of many geodetic programs.

Helmert transformation

Main article: Helmert transformation

Use of the Helmert transform in the transformation from geodetic coordinates of datum $A$ to geodetic coordinates of datum $B$ occurs in the context of a three step process:

1. Convert from geodetic coordinates to ECEF coordinates for datum $A$
2. Apply the Helmert transform, with the appropriate $A\to B$ transform parameters, to transform from datum $A$ ECEF coordinates to datum $B$ ECEF coordinates
3. Convert from ECEF coordinates to geodetic coordinates for datum $B$

In terms of ECEF XYZ vectors, the Helmert transform has the form

${\begin{bmatrix}X_{B}\\Y_{B}\\Z_{B}\end{bmatrix}}={\begin{bmatrix}c_{x}\\c_{y}\\c_{z}\end{bmatrix}}+(1+s\times 10^{-6})\cdot {\begin{bmatrix}1&-r_{z}&r_{y}\\r_{z}&1&-r_{x}\\-r_{y}&r_{x}&1\end{bmatrix}}\cdot {\begin{bmatrix}X_{A}\\Y_{A}\\Z_{A}\end{bmatrix}}.$

The Helmert transform is a seven-parameter transform with three translation (shift) parameters $c_{x},c_{y},c_{z}$, three rotation parameters $r_{x},r_{y},r_{z}$ and one scaling (dilation) parameter $s$. The Helmert transform is an approximate method that is accurate when the transformation parameters are small relative to the magnitudes of theECEF vectors. Under these conditions, the transform is considered reversible.

A fourteen-parameter Helmert transform, with linear time dependence for each parameter, can be used to capture the time evolution of geographic coordinates dues to geomorphic processes, such as continental drift. This has been incorporated into software, such as the Horizontal Time Dependent Positioning (HTDP) tool from the U.S. NGS.

Molodensky-Badekas transformation

To eliminate the coupling between the rotations and translations of the Helmert transformation, three additional parameters can be introduced to give a new XYZ center of rotation closer to coordinates being transformed. This ten-parameter model is called the Molodensky-Badekas transformation and should not be confused with the more basic Molodensky transformation.

Like the Helmert transform, using the Molodensky-Badekas transformation is a three step process:

1. Convert from geodetic coordinates to ECEF coordinates for datum $A$
2. Apply the Molodensky-Badekas transform, with the appropriate $A\to B$ transform parameters, to transform from datum $A$ ECEF coordinates to datum $B$ ECEF coordinates
3. Convert from ECEF coordinates to geodetic coordinates for datum $B$

The transform has the form

${\begin{bmatrix}X_{B}\\Y_{B}\\Z_{B}\end{bmatrix}}={\begin{bmatrix}X_{A}\\Y_{A}\\Z_{A}\end{bmatrix}}+{\begin{bmatrix}\Delta X_{A}\\\Delta Y_{A}\\\Delta Z_{A}\end{bmatrix}}+{\begin{bmatrix}1&-r_{z}&r_{y}\\r_{z}&1&-r_{x}\\-r_{y}&r_{x}&1\end{bmatrix}}\cdot {\begin{bmatrix}X_{A}-X_{A}^{0}\\Y_{A}-Y_{A}^{0}\\Z_{A}-Z_{A}^{0}\end{bmatrix}}+\Delta S{\begin{bmatrix}X_{A}-X_{A}^{0}\\Y_{A}-Y_{A}^{0}\\Z_{A}-Z_{A}^{0}\end{bmatrix}}.$

where $(X_{A}^{0},Y_{A}^{0},Z_{A}^{0})$ is the origin for the rotation and scaling transforms and $\Delta S$ is the scaling factor.

The Molodensky-Badekas transformation is used to transform local geodetic datums to a global geodetic datum, such as WGS 84. Unlike the Helmert transformation, the Molodensky-Badekas transformation is not reversible due to the rotational origin being associated with the original datum.

Multiple regression equations

Datum transformations through the use of empirical multiple regression methods were create to achieve higher accuracy results over small regions than the standard Molodensky transformations. MRE transforms are used to transform local datums over continent-sized or smaller regions to global datums, such as WGS 84. The standard NIMA TM 8350.2, Appendix D, lists MRE transforms from several local datums to WGS 84, with accuracies of generally 2 meters.

The MREs are a direct transformation of geodetic coordinates. Geodetic coordinates $\phi _{B},\lambda _{B},h_{B}$ in the new datum $B$ are modeled as low-order polynomials of the geodetic coordinates $\phi _{A},\lambda _{A},h_{A}$ of the original datum $A$. Given a number of $(A,B)$ coordinate pairs for landmarks in both datums, multiple regression methods are used to fit the parameters of these polynomials.

Grid-based method

Grid-based transformations directly convert map coordinates from one (map-projection, geodetic datum) pair to map coordinates of another (map-projection, geodetic datum) pair. An example is the NADCON method for transforming from the North American Datum (NAD) 1927 to the NAD 1983 datum. The High Accuracy Reference Network (HARN), a high accuracy version of the NADCON transforms, have an accuracy of approximately 5 centimeters. The National Transformation version 2 (NTv2) is a Canadian version of NADCON for transforming between NAD 1927 and NAD 1983. HARNs are also known as NAD 83/91 and High Precision Grid Networks (HPGN). Like the multiple regression equation transform , grid-based methods use a low-order interpolation for convert map coordinates, but in two dimensions instead of three. The NOAA provides a software tool (as part of the NGS Geodetic Toolkit) for performing NADCON transformations.