GPS to Cartesian

I'm interested in turning GPS position into Cartesian, I want to then transform it so the starting point is the origin.

Anyone have experience with this?

• I've been using the geodesic and distance functions in geopy in python geopy.readthedocs.io/en/stable/#module-geopy.distance, where two calls to geodesic() get x and y offsets in meters from starting lat long location (and decomposing x y meters into two distance() calls would go in reverse to get lat long)- if this seems promising I can expand to full answer, but I'm not super confident how good this method is even over small distances. Some additional code here: github.com/geopy/geopy/discussions/506 Dec 3, 2022 at 1:49
• Let us know if you want a more complete answer, but GPS is usually given in latitude and longitude. It can be given in MGRS (military way of doing the same thing). After accounting for where you are on earth, both can be used for short distances w/o too much regard to error as X & Y offsets. But if you travel far, you need to consider using the Great Circle Distance Formula. Dec 4, 2022 at 17:16
• I would like a more complete answer, GPS should also give altitude so you get a 3d coordinate? Dec 4, 2022 at 20:23

It sounds like maybe you're looking for UTM projections. From the Wiki article:

The Universal Transverse Mercator (UTM) is a map projection system for assigning coordinates to locations on the surface of the Earth.

You're not clear with your purpose here so it's not really possible to give you more context related to your use case. Generally speaking, though, you convert latitude and longitude to a UTM easting (x) and northing (y) value. Your altitude can stay as elevation, and you've got x/y/z values.

This is a projection, so like with any project there is distortion. Again quoting the Wiki entry:

By using narrow zones of 6° of longitude (up to 668 km) in width, and reducing the scale factor along the central meridian to 0.9996 (a reduction of 1:2500), the amount of distortion is held below 1 part in 1,000 inside each zone. Distortion of scale increases to 1.0010 at the zone boundaries along the equator.

You can get into trouble when you cross UTM zones, but then it's your choice on how you want to handle that, if it's swapping zones or extending the origin zone's projection to the outlier point(s).

This is a geopy https://geopy.readthedocs.io/en/stable/#module-geopy.distance method for going back and forth from a coordinate and a reference coordinate to relative meters in north and east:

#!/usr/bin/env python
# convert latitude-longitude-altitude to meters relative to reference point

import geopy
# geopy.distance doesn't work unless this is also done
from geopy import distance

# x is north and y is east
def coords_to_relative_meters(lat0, lng0, lat1, lng1, ellipsoid="WGS-84", north_then_east=False):
xa = distance.geodesic((lat0, lng0), (lat1, lng0), ellipsoid=ellipsoid).m
ya = distance.geodesic((lat1, lng0), (lat1, lng1), ellipsoid=ellipsoid).m
print(f"north-then-east relative x north {xa}m, y east {ya}m")

yb = distance.geodesic((lat0, lng0), (lat0, lng1), ellipsoid=ellipsoid).m
xb = distance.geodesic((lat0, lng1), (lat1, lng1), ellipsoid=ellipsoid).m
print(f"east-then-north relative x north {xb}m, y east {yb}m")

print(f"north-then-east vs. east-then-north error {(xb - xa):0.6f}m {(yb - ya):0.6f}m")
if north_then_east:
return xa, ya
else:
return xb, yb

def relative_meters_to_coords(lat0, lng0, rel_x, rel_y, ellipsoid="WGS-84", north_then_east=False):
# 0 is north (rel_x), 90 is east (rel_y)
ll_na = distance.geodesic(meters=rel_x, ellipsoid=ellipsoid).destination((lat0, lng0), bearing=0.0)
ll_ea = distance.geodesic(meters=rel_y, ellipsoid=ellipsoid).destination(ll_na, bearing=90.0)
print(f"north-then-east relative x north {ll_ea.latitude}° {ll_ea.longitude}°")

ll_eb = distance.geodesic(meters=rel_y, ellipsoid=ellipsoid).destination((lat0, lng0), bearing=90.0)
ll_nb = distance.geodesic(meters=rel_x, ellipsoid=ellipsoid).destination(ll_eb, bearing=0.0)
print(f"east-then-north relative x north {ll_nb.latitude}° {ll_ea.longitude}°")

text = "north-then-east vs. east-then-north error"
text += f" lat {(ll_ea.latitude - ll_nb.latitude):0.6f}°"
text += f", long {(ll_ea.longitude - ll_nb.longitude):0.6f}°"
print(text)

if north_then_east:
rv = ll_ea
else:
rv = ll_nb

return rv.latitude, rv.longitude

if __name__ == "__main__":
print(geopy.__version__)

lat0 = 46.8454439
lng0 = -121.7649698
print(f"reference {lat0}° {lng0}°")

lat1 = lat0 + 0.01
lng1 = lng0 + 0.01
print(f"target    {lat1}° {lng1}°")

# if this is set false then round trip error much lower because of the way the conversion
# back to lat to long works
north_then_east = True
print(f"relative coordinates {lat1 - lat0}° {lng1 - lng0}°")
for ellipsoid in distance.ELLIPSOIDS.keys():
print(f"\n{ellipsoid}")
rel_x, rel_y = coords_to_relative_meters(lat0, lng0, lat1, lng1,
ellipsoid=ellipsoid, north_then_east=north_then_east)
lat1b, lng1b = relative_meters_to_coords(lat0, lng0, rel_x, rel_y,
ellipsoid=ellipsoid, north_then_east=north_then_east)
print(f"round trip conversion error lat {(lat1b - lat1)}°, lng {(lng1b - lng1)}°")


and the output for WGS-84:

reference 46.8454439° -121.7649698°
target    46.8554439° -121.7549698°
relative coordinates 0.00999999999999801° 0.010000000000005116°

WGS-84
north-then-east relative x north 1111.6792380437344m, y east 762.6088552752285m
east-then-north relative x north 1111.6792380437344m, y east 762.7504100234593m
north-then-east vs. east-then-north error 0.000000m 0.141555m

north-then-east relative x north 46.85544346321024° -121.7549698000608°
east-then-north relative x north 46.85544346336275° -121.7549698000608°
north-then-east vs. east-then-north error lat -0.000000°, long 0.000002°
round trip conversion error lat -4.367897545876076e-07°, lng -6.080824732634937e-11°


The bigger the change in latitude the more error the relative y/east meters has- 0.15 meters when for example about a kilometer is traversed moving north from the reference point to the target point using the above output. As long as everything has good lat/long coordinates to start and they share the reference point and use consistent coordinate conversion functions it may not matter much (it's a little useless if a gnss sensor sdk has the option to output relative meters but doesn't expose a call to convert any coordinate using the same function, or document what they did to do the conversion, you wouldn't want to commingle two different sensors with two conversion methods).

If you had really good relative distance measurements and wanted to convert in the other direction, get latitude longitude given distance traveled from a reference point it would be better to break the path up into segments than take the endpoints and ignore the path between them if the amount of error is a concern.