Converting from longitude\latitude to Cartesian coordinates

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,

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 

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) 

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)) 

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]) 

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).

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.

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; } 

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);