[GIS] Converting from ECEF to geodetic coordinates

coordinate systemmathematicsvector

I am trying to convert an ECEF coordinate to a geodetic coordinate using this procedure, but I don't fully understand the process. It is prefaced by stating the geodetic parameters {a,b,e,e'} are assumed to be known, but does not state what they reference to.

I assume that a and b are the equatorial and polar semi-axes, and that e is just Euler's number. Am I correct in this? And what is e' supposed to be?

Best Answer

You can use this function to perform the conversion of ECEF coordinates to Geodetic coordinates.

The function is implemented in python but can be easily written in java.

Example:

import math

def xyz2llh(x,y,z):
    '''
    Function to convert xyz ECEF to llh
    convert cartesian coordinate into geographic coordinate
    ellipsoid definition: WGS84
      a= 6,378,137m
      f= 1/298.257

    Input
      x: coordinate X meters
      y: coordinate y meters
      z: coordinate z meters
    Output
      lat: latitude rad
      lon: longitude rad
      h: height meters
    '''
    # --- WGS84 constants
    a = 6378137.0
    f = 1.0 / 298.257223563
    # --- derived constants
    b = a - f*a
    e = math.sqrt(math.pow(a,2.0)-math.pow(b,2.0))/a
    clambda = math.atan2(y,x)
    p = math.sqrt(pow(x,2.0)+pow(y,2))
    h_old = 0.0
    # first guess with h=0 meters
    theta = math.atan2(z,p*(1.0-math.pow(e,2.0)))
    cs = math.cos(theta)
    sn = math.sin(theta)
    N = math.pow(a,2.0)/math.sqrt(math.pow(a*cs,2.0)+math.pow(b*sn,2.0))
    h = p/cs - N
    while abs(h-h_old) > 1.0e-6:
        h_old = h
        theta = math.atan2(z,p*(1.0-math.pow(e,2.0)*N/(N+h)))
        cs = math.cos(theta)
        sn = math.sin(theta)
        N = math.pow(a,2.0)/math.sqrt(math.pow(a*cs,2.0)+math.pow(b*sn,2.0))
        h = p/cs - N
    llh = {'lon':clambda, 'lat':theta, 'height': h}
    return llh
Related Question