3

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

4

2 Answers 2

7

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
2
  • 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.
    – Ethan T
    Mar 24, 2022 at 2:27
  • The pyproj code seems to have a typo, saying py2 instead of pt2.
    – Matt F.
    Sep 1 at 21:35
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.