[GIS] Converting ECEF coordinates into geodetic LLA

coordinate system

I am trying to convert ECEF coordinates into geodetic LLA.

This is the coordinates:

X = 1199816.63380
Y = -4812201.6785
Z = 4016140.62425

I checked the coordinates using this and it returns correct values:

Lat: 39.111677542
Long: -76
Alt: 12068.6198

But using my code, I get these values:

Lat: 179.79495170382071
Long: -75.99999999999999
Alt: NaN

And this is my code:

const RADIUS_EQUATOR = 6378137.0;
const RADIUS_POLES = 6356752.3;
const e = Math.sqrt((RADIUS_EQUATOR*RADIUS_EQUATOR / RADIUS_POLES*RADIUS_POLES) / RADIUS_EQUATOR*RADIUS_EQUATOR);
const eprime = Math.sqrt((RADIUS_EQUATOR*RADIUS_EQUATOR - RADIUS_POLES*RADIUS_POLES) / RADIUS_POLES*RADIUS_POLES);

function radianToDegree(r) {
  return r * (180 / Math.PI);
}

function ECRtoLLA(ECR) {
  let X = ECR.x;
  let Y = ECR.y;
  let Z = ECR.z;

  let a = RADIUS_EQUATOR;
  let b = RADIUS_POLES;

  //Auxiliary values first
  let p = Math.sqrt(X*X + Y*Y);
  let theta = Math.atan((Z*a) / (p*b));

  let sintheta = Math.sin(theta);
  let costheta = Math.cos(theta);

  let num = Z + eprime * eprime * b * sintheta * sintheta * sintheta;
  let denom = p - e * e * a * costheta * costheta * costheta;

  //Now calculate LLA
  let latitude  = Math.atan2(num, denom);
  let longitude = Math.atan2(Y, X);
  let N = getN(latitude, a);
  let altitude  = (p / Math.cos(latitude)) - N;

  if (X < 0 && Y < 0) {
    longitude = longitude - Math.PI;
  }

  if (X < 0 && Y > 0) {
    longitude = longitude + Math.PI;
  }

  latitude = radianToDegree(latitude);
  longitude = radianToDegree(longitude);

  return {
    latitude: latitude,
    longitude: longitude,
    altitude: altitude
  };
}

console.info(ECRtoLLA({
  x: 1199816.63380,
  y: -4812201.6785,
  z: 4016140.62425
}));

Any pointer or help?

Best Answer

Seems like I calculate wrong e and eprime.

After fixing them, I got the correct results. Sorry this is the correction

const RADIUS_EQUATOR = 6378137.0;
const a = RADIUS_EQUATOR;
const asqr = a*a;

const FLATTENING_DENOM = 298.257223563
const FLATTENING = 1 / FLATTENING_DENOM;
const RADIUS_POLES = RADIUS_EQUATOR*(1-FLATTENING);
const b = RADIUS_POLES;
const bsqr = b*b;

const e = Math.sqrt((asqr - bsqr) / asqr);
const eprime = Math.sqrt((asqr - bsqr) / bsqr);
Related Question