Converting from longitude\latitude to Cartesian coordinates

Question

I have some earth-centered coordinate points given as latitude and longitude (WGS-84).

How can i convert them to Cartesian coordinates (x,y,z) with the origin at the center of the earth?

Bromoform

Have you managed to convert WGS-84 longitude and latitude into Cartesian coordinates ?. I also have elevation. I have tried the accepted answer here, but it doesn't give me the right answer. I compared my results with this website: apsalin.com/convert-geodetic-to-cartesian.aspx. – Antecedence

Answer 1

I have recently done something similar to this using the "Haversine Formula" on WGS-84 data, which is a derivative of the "Law of Haversines" with very satisfying results.

Yes, WGS-84 assumes the Earth is an ellipsoid, but I believe you only get about a 0.5% average error using an approach like the "Haversine Formula", which may be an acceptable amount of error in your case. You will always have some amount of error unless you're talking about a distance of a few feet and even then there is theoretically curvature of the Earth... If you require a more rigidly WGS-84 compatible approach checkout the "Vincenty Formula."

I understand where starblue is coming from, but good software engineering is often about trade-offs, so it all depends on the accuracy you require for what you are doing. For example, the result calculated from "Manhattan Distance Formula" versus the result from the "Distance Formula" can be better for certain situations as it is computationally less expensive. Think "which point is closest?" scenarios where you don't need a precise distance measurement.

Regarding, the "Haversine Formula" it is easy to implement and is nice because it is using "Spherical Trigonometry" instead of a "Law of Cosines" based approach which is based on two-dimensional trigonometry, therefore you get a nice balance of accuracy over complexity.

A gentleman by the name of Chris Veness has a great website that explains some of the concepts you are interested in and demonstrates various programmatic implementations; this should answer your x/y conversion question as well.

Pyrometallurgy

0.5% error - 0.5% of what? In the context of this question it could be the radius of the earth, so 0.5% could be 30 km :) – Condense

Checked out your link. The 0.5% quote is for error in the great-circle distance between two points so not strictly relevant to this question. I would think when converting lat-long to Cartesian coordinates with the origin at the centre of the earth, the errors from assuming a spherical earth could be significant. It's not clear what the questionner wants to do with the Cartesian coordinates. Either it's just more convenient to work in them for some bizarre reason, or perhaps it's some requirement for data export? If the latter, accuracy would be important. – Condense

Answer 2

Here's the answer I found:

Just to make the definition complete, in the Cartesian coordinate system:

• the x-axis goes through long,lat (0,0), so longitude 0 meets the equator;
• the y-axis goes through (0,90);
• and the z-axis goes through the poles.

The conversion is:

x = R * cos(lat) * cos(lon)

y = R * cos(lat) * sin(lon)

z = R *sin(lat)

Where R is the approximate radius of earth (e.g. 6371 km).

If your trigonometric functions expect radians (which they probably do), you will need to convert your longitude and latitude to radians first. You obviously need a decimal representation, not degrees\minutes\seconds (see e.g. here about conversion).

The formula for back conversion:

lat = asin(z / R)
lon = atan2(y, x)

asin is of course arc sine. read about atan2 in wikipedia. Don’t forget to convert back from radians to degrees.

This page gives c# code for this (note that it is very different from the formulas), and also some explanation and nice diagram of why this is correct,

Bromoform

-1 This is wrong. You are assuming the earth is a sphere, while WGS-84 assumes an ellipsoid. – Sefton

@starblue: I'm not sure you are in a position to label the given answer "right" or "wrong". The spherical approximation (to get ECEF-style x,y,z coords) using available lat/lngs (that are referenced to WGS-84) is either "adequate" for the original poster's needs, or "not adequate". For distance and bearing estimates, I'd bet this simple conversion is fine. If he is launching satellites, maybe not. After all, WGS-84 itself is "wrong"... in that it is not a perfect model of the earth's surface; all ellipsoidal models are approximations. Too bad the OP didn't tell us what he was trying to do. – Theatricalize

@Dan H The question asks for WGS-84, and if you answer something else you should at least discuss the differences/error, which this answer doesn't. – Sefton

@daphna-shezaf can't make a back conversion... I have done also the back from radians to degrees, but result is not the same... – Efferent

thank you, spend hour figuring out why it doesn't work, turns out i swapped some cos(lat) and sin(lat) – Alienable

Is it possible to do the conversion from cartesian to regular without using z? if you only have x and y from the start? – Dismissal

What about the altitude? R = R + h? – Blanche

Answer 3

In python3.x it can be done using :

# Converting lat/long to cartesian
import numpy as np

def get_cartesian(lat=None,lon=None):
    lat, lon = np.deg2rad(lat), np.deg2rad(lon)
    R = 6371 # radius of the earth
    x = R * np.cos(lat) * np.cos(lon)
    y = R * np.cos(lat) * np.sin(lon)
    z = R *np.sin(lat)
    return x,y,z

Degeneration

If the value of lat, lon are in decimal representation (i.e. meters unit), do we need to convert R value from Km to m to calculate the x and y? – Carbonaceous

I used this approach to plot points onto the globe using CesiumJS, which uses cartesian as opposed to lat/long/altitude. I also needed to convert km to meters...simple! Just add * 1000 at the end of each line where the x, y, and z variables are used – Hom

Answer 4

Theory for convert GPS（WGS84) to Cartesian coordinates https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates

The following is what I am using:

• Longitude in GPS（WGS84） and Cartesian coordinates are the same.
• Latitude need be converted by WGS 84 ellipsoid parameters semi-major axis is 6378137 m, and
• Reciprocal of flattening is 298.257223563.

I attached a VB code I wrote:

Imports System.Math

'Input GPSLatitude is WGS84 Latitude,h is altitude above the WGS 84 ellipsoid

Public Function GetSphericalLatitude(ByVal GPSLatitude As Double, ByVal h As Double) As Double

        Dim A As Double = 6378137 'semi-major axis 
        Dim f As Double = 1 / 298.257223563  '1/f Reciprocal of flattening
        Dim e2 As Double = f * (2 - f)
        Dim Rc As Double = A / (Sqrt(1 - e2 * (Sin(GPSLatitude * PI / 180) ^ 2)))
        Dim p As Double = (Rc + h) * Cos(GPSLatitude * PI / 180)
        Dim z As Double = (Rc * (1 - e2) + h) * Sin(GPSLatitude * PI / 180)
        Dim r As Double = Sqrt(p ^ 2 + z ^ 2)
        Dim SphericalLatitude As Double =  Asin(z / r) * 180 / PI
        Return SphericalLatitude
End Function

Please notice that the h is altitude above the WGS 84 ellipsoid.

Usually GPS will give us H of above MSL height. The MSL height has to be converted to height h above the WGS 84 ellipsoid by using the geopotential model EGM96 (Lemoine et al, 1998).
This is done by interpolating a grid of the geoid height file with a spatial resolution of 15 arc-minutes.

Or if you have some level professional GPS has Altitude H (msl,heigh above mean sea level) and UNDULATION,the relationship between the geoid and the ellipsoid (m) of the chosen datum output from internal table. you can get h = H(msl) + undulation

To XYZ by Cartesian coordinates:

x = R * cos(lat) * cos(lon)

y = R * cos(lat) * sin(lon)

z = R *sin(lat)

Inroad

What is the value of R? – Raconteur

I guess it's the radius of the sphere, which is 6371km for the earth. – Multistage

Answer 5

The proj.4 software provides a command line program that can do the conversion, e.g.

LAT=40
LON=-110
echo $LON $LAT | cs2cs +proj=latlong +datum=WGS84 +to +proj=geocent +datum=WGS84

It also provides a C API. The C equivalent of the above would be something like

#include <math.h>
#include <proj.h>
#include <stdio.h>

int main(int argc, char* argv[])
{

    PJ* P = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
        "+proj=latlong +datum=WGS84",
        "+proj=geocent +datum=WGS84", NULL);

    if (P == NULL) {
        return 1;
    }

    PJ_COORD c;
    c.lpzt.lam = -110.0;  /* Longitude in degrees */
    c.lpzt.phi = 40.0;    /* Latitude in degrees */
    c.lpzt.z = 0.0;       /* Height in meters */
    c.lpzt.t = HUGE_VAL;  /* Time (not used here) */

    PJ_COORD c_out = proj_trans(P, PJ_FWD, c);

    printf("X = %.17g m\n", c_out.xyzt.x);
    printf("Y = %.17g m\n", c_out.xyzt.y);
    printf("Z = %.17g m\n", c_out.xyzt.z);

    return 0;
}

Resulting in the output

X = -1673404.5546274509 m
Y = -4597641.2274514409 m
Z = 4077985.5722003761 m

There's also a Python package available:

from pyproj import Transformer
lon, lat, hgt = -110., 40., 0.
llh_to_xyz = Transformer.from_crs("+proj=latlong +datum=WGS84", "+proj=geocent +datum=WGS84")
print(llh_to_xyz.transform(lon, lat, hgt))

Dittany

It looks like pj_geodetic_to_geocentric() is no longer in the API. – Overblown

It seems that pj_geodetic_to_geocentric was deprecated in 2019 with the release of PROJ version 6. I've updated the answer with a C example that works with PROJ version 9.2.0. – Dittany

Answer 6

If you care about getting coordinates based on an ellipsoid rather than a sphere, take a look at Geographic_coordinate_conversion - it gives the formulae. GEodetic Datum has the WGS84 constants you need for the conversion.

The formulae there also take into account the altitude relative to the reference ellipsoid surface (useful if you are getting altitude data from a GPS device).

Smilacaceous

Answer 7

I made a function in Python that takes into account the fact that the earth is not a perfect sphere. References are in the comments:

    # this function converts latitude,longitude and height above sea level 
    # to earthcentered xyx coordinates in wgs84, lat and lon in decimal degrees 
    # e.g. 52.724156(West and South are negative), heigth in meters
    # for algoritm see https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates
    # for values of a and b see https://en.wikipedia.org/wiki/Earth_radius#Radius_of_curvature

from math import *

def latlonhtoxyzwgs84(lat,lon,h):


    a=6378137.0             #radius a of earth in meters cfr WGS84
    b=6356752.3             #radius b of earth in meters cfr WGS84
    e2=1-(b**2/a**2)
    latr=lat/90*0.5*pi      #latitude in radians
    lonr=lon/180*pi         #longituede in radians
    Nphi=a/sqrt(1-e2*sin(latr)**2)
    x=(Nphi+h)*cos(latr)*cos(lonr)
    y=(Nphi+h)*cos(latr)*sin(lonr)
    z=(b**2/a**2*Nphi+h)*sin(latr)
    return([x,y,z])

Hiller

Could you specify what is the h parameter. Is it the altitude of the point above sea level or the distance to the center of Earth ? – Perfect

It is the height above sea level – Hiller

Answer 8

Why implement something which has already been implemented and test-proven?

C#, for one, has the NetTopologySuite which is the .NET port of the JTS Topology Suite.

Specifically, you have a severe flaw in your calculation. The earth is not a perfect sphere, and the approximation of the earth's radius might not cut it for precise measurements.

If in some cases it's acceptable to use homebrew functions, GIS is a good example of a field in which it is much preferred to use a reliable, test-proven library.

Ingratitude

+1. Using a reliable library is more accurate than a homebrew function and also easier. – Condense

How does NetTopologySuite convert from long/late to cartesion? – Ferrous

NTS doesn't include coordinate conversion capabilities, maybe you need Proj.NET projnet.codeplex.com – Contraception

Ridiculous, answer doesn't even provide conversion capability. – Devi

Answer 9

Coordinate[] coordinates = new Coordinate[3];
coordinates[0] = new Coordinate(102, 26);
coordinates[1] = new Coordinate(103, 25.12);
coordinates[2] = new Coordinate(104, 16.11);
CoordinateSequence coordinateSequence = new CoordinateArraySequence(coordinates);

Geometry geo = new LineString(coordinateSequence, geometryFactory);

CoordinateReferenceSystem wgs84 = DefaultGeographicCRS.WGS84;
CoordinateReferenceSystem cartesinaCrs = DefaultGeocentricCRS.CARTESIAN;

MathTransform mathTransform = CRS.findMathTransform(wgs84, cartesinaCrs, true);

Geometry geo1 = JTS.transform(geo, mathTransform);

Dalia

Would you be able to elaborate? I created a simple application that tiers to transform a single coordinate using your approach. It always fails though as the dimensions of the source (2) and the dimensions of the target (3) differ, resulting in an exception java.lang.IllegalArgumentException: dimension must be <= 3 – Glamorous

Hmmm... I've looked at JTS for a bit. The lines upto and including the new LineString() look like JTS. But I don't see the CRS and Transform stuff in JTS. So: are they there and I'm missing them? Were there, and removed in 1.12? Or: is that a different library? – Theatricalize

Answer 10

You can do it this way on Java.

public List<Double> convertGpsToECEF(double lat, double longi, float alt) {
    
    double a=6378.1;
    double b=6356.8;
    double N;
    double e= 1-(Math.pow(b, 2)/Math.pow(a, 2));
    N= a/(Math.sqrt(1.0-(e*Math.pow(Math.sin(Math.toRadians(lat)), 2))));
    double cosLatRad=Math.cos(Math.toRadians(lat));
    double cosLongiRad=Math.cos(Math.toRadians(longi));
    double sinLatRad=Math.sin(Math.toRadians(lat));
    double sinLongiRad=Math.sin(Math.toRadians(longi));
    double x =(N+0.001*alt)*cosLatRad*cosLongiRad;
    double y =(N+0.001*alt)*cosLatRad*sinLongiRad;
    double z =((Math.pow(b, 2)/Math.pow(a, 2))*N+0.001*alt)*sinLatRad;
    
    List<Double> ecef= new ArrayList<>();
    ecef.add(x);
    ecef.add(y);
    ecef.add(z);
    
    return ecef;
    
    
}

Abby

what is parameter alt ? – Mcferren

altitude, what are you even doing here if you don't know how GPS works ;) – Wolfram

The question does not mention GPS or altitude at all. – Wenger