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 is composed of a number of different types of conversion: format change of geographic coordinates, conversion of coordinate systems, or transformation to different geodetic datums. Geographic coordinate conversion has applications in cartography, surveying, navigation 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.[1] A geographic coordinate transformation is a translation among different geodetic datums. Both geographic coordinate conversion and transformation will be considered in this article.
This article assumes readers are already familiar with the content in the articles geographic coordinate system and geodetic datum.
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 units or formats:[2]
There are 60 minutes in a degree and 60 seconds in a minute. Therefore, to convert from a degrees minutes seconds format to a decimal degrees format, one may use the formula
\rm{decimal degrees}=\rm{degrees}+
\rm{minutes | |
To convert back from decimal degree format to degrees minutes seconds format,
\begin{align} \rm{absDegrees}&=|\rm{decimal degrees}|\\ \rm{floorAbsDegrees}&=\lfloor\rm{absDegrees}\rfloor\\ \rm{degrees}&=sgn(\rm{decimal degrees}) x \rm{floorAbsDegrees}\\ \rm{minutes}&=\lfloor60 x (\rm{absDegrees}-\rm{floorAbsDegrees})\rfloor\\ \rm{seconds}&=3600 x (\rm{absDegrees}-\rm{floorAbsDegrees})-60 x \rm{minutes}\\ \end{align}
\rm{absDegrees}
\rm{floorAbsDegrees}
A coordinate system conversion is a conversion from one coordinate system to another, with both coordinate systems based on the same geodetic datum. Common conversion tasks include conversion between geodetic and earth-centered, earth-fixed (ECEF) coordinates and conversion from one type of map projection to another.
Geodetic coordinates (latitude
\phi
λ
h
\begin{align} X&=\left(N(\phi)+h\right)\cos{\phi}\cos{λ}\\ Y&=\left(N(\phi)+h\right)\cos{\phi}\sin{λ}\\ Z&=\left(
b2 | |
a2 |
N(\phi)+h\right)\sin{\phi}\\ &=\left((1-e2)N(\phi)+h\right)\sin{\phi}\\ &=\left((1-f)2N(\phi)+h\right)\sin{\phi} \end{align}
where
N(\phi)=
a2 | |
\sqrt{a2\cos2\phi+b2\sin2\phi |
and
a
b
e2=1-
b2 | |
a2 |
f=1-
b | |
a |
N(\phi)
The following condition holds for the longitude in the same way as in the geocentric coordinates system:
X | |
\cosλ |
-
Y | |
\sinλ |
=0.
And the following holds for the latitude:
p | |
\cos\phi |
-
Z | |
\sin\phi |
-e2N(\phi)=0,
where
p=\sqrt{X2+Y2}
h
p | |
\cos\phi |
=N+h
and
Z | |
\sin\phi |
=
b2 | |
a2 |
N+h.
The following holds furthermore, derived from dividing above equations:
Z | |
p |
\cot\phi=1-
e2N | |
N+h |
.
The orthogonality of the coordinates is confirmed via differentiation:
\begin{align} \begin{pmatrix}dX\ dY\ dZ\end{pmatrix}&=\begin{pmatrix} -\sinλ&-\sin\phi\cosλ&\cos\phi\cosλ\\ \cosλ&-\sin\phi\sinλ&\cos\phi\sinλ\\ 0&\cos\phi&\sin\phi\\ \end{pmatrix}\begin{pmatrix}dE\ dN\ dU\end{pmatrix},\\[3pt] \begin{pmatrix}dE\ dN\ dU\end{pmatrix}&= \begin{pmatrix} \left(N(\phi)+h\right)\cos\phi&0&0\\ 0&M(\phi)+h&0\\ 0&0&1\\ \end{pmatrix} \begin{pmatrix}dλ\ d\phi\ dh\end{pmatrix}, \end{align}
where
M(\phi)=
a\left(1-e2\right) | |||||||||
|
=N(\phi)
1-e2 | |
1-e2\sin2\phi |
(see also "Meridian arc on the ellipsoid").
The conversion of ECEF coordinates to longitude is:
λ=\operatorname{atan2}(Y,X)
The conversion for the latitude and height involves a circular relationship involving N, which is a function of latitude:
Z | |
p |
\cot\phi=1-
e2N | |
N+h |
h= | p |
\cos\phi |
-N
N
h
The following Bowring's irrational geodetic-latitude equation,[8] derived simply from the above properties, is efficient to be solved by Newton–Raphson iteration method:[9] [10]
\kappa-1-
e2a\kappa | |
\sqrt{p2+\left(1-e2\right)Z2\kappa2 |
where
\kappa=
p | |
Z |
\tan\phi
p=\sqrt{X2+Y2}
\begin{align} h&=e-2\left(\kappa-1-
-1 | |
{\kappa | |
0} |
\right)\sqrt{p2+Z2\kappa2},\\ \kappa0&\triangleq\left(1-e2\right)-1. \end{align}
The iteration can be transformed into the following calculation:
\kappai+1=
| |||||||||||||
ci-p2 |
=1+
| |||||||||||||
ci-p2 |
,
where
ci=
| ||||||||||||||||||||
ae2 |
.
The constant
\kappa0
h ≈ 0
The quartic equation of
\kappa
\begin{align} \zeta&=\left(1-
| ||||
e |
,\\[4pt] \rho&=
1 | \left( | |
6 |
p2 | |
a2 |
+\zeta-e4\right),\\[4pt] s&=
e4\zetap2 | |
4\rho3a2 |
,\\[4pt] t&=\sqrt[3]{1+s+\sqrt{s(s+2)}},\\[4pt] u&=\rho\left(t+1+
1 | |
t |
\right),\\[4pt] v&=\sqrt{u2+e4\zeta},\\[4pt] w&=e2
u+v-\zeta | |
2v |
,\\[4pt] \kappa&=1+e2
\sqrt{u+v+w2 | |
+ |
w}{u+v}. \end{align}
A number of techniques and algorithms are available but the most accurate, according to Zhu,[13] is the following procedure established by Heikkinen,[14] as cited by Zhu. This overlaps with above. It is assumed that geodetic parameters
\{a,b,e\}
\begin{align} a&=6378137.0m.EarthEquatorialRadius\\[3pt] b&=6356752.3142m.EarthPolarRadius\\[3pt] e2&=
a2-b2 | |
a2 |
\\[3pt] e'2&=
a2-b2 | |
b2 |
\\[3pt] p&=\sqrt{X2+Y2}\\[3pt] F&=54b2Z2\\[3pt] G&=p2+\left(1-e2\right)Z2-e2\left(a2-b2\right)\\[3pt] c&=
e4Fp2 | |
G3 |
\\[3pt] s&=\sqrt[3]{1+c+\sqrt{c2+2c}}\\[3pt] k&=s+1+
1 | |
s |
\\[3pt] P&=
F | |
3k2G2 |
\\[3pt] Q&=\sqrt{1+2e4P}\\[3pt] r0&=
-Pe2p | |
1+Q |
+\sqrt{
1 | |
2 |
a2\left(1+
1 | |
Q |
\right)-
P\left(1-e2\right)Z2 | |
Q(1+Q) |
-
1 | |
2 |
Pp2}\\[3pt] U&=\sqrt{\left(p-e2
2 | |
r | |
0\right) |
+Z2}\\[3pt] V&=\sqrt{\left(p-e2
2 | |
r | |
0\right) |
+\left(1-e2\right)Z2}\\[3pt] z0&=
b2Z | |
aV |
\\[3pt] h&=U\left(1-
b2 | |
aV |
\right)\\[3pt] \phi&=\arctan\left[
Z+e'2z0 | |
p |
\right]\\[3pt] λ&=\operatorname{arctan2}[Y,X] \end{align}
Note: arctan2[Y, X] is the four-quadrant inverse tangent function.
For small the power series
\kappa=\sumi\ge\alphaie2i
starts with
\begin{align} \alpha0&=1;\\ \alpha1&=
a | |
\sqrt{Z2+p2 |
To convert from geodetic coordinates to local tangent plane (ENU) coordinates is a two-stage process:
To transform from ECEF coordinates to the local coordinates we need a local reference point. Typically, this might be the location of a radar. If a radar is located at
\left\{Xr,Yr,Zr\right\}
\left\{Xp,Yp,Zp\right\}
\begin{bmatrix}x\ y\ z\end{bmatrix}= \begin{bmatrix} -\sinλr&\cosλr&0\\ -\sin\phir\cosλr&-\sin\phir\sinλr&\cos\phir\\ \cos\phir\cosλr&\cos\phir\sinλr&\sin\phir \end{bmatrix} \begin{bmatrix} Xp-Xr\\ Yp-Yr\\ Zp-Zr \end{bmatrix}
Note:
\phi
This is just the inversion of the ECEF to ENU transformation so
\begin{bmatrix}Xp\ Yp\ Zp\end{bmatrix}= \begin{bmatrix} -\sinλr&-\sin\phir\cosλr&\cos\phir\cosλr\\ \cosλr&-\sin\phir\sinλr&\cos\phir\sinλr\\ 0&\cos\phir&\sin\phir \end{bmatrix} \begin{bmatrix}x\ y\ z\end{bmatrix}+ \begin{bmatrix}Xr\ Yr\ Zr\end{bmatrix}
Conversion of coordinates and map positions among different map projections reference to the same datum may be accomplished either through direct translation formulas from one projection to another, or by first converting from a projection
A
B
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 another, then transform ECEF coordinates 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.
See main article: Helmert transformation.
Use of the Helmert transform in the transformation from geodetic coordinates of datum
A
B
A
A\toB
A
B
B
In terms of ECEF XYZ vectors, the Helmert transform has the form (position vector transformation convention and very small rotation angles simplification)[18]
\begin{bmatrix}XB\ YB\ ZB\end{bmatrix}= \begin{bmatrix}cx\ cy\ cz\end{bmatrix}+\left(1+s x 10-6\right) \begin{bmatrix} 1&-rz&ry\\ rz&1&-rx\\ -ry&rx&1 \end{bmatrix}\begin{bmatrix}XA\ YA\ ZA\end{bmatrix}.
The Helmert transform is a seven-parameter transform with three translation (shift) parameters
cx,cy,cz
rx,ry,rz
s
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[20] and earthquakes.[21] This has been incorporated into software, such as the Horizontal Time Dependent Positioning (HTDP) tool from the U.S. NGS.[22]
To eliminate the coupling between the rotations and translations of the Helmert transform, 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 transform.
Like the Helmert transform, using the Molodensky-Badekas transform is a three-step process:
A
A\toB
A
B
B
The transform has the form[23]
\begin{bmatrix}XB\ YB\ ZB\end{bmatrix}= \begin{bmatrix}XA\ YA\ ZA\end{bmatrix}+ \begin{bmatrix}\DeltaXA\ \DeltaYA\ \DeltaZA\end{bmatrix}+ \begin{bmatrix} 1&-rz&ry\\ rz&1&-rx\\ -ry&rx&1 \end{bmatrix} \begin{bmatrix}XA-
0 | |
X | |
A |
\ YA-
0 | |
Y | |
A |
\ ZA-
0 | |
Z | |
A |
\end{bmatrix}+ \DeltaS\begin{bmatrix}XA-
0 | |
X | |
A |
\ YA-
0 | |
Y | |
A |
\ ZA-
0 | |
Z | |
A |
\end{bmatrix}.
where
0 | |
\left(X | |
A, |
0 | |
Y | |
A, |
0 | |
Z | |
A\right) |
\DeltaS
The Molodensky-Badekas transform is used to transform local geodetic datums to a global geodetic datum, such as WGS 84. Unlike the Helmert transform, the Molodensky-Badekas transform is not reversible due to the rotational origin being associated with the original datum.
The Molodensky transformation converts directly between geodetic coordinate systems of different datums without the intermediate step of converting to geocentric coordinates (ECEF).[24] It requires the three shifts between the datum centers and the differences between the reference 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.[25] The Molodensky method was popular before the advent of modern computers and the method is part of many geodetic programs.
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.[26] 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).[27] Subsequently, Australia and New Zealand adopted the NTv2 format to create grid-based methods for transforming among their own local datums.
Like the multiple regression equation transform, grid-based methods use a low-order interpolation method for converting 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.[28] [29]
Datum transformations through the use of empirical multiple regression methods were created to achieve higher accuracy results over small geographic 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.[30] The standard NIMA TM 8350.2, Appendix D,[31] lists MRE transforms from several local datums to WGS 84, with accuracies of about 2 meters.[32]
The MREs are a direct transformation of geodetic coordinates with no intermediate ECEF step. Geodetic coordinates
\phiB,λB,hB
B
\phiA,λA,hA
A
\phiB
\Delta\phi=a0+a1U+a2V+a3U2+a4UV+a5V2+ …
where
ai,
\begin{align} U&=K(\phiA-\phim)\\ V&=K(λA-λm)\\ \end{align}
K,
\phim,λm,
A.
with similar equations for
\Deltaλ
\Deltah
(A,B)