Converting between 3D Cartesian and ellipsoidal latitude, longitude and height coordinates

All the formula here are coded into a spreadsheet available from the OS web site

Converting latitude, longitude and ellipsoid height to 3D Cartesian coordinates

These formulae are used to convert the format of coordinates between 3D Earth Centred Earth Fixed (ECEF) Cartesian coordinates and latitude, longitude and ellipsoidal height. This is NOT a datum transformation.

Values are required for the following ellipsoid constants: the semi-major axis length a and eccentricity squared . The latter can be calculated from a and b or the flattening f by

$$e^2 = \frac{a^2 - b^2}{a^2}= 2f -f^2$$

The Cartesian coordinates x y and z of a point are obtained from the latitude φ, longitude λ and ellipsoid height H by

$$\nu = \frac{a}{\sqrt{1-e^2sin^2\phi}}\\ x = (v+H)cos\phi cos\lambda\\ y = (v+H)cos\phi sin\lambda \\ z = ((1-e^2)v+H)\sin\phi$$

$$x = (v+H)cos\phi\lambda\\$$

$$y = (v+H)cos\phi\lambda\\$$

$$z = ((1-e^2)v+H\sin\phi\\$$

Here’s a worked example using the GRS80 ellipsoid. Intermediate values are shown here to 10 decimal places. Compute all values using double-precision arithmetic.

Latitude, φ 53° 36′ 43.1653″ N Longitude, λ 001° 39′ 51.9920″ W Ellipsoidal height, H 299.800 m

e² 6.6943800355E-03 ν 6.3920173768E+06

x 3790644.900 m y -110149.210 m z 5111482.970 m

Converting 3D Cartesian coordinates to latitude, longitude and ellipsoid height

Again, we need the defining constants of the ellipsoid. Longitude λ is easily computed from Cartesian coordinates, remembering to be careful about the quadrant of the resulting angle:

$$\lambda = arctan (\frac{y}{x})$$

The latitude φ is obtained by an iterative procedure. The initial value of latitude is given by

$\lambda = arctan (\frac{z}{p(1-e^2})$

where

$$p = \sqrt{x^2 + Y^2}$$

φ is iteratively improved by repeatedly computing ν from equation:

$$\nu = \frac{a}{\sqrt{1-e^2sin^2\phi}}\\$$

(using the latest value of φ) and then a new value for φ by

$$\phi = arctan[\frac{z+e^2vsin\phi}{p}]$$

until the change between two successive values of φ is smaller than the precision to which you want to calculate the latitude. Ellipsoid height H is then given by:

$$H = \frac{p}{cos\phi} - v$$

Here’s a worked example using the GRS80 ellipsoid. Intermediate values are shown here to 10 decimal places. Compute all values using double-precision arithmetic.

x 3790644.900 m y -110149.210 m z 5111482.970 m

e² 6.6943800355E-03 Initial φ 9.3570590125E-01 Initial ν 6.3920173799E+06 φ#2 9.3570575065E-01 ν#2 6.3920173768E+06 initial φ – φ#2 1.5059657843E-07 φ#3 9.3570575030E-01 ν#3 6.3920173768E+06 φ#2 – φ#3 3.5634395434E-10 φ#4 9.3570575029E-01 ν#4 6.3920173768E+06

φ#3 – φ#4 8.4310336490E-13 φ#5 9.3570575029E-01 ν#5 6.3920173768E+06 φ#4 – φ#5 1.9984014443E-15 φ#6 9.3570575029E-01 ν#6 6.3920173768E+06 φ#5 – φ#6 0.0000000000E+00

Latitude, φ 53° 36′ 43.1653″ N Longitude, λ 001° 39′ 51.9920″ W Ellipsoidal height, H 299.800 m