Calculate distance between two lat/lon/alt points in Python

I'm trying to find out how to calculate the distance between two points with latitude, longitude, and elevation data.

I found that `geographiclib.geodesic.Geodesic.WGS84.Inverse()` outputs a `s12` parameter that has the distance data according to the geographiclib documentation. I can also see a Wikipedia entry on altitude correction.

Is there any pre-built function that can take all three dimentions and output the correct distance or would the solution be to run a correction on the 2d method output?

I'm very new to this, so I'm also trying to learn and am open to using a different library

2 Answers

As outlined in Calculating distance between two points using latitude longitude and altitude (elevation), Geopy: Calculating Distance, [GPS] Use pyproj to calculate distance, azimuth, and elevation from GPS latitude and longitude or Haversine formula that includes an altitude parameter?, you need to calculate first the 2D great-circle distance or the geodesic distance between the 2 points (distance between flat coordinates) and combine it with the difference in altitudes inside the Euclidean Formula for the 3D distance calculation. The solution assumes that the altitude is in meters and thus the great_circle's or geodesic's distance needs to be in meters

You can use GeoPy, pyproj or geographiclib for example.

With GeoPy

``````import numpy as np
# points with latitude,longitude units in degrees ,altitude unit in meters
pt1 =[35.3524,135.0302, 100]
pt2 =[35.3532,135.0305,500]
# 2D geodesic distance in meters
from geopy import distance
distance_2d= distance.distance(pt1[:2], pt2[:2]).m
print(distance_2d)
92.85194331754518
# 3D euclidean distance
distance_3d = np.sqrt(distance_2d**2 + (pt2[2] - pt1[2])**2)
print(distance_3d)
410.6354628838632
``````

With pyproj

``````from pyproj import Geod
g = Geod(ellps='WGS84')
# 2D distance in meters with longitude, latitude of the points
azimuth1, azimuth2, distance_2d = g.inv(pt1[1], pt1[0], pt2[1], pt2[0])
print(distance_2d)
92.85194331754519
# 3D euclidean distance
distance_3d = np.hypot(distance,pt2[2]-pt1[2])
print(distance_3d)
410.63546288386
``````

With geographiclib

``````from geographiclib.geodesic import Geodesic
geod = Geodesic.WGS84
g = geod.Inverse(pt1[0], pt1[1],pt2[0],pt2[1])
#2D geodesic distance in meters
print(g['s12'])
92.85194331754519
# 3D euclidean distance
distance_3d = np.hypot(g['s12'],pt2[2]-pt1[2])
print(distance_3d)
410.63546288386
``````
• This is probably good enough for most situations, but technically it's not quite perfect because the 2D distance forms an arc, not a line segment, so using the Pythagorean Theorem to find the 3D distance will introduce a small amount of error. Over long distances this might matter in some use cases. Commented Mar 24, 2022 at 2:27
• The pyproj code seems to have a typo, saying `py2` instead of `pt2`. Commented Sep 1, 2023 at 21:35

This discards any topographical differences between the two points even for the 3D. Even in short distances (few hundred of meters and more) this induces significant errors in distance, if you're looking for the measured distance (on land) and not the horizontal one. A fix for this could be to use a GIS map provider (such as SRTM, ESRI, BING...)