Editing Geographic coordinate conversion (section)

==Coordinate system conversion==
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.

=== From geodetic to ECEF coordinates ===
[[Image:Geodetic latitude and the length of Normal.svg|thumb|The length PQ, called the ''prime vertical radius'', is <math>N(\phi)</math>. The length IQ is equal to <math>\, e^2 N(\phi) </math>. <math>R = (X,\, Y,\, Z)</math>.]]

[[Geodetic coordinates]] (latitude <math>\ \phi</math>, longitude <math>\ \lambda</math>, height <math>h</math>) can be converted into [[ECEF]] coordinates using the following equation:<ref name="gps-chap10">{{cite book|title=GPS - theory and practice|author1=B. Hofmann-Wellenhof |author2=H. Lichtenegger |author3=J. Collins |isbn=3-211-82839-7|page=282|others=Section 10.2.1|year=1997 }}</ref>

: <math>\begin{align}
      X & = \left( N(\phi) + h\right)\cos{\phi}\cos{\lambda} \\
      Y & = \left( N(\phi) + h\right)\cos{\phi}\sin{\lambda} \\
      Z & = \left( \frac{b^2}{a^2} N(\phi) + h\right)\sin{\phi} \\
        & = \left( (1 - e^2)       N(\phi) + h\right)\sin{\phi} \\
        & = \left( (1 - f)^2       N(\phi) + h\right)\sin{\phi}
\end{align}</math>

where
: <math>
  N(\phi) = \frac{a^2}{\sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi }}
          = \frac{a}{\sqrt{1 - e^2\sin^2\phi}} = \frac{a}{\sqrt{1 - \frac{e^2}{1 + \cot^2 \phi}}},
</math>

and <math>a</math> and <math>b</math> are the equatorial radius ([[semi-major axis]]) and the polar radius ([[semi-minor axis]]), respectively. <math>e^2 = 1 - \frac{b^2}{a^2}</math> is the square of the first numerical eccentricity of the ellipsoid. <math>f = 1 - \frac{b}{a}</math> is the flattening of the ellipsoid. The ''[[prime vertical radius of curvature]]'' <math>\, N(\phi) </math> is the distance from the surface to the Z-axis along the ellipsoid normal.

====Properties====
The following condition holds for the longitude in the same way as in the geocentric coordinates system:
:<math>\frac{X}{\cos\lambda} - \frac{Y}{\sin\lambda} = 0.</math>

And the following holds for the latitude:
:<math>\frac{p}{\cos\phi} - \frac{Z}{\sin\phi} - e^2 N(\phi) = 0,</math>

where <math>p = \sqrt{X^2 + Y^2}</math>, as the parameter <math>h</math> is eliminated by subtracting
:<math>\frac{p}{\cos\phi} = N + h</math>

and
:<math>\frac{Z}{\sin\phi} = \frac{b^2}{a^2}N + h.</math>

The following holds furthermore, derived from dividing above equations:
:<math>\frac{Z}{p} \cot \phi = 1 - \frac{e^2 N}{N + h}.</math>

====Orthogonality====
The [[Orthogonal coordinates|orthogonality]] of the coordinates is confirmed via differentiation:
:<math>\begin{align}
  \begin{pmatrix} dX \\ dY \\ dZ \end{pmatrix} &= 
    \begin{pmatrix}
      -\sin\lambda & -\sin\phi \cos\lambda & \cos\phi \cos\lambda \\
       \cos\lambda & -\sin\phi \sin\lambda & \cos\phi \sin\lambda \\
                  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\lambda \\ d\phi \\ dh \end{pmatrix},
\end{align}</math>
<!--
: <math>\begin{align}
      & \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{align}</math>
-->

where
:<math>
 M(\phi) = \frac{a\left(1 - e^2\right)}{\left(1 - e^2 \sin^2 \phi\right)^\frac{3}{2}} = N(\phi) \frac{1 - e^2}{1 - e^2\sin^2\phi}
</math>

(see also "[[Meridian arc#Definition|Meridian arc on the ellipsoid]]").
<!--
The infinitesimal length caused by latitude and longitude is calculated as follows (see also "[[Meridian arc#Meridian distance on the ellipsoid|Meridian arc on the ellipsoid]]"):
: <math>
 ds^2 = \left(\frac{a\left(1 - e^2\right)}{\left(1 - e^2 \sin^2\phi\right)^\frac{3}{2}} + h\right)^2 d\phi^2 + \left(\frac{a}{\sqrt{1 - e^2 \sin^2\phi}} + h\right)^2 \cos^2\phi\, d\lambda^2 .
</math>
-->

=== From ECEF to geodetic coordinates ===
====Conversion for the longitude====
The conversion of ECEF coordinates to longitude is:
: <math>\lambda = \operatorname{atan2}(Y,X)</math>. 
where [[atan2]] is the quadrant-resolving arc-tangent function.
The geocentric longitude and geodetic longitude have the same value; this is true for Earth and other similar shaped planets because they have a large amount of rotational symmetry around their spin axis (see [[triaxial ellipsoidal longitude]] for a generalization).

==== Simple iterative conversion for latitude and height ====
The conversion for the latitude and height involves a circular relationship involving ''N'', which is a function of latitude:
:<math>\frac{Z}{p} \cot \phi = 1 - \frac{e^2 N}{N + h}</math>,
:<math>h=\frac{p}{\cos\phi} - N</math>.
It can be solved iteratively,<ref name=osgb>A guide to coordinate systems in Great Britain. This is available as a pdf document at
{{cite web|url=http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents |title=ordnancesurvey.co.uk |access-date=2012-01-11 |url-status=dead |archive-url=https://web.archive.org/web/20120211075826/http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/ |archive-date=2012-02-11 }} Appendices B1, B2</ref><ref name=osborne>Osborne, P (2008). [http://mercator.myzen.co.uk/mercator.pdf The Mercator Projections] {{webarchive|url=https://web.archive.org/web/20120118224152/http://mercator.myzen.co.uk/mercator.pdf |date=2012-01-18 }} Section 5.4</ref> for example, starting with a first guess ''h''≈0 then updating ''N''.
More elaborate methods are shown below. 
The procedure is, however, sensitive to small accuracy due to <math>N</math> and <math>h</math> being maybe 10{{sup|6}} apart.<ref>[https://web.archive.org/web/20080920155754/http://www.ferris.edu/faculty/burtchr/papers/cartesian_to_geodetic.pdf R. Burtch, A Comparison of Methods Used in Rectangular to Geodetic Coordinate Transformations.]</ref><ref>{{cite journal |last1=Featherstone |first1=W. E. |last2=Claessens |first2=S. J. |title=Closed-Form Transformation between Geodetic and Ellipsoidal Coordinates |journal=Stud. Geophys. Geod. |volume=52 |issue=1 |pages=1–18 |year=2008 |doi=10.1007/s11200-008-0002-6 |bibcode=2008StGG...52....1F |hdl=20.500.11937/11589 |s2cid=59401014 |hdl-access=free }}</ref>
<!-- There are several methods that solve the equation; two are shown. -->

==== Newton–Raphson method ====
The following Bowring's irrational geodetic-latitude equation,<ref>{{cite journal |last=Bowring |first=B. R. |title=Transformation from Spatial to Geographical Coordinates |journal=Surv. Rev. |volume=23 |issue=181 |pages=323–327 |year=1976 |doi=10.1179/003962676791280626 }}</ref> derived simply from the above properties, is efficient to be solved by [[Newton–Raphson]] iteration method:<ref>{{cite journal |last=Fukushima |first=T. |title=Fast Transform from Geocentric to Geodetic Coordinates |journal=J. Geod. |volume=73 |issue=11 |pages=603–610 |year=1999 |doi=10.1007/s001900050271 |bibcode=1999JGeod..73..603F |s2cid=121816294 }} (Appendix B)</ref><ref>{{cite book|first1=J. J. |title=Proceedings of the IEEE 1997 National Aerospace and Electronics Conference. NAECON 1997|volume=2|pages=646–650|last1=Sudano|doi=10.1109/NAECON.1997.622711|chapter=An exact conversion from an earth-centered coordinate system to latitude, longitude and altitude|year=1997|isbn=0-7803-3725-5|s2cid=111028929 }}</ref>
: <math>\kappa - 1 - \frac{e^2 a\kappa}{\sqrt{p^2 + \left(1 - e^2\right) Z^2 \kappa^2}} = 0,</math>

where <math>\kappa = \frac{p}{Z} \tan \phi</math> and <math>p = \sqrt{X^2 + Y^2}</math> as before. The height is calculated as:
: <math>\begin{align}
         h &= e^{-2} \left(\kappa^{-1} - {\kappa_0}^{-1}\right) \sqrt{p^2 + Z^2 \kappa^2}, \\
  \kappa_0 &\triangleq \left(1 - e^2\right)^{-1}.
\end{align}</math>

The iteration can be transformed into the following calculation:
: <math>\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},</math>

where <math>c_i = \frac{\left(p^2 + \left(1 - e^2\right) Z^2 \kappa_i ^2\right)^\frac{3}{2}}{ae^2} .</math>

The constant <math>\,\kappa_0</math> is a good starter value for the iteration when <math>h \approx 0</math>. Bowring showed that the single iteration produces a sufficiently accurate solution. He used extra trigonometric functions in his original formulation.
<!--
: <math>\kappa \approx \kappa_1 = \left(c + \frac{z^2}{1 - e^2 }\right)/\left(c - \left(1 - e^2\right)\left(x^2 + y^2\right)\right),</math>

where
: <math>c = \frac{\left(\left(1 - e^2\right)\left(x^2 + y^2\right) + z^2\right)^\frac{3}{2}}{ae^2 \sqrt{1 - e^2}}.</math> --><!--
For <math>h = 0</math>, <math>\kappa = \frac{1}{1 - e^2}</math>, which is a good starter for the iteration. Bowring showed that the single iteration produces the sufficiently accurate solution under the condition of <math>h \approx 0</math>.
-->

==== Ferrari's solution ====
The quartic equation of <math>\kappa</math>, derived from the above, <!--for this transformation--> can be solved by [[Quartic equation#Ferrari.27s solution|Ferrari's solution]]<ref>{{cite journal|first1=H. |last1=Vermeille, H.|title=Direct Transformation from Geocentric to Geodetic Coordinates|journal= J. Geod.|volume=76|number=8|pages=451–454
|year= 2002|doi=10.1007/s00190-002-0273-6|s2cid=120075409 }}</ref><ref>{{cite journal|first1=Laureano|last1=Gonzalez-Vega|first2=Irene|last2=PoloBlanco|title=A symbolic analysis of Vermeille and Borkowski polynomials for transforming 3D Cartesian to geodetic coordinates|journal=J. Geod.|volume=83|number=11|pages=1071–1081|doi=10.1007/s00190-009-0325-2|year=2009|bibcode=2009JGeod..83.1071G |s2cid=120864969 }}</ref> to yield:
: <math>
\begin{align}
   \zeta &= \left(1 - e^2\right)\frac{z^2}{a^2} ,\\[4pt]
    \rho &= \frac{1}{6}\left(\frac{p^2}{a^2} + \zeta - e^4\right) ,\\[4pt]
       s &= \frac{e^4 \zeta p^2}{4\rho^3 a^2} ,\\[4pt]
       t &= \sqrt[3]{1 + s + \sqrt{s(s + 2)}} ,\\[4pt]
       u &= \rho \left(t + 1 + \frac{1}{t}\right) ,\\[4pt]
       v &= \sqrt{u^2 + e^4 \zeta} ,\\[4pt]
       w &= e^2 \frac{u + v - \zeta}{2v} ,\\[4pt]
  \kappa &= 1 + e^2 \frac{\sqrt{u + v + w^2} + w}{u + v}.
\end{align}
</math>

===== The application of Ferrari's solution =====
A number of techniques and algorithms are available but the most accurate, according to Zhu,<ref>{{cite journal|first1=J.|last1=Zhu|title=Conversion of Earth-centered Earth-fixed coordinates to geodetic coordinates|journal=IEEE Transactions on Aerospace and Electronic Systems|volume=30|issue=3|year=1994|pages=957–961|doi=10.1109/7.303772|bibcode=1994ITAES..30..957Z }}</ref> is the following procedure established by Heikkinen,<ref>{{cite journal|first1=M.|last1=Heikkinen|title=Geschlossene formeln zur berechnung räumlicher geodätischer koordinaten aus rechtwinkligen koordinaten.|journal=Z. Vermess.|volume=107|year=1982|pages=207–211|language=de}}</ref> as cited by Zhu. This overlaps with above. It is assumed that geodetic parameters <math>\{a,\, b,\, e\}</math> are known

: <math>\begin{align}
        a &= 6378137.0 \text{     m. Earth Equatorial Radius} \\[3pt]
        b &= 6356752.3142 \text{      m. Earth Polar Radius} \\[3pt]
      e^2 &= \frac{a^2-b^2}{a^2} \\[3pt]
     e'^2 &= \frac{a^2 - b^2}{b^2} \\[3pt]
        p &= \sqrt{X^2 + Y^2} \\[3pt]
        F &= 54b^2 Z^2 \\[3pt]
        G &= p^2 + \left(1 - e^2\right)Z^2 - e^2\left(a^2 - b^2\right) \\[3pt]
        c &= \frac{e^4 Fp^2}{G^3} \\[3pt]
        s &= \sqrt[3]{1 + c + \sqrt{c^2 + 2c}} \\[3pt]
        k &= s + 1 + \frac{1}{s}\\[3pt]
        P &= \frac{F}{3 k^2 G^2} \\[3pt]
        Q &= \sqrt{1 + 2e^4 P} \\[3pt]
      r_0 &= \frac{-Pe^2 p}{1 + Q} + \sqrt{\frac{1}{2} a^2\left(1 + \frac{1}{Q}\right) - \frac{P\left(1 - e^2\right)Z^2}{Q(1 + Q)} - \frac{1}{2}Pp^2} \\[3pt]
        U &= \sqrt{\left(p - e^2 r_0\right)^2 + Z^2} \\[3pt]
        V &= \sqrt{\left(p - e^2 r_0\right)^2 + \left(1 - e^2\right)Z^2} \\[3pt]
      z_0 &= \frac{b^2 Z}{aV} \\[3pt]
        h &= U\left(1 - \frac{b^2}{aV}\right) \\[3pt]
     \phi &= \arctan\left[\frac{Z + e'^2 z_0}{p}\right] \\[3pt]
  \lambda &= \operatorname{arctan2}[Y,\, X]
\end{align}</math>

Note: [[atan2|arctan2]][Y, X] is the four-quadrant inverse tangent function.

==== Power series ====
For small {{math|e<sup>2</sup>}} the power series
:<math>\kappa = \sum_{i\ge 0} \alpha_i e^{2i}</math>

starts with
:<math>\begin{align}
  \alpha_0 &= 1; \\
  \alpha_1 &= \frac{a}{\sqrt{Z^2 + p^2}}; \\
  \alpha_2 &= \frac{aZ^2\sqrt{Z^2 + p^2} + 2a^2 p^2}{2\left(Z^2 + p^2\right)^2}.
\end{align}</math>

=== Geodetic to/from ENU coordinates ===
To convert from geodetic coordinates to [[local tangent plane]] ([[Axes conventions#Ground reference frames: ENU and NED|ENU]]) coordinates is a two-stage process:

# Convert geodetic coordinates to ECEF coordinates
# Convert ECEF coordinates to local ENU coordinates

==== From ECEF to ENU ====

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 <math>\left\{X_r,\, Y_r,\, Z_r\right\}</math> and an aircraft at <math>\left\{X_p,\, Y_p,\, Z_p\right\}</math>, then the vector pointing from the radar to the aircraft in the ENU frame is

: <math>
  \begin{bmatrix}x \\ y \\ z\end{bmatrix} =
  \begin{bmatrix}
    -\sin\lambda_r &            \cos\lambda_r           &         0 \\
    -\sin\phi_r\cos\lambda_r & -\sin\phi_r\sin\lambda_r & \cos\phi_r \\
     \cos\phi_r\cos\lambda_r &  \cos\phi_r\sin\lambda_r & \sin\phi_r
  \end{bmatrix}
  \begin{bmatrix}
    X_p - X_r \\
    Y_p - Y_r \\
    Z_p - Z_r
  \end{bmatrix}
</math>

Note: <math>\ \phi</math> is the ''[[geodetic latitude]]''; the ''[[geocentric latitude]]'' is inappropriate for representing [[vertical direction]] for the local tangent plane and must be [[geocentric latitude|converted]] if necessary.

==== From ENU to ECEF ====

This is just the inversion of the ECEF to ENU transformation so

: <math>
  \begin{bmatrix}X_p \\ Y_p \\ Z_p\end{bmatrix} =
  \begin{bmatrix}
    -\sin\lambda_r &  -\sin\phi_r\cos\lambda_r & \cos\phi_r\cos\lambda_r \\
     \cos\lambda_r &  -\sin\phi_r\sin\lambda_r & \cos\phi_r\sin\lambda_r \\
               0 &   \cos\phi_r            & \sin\phi_r
  \end{bmatrix}
  \begin{bmatrix}x \\ y \\ z\end{bmatrix} +
  \begin{bmatrix}X_r \\ Y_r \\ Z_r\end{bmatrix}
</math>

=== Conversion across map projections ===

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 <math>A</math> to an intermediate coordinate system, such as ECEF, then converting from ECEF to projection <math>B</math>. 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''<ref name=TM8358.2>{{cite web|title=TM8358.2: The Universal Grids: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)|url=http://earth-info.nga.mil/GandG/publications/tm8358.2/TM8358_2.pdf|publisher=National Geospatial-Intelligence Agency|access-date=4 March 2014|archive-date=3 March 2020|archive-url=https://web.archive.org/web/20200303184711/https://earth-info.nga.mil/GandG/publications/tm8358.2/TM8358_2.pdf|url-status=live}}</ref> and the USGS paper ''Map Projections: A Working Manual''<ref name=Snyder1987>{{cite book|last=Snyder|first=John P.|title=Map Projections: A Working Manual|year=1987|publisher=USGS Professional Paper: 1395|url=https://pubs.er.usgs.gov/publication/pp1395|access-date=2017-08-28|archive-date=2011-05-17|archive-url=https://web.archive.org/web/20110517082057/http://pubs.er.usgs.gov/publication/pp1395|url-status=live}}</ref> contain formulas for conversion of map projections. It is  common to use computer programs to perform coordinate conversion tasks, such as with the DoD and NGA supported GEOTRANS program.<ref name=GEOTRANS_NGA>{{cite web|title=MSP GEOTRANS 3.3 (Geographic Translator)|url=http://earth-info.nga.mil/GandG/geotrans/|publisher=NGA: Coordinate Systems Analysis Branch|access-date=4 March 2014|archive-date=15 March 2014|archive-url=https://web.archive.org/web/20140315075748/http://earth-info.nga.mil/GandG/geotrans/|url-status=live}}</ref>