Converting from longitude\latitude to Cartesian coordinates
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.
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
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,