// Converts cartesian coordinates (x,y,z) to geodetic coordinates (lat,lon,h)
// with standard deviations (slat,slon,sh), covariance (covl) and the WGS84 ellipsoid
// Assumes x-axis points towards 0° latitude and 0° longitude, and y-axis points towards 0° latitude and 90°E longitude
// Input units are meters, output units are radians (lat, lon) and meters (h, slat, slon, sh)
// If fast is true, a faster approximation is used (WGS84 ellipsoid only)
function xyz2latlon(x,y,z,sx,sy,sz,sxy,syz,sxz,fast) {
const a = 6378137; // WGS84 semimajor axis (m)
const finv = 298.257223563; // WGS84 inverse flattening
const b = a * (finv - 1) / finv; // WGS84 semiminor axis (m)
const e2 = (a**2 - b**2) / a**2; // First eccentricity squared
const ep2 = (a**2 - b**2) / b**2; // Second eccentricity squared
// Quick approximation (faster but less accurate)
if (fast) {
const lon = Math.atan2(y,x); // Longitude (rad)
const p = Math.sqrt(x**2 + y**2); // Distance from z-axis (m)
const lat = Math.atan2(z,p*(1-e2)); // Latitude (rad)
const sinlat = Math.sin(lat);
const coslat = Math.cos(lat);
const N = a / Math.sqrt(1 - e2*sinlat**2); // Prime vertical radius of curvature
const M = a*(1-e2) / (1 - e2*sinlat**2)**(3/2); // Meridian radius of curvature
const slon = Math.sqrt(sx**2 + sy**2) / (p*coslat); // Standart deviation of longitude (rad)
const slat = Math.sqrt(sz**2*(p**2*coslat**2 + z**2*sinlat**2) / (p**2*(1 - e2)*coslat**2 + z**2*(1 - e2*sinlat**2))) / (N*(1 - e2*sinlat**2)); // Standart deviation of latitude (rad)
const sh = Math.sqrt(sz**2) / (N*(1 - e2*sinlat**2)); // Standart deviation of height (m)
covl = sxy / (p*coslat*slon*slat); // Covariance between latitude and longitude (rad^2)
return [lat, lon, z, slat, slon, sz, covl];
}
// Full algorithm (slower but more accurate)
const p = Math.sqrt(x**2 + y**2); // Distance from z-axis (m)
const theta = Math.atan2(z*a, p*b); // Geocentric latitude (rad)
const sintheta = Math.sin(theta);
const costheta = Math.cos(theta);
const lon = Math.atan2(y,x); // Longitude (rad)
const N = a / Math.sqrt(1 - e2*sintheta**2); // Prime vertical radius of curvature
const h = p/costheta - N; // Height above ellipsoid (m)
const tantheta = (z + ep2*N*sintheta) / p; // Tangent of geocentric latitude
const lat = Math.atan2(z + ep2*N*sintheta**3, p - e2*a*N*costheta**3); // Geodetic latitude (rad)
const sinlat = Math.sin(lat);
const coslat = Math.cos(lat);
const slon = Math.sqrt(sx**2 + sy**2) / (p*coslat); // Standart deviation of longitude (rad)
const slat = Math.sqrt((sz**2*(p**2 + (z+ep2*N*sintheta)**2*ep2*coslat**2*sintheta**2)) / (p**2*(1 - e2*sinlat**2)**2*coslat**2)) / N; // Standart deviation of latitude (rad)
const sh = Math.sqrt(sz**2) / (sintheta*N); // Standart deviation of height (m)
const D = -sintheta*coslat*(sx**2*sinlon**2 - 2*sxy*sinlon*coslon + sy**2*coslon**2) + p*coslat*(sintheta*coslat*sxz - costheta*coslat*syz - sintheta*sinlat*sxy);
const covll = D / ((p*N*coslat)**2 * slon**2 * slat**2); // Covariance between latitude and longitude (rad^2)
const covhh = sz**2 / sh**4; // Covariance between height and height (m^2)
const covlh = -sintheta*(syz*costheta*slat*coslat*sinlon**2 + sxz*costheta*slat*coslat*sinlon*coslon - sxy*sintheta*slat*sinlat*coslon - sxy*sintheta*slat*sinlon*coslat + syz*costheta*salt*coslat*coslon - sxz*costheta*clat*sinlat*sinlon + sxy*sintheta*clat*coslat*coslon - syz*costheta*clat*coslat*sinlon) / (N**2 * slat**2 * sh**2); // Covariance between latitude and height (m*rad)
covl = covll + covlh; // Total covariance between latitude and longitude (rad^2)
return [lat, lon, h, slat, slon, sh, covl];
}
Note that this function assumes the WGS84 ellipsoid and the axes orientation described in the code comments. If you need to use a different ellipsoid or axis orientation, you will need to modify the constants and formulas accordingly.