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

Question

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

Answer 1

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

Answer 2

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