How to calculate 3D distance (including altitude) between two points in GeoDjango

Question

Prologue:

This is a question arising often in SO:

• 3d distance calculations with GeoDjango
• Calculating distance between two points using latitude longitude and altitude (elevation)
• Distance between two 3D point in geodjango (postgis)

I wanted to compose an example on SO Documentation but the geodjango chapter never took off and since the Documentation got shut down on August 8, 2017, I will follow the suggestion of this widely upvoted and discussed meta answer and write my example as a self-answered post.

Of course, I would be more than happy to see any different approach as well!!

---

Question:

Assume the model:

class MyModel(models.Model):
    name = models.CharField()
    coordinates = models.PointField()

Where I store the point in the coordinate variable as a lan, lng, alt point:

MyModel.objects.create(
    name='point_name', 
    coordinates='SRID=3857;POINT Z (100.00 10.00 150)')

I am trying to calculate the 3D distance between two such points:

p1 = MyModel.objects.get(name='point_1').coordinates
p2 = MyModel.objects.get(name='point_2').coordinates

d = Distance(m=p1.distance(p2))

Now d=X in meters.

If I change only the altitude of one of the points in question:

For example:

p1.coordinates = 'SRID=3857;POINT Z (100.00 10.00 200)'

from 150 previously, the calculation:

d = Distance(m=p1.distance(p2))

returns d=X again, like the elevation is ignored.
How can I calculate the 3D distance between my points?

Answer 1

Reading from the documentation on the GEOSGeometry.distance method:

Returns the distance between the closest points on this geometry and the given geom (another GEOSGeometry object).

Note

GEOS distance calculations are linear – in other words, GEOS does not perform a spherical calculation even if the SRID specifies a geographic coordinate system.

Therefore we need to implement a method to calculate a more accurate 2D distance between 2 points and then we can try to apply the altitude (Z) difference between those points.

1. Great-Circle 2D distance calculation (Take a look at the 2022 UPDATE below the explanation for a better approach using geopy):

The most common way to calculate the distance between 2 points on the surface of a sphere (as the Earth is simplistically but usually modeled) is the Haversine formula:

The haversine formula determines the great-circle distance between two points on a sphere given their longitudes and latitudes.

Although from the great-circle distance wiki page we read:

Although this formula is accurate for most distances on a sphere, it too suffers from rounding errors for the special (and somewhat unusual) case of antipodal points (on opposite ends of the sphere). A formula that is accurate for all distances is the following special case of the Vincenty formula for an ellipsoid with equal major and minor axes.

We can create our own implementation of the Haversine or the Vincenty formula (as shown here for Haversine: Haversine Formula in Python (Bearing and Distance between two GPS points)) or we can use one of the already implemented methods contained in geopy:

• geopy.distance.great_circle (Haversine):

      from geopy.distance import great_circle
      newport_ri = (41.49008, -71.312796)
      cleveland_oh = (41.499498, -81.695391)
  
      # This call will result in 536.997990696 miles
      great_circle(newport_ri, cleveland_oh).miles) 
  

• geopy.distance.vincenty (Vincenty):

      from geopy.distance import vincenty
      newport_ri = (41.49008, -71.312796)
      cleveland_oh = (41.499498, -81.695391)
  
      # This call will result in 536.997990696 miles
      vincenty(newport_ri, cleveland_oh).miles
  

!!!2022 UPDATE: On 2D distance calculation using geopy:

GeoPy discourages the use of Vincenty as of version 1.14.0. Changelog states:

CHANGED: Vincenty usage now issues a warning. Geodesic should be used instead. Vincenty is planned to be removed in geopy 2.0. (#293)

So (especially if we are going to apply the calculation on a WGS84 ellipsoid) we should use geodesic distance instead:

from geopy.distance import geodesic
newport_ri = (41.49008, -71.312796)
cleveland_oh = (41.499498, -81.695391)

# This call will result in 538.390445368 miles
geodesic(newport_ri, cleveland_oh).miles

2. Adding altitude to the mix:

As mentioned, each of the above calculations yields a great circle distance between 2 points. That distance is also called "as the crow flies", assuming that the "crow" flies without changing altitude and as straight as possible from point A to point B.

We can have a better estimation of the "walking/driving" ("as the crow walks"??) distance by combining the result of one of the previous methods with the difference (delta) in altitude between point A and point B, inside the Euclidean Formula for distance calculation:

acw_dist = sqrt(great_circle(p1, p2).m**2 + (p1.z - p2.z)**2)

---

The previous solution is prone to errors especially the longer the real distance between the points is.
I leave it here for comment continuation reasons.

GeoDjango Distance calculates the 2D distance between two points and doesn't take into consideration the altitude differences.
In order to get the 3D calculation, we need to create a distance function that will consider altitude differences in the calculation:

Theory:

The latitude, longitude and altitude are Polar coordinates and we need to translate them to Cartesian coordinates (x, y, z) in order to apply the Euclidean Formula on them and calculate their 3D distance.

• Assume:
  polar_point_1 = (long_1, lat_1, alt_1)
  and
  polar_point_2 = (long_2, lat_2, alt_2)

• Translate each point to it's Cartesian equivalent by utilizing this formula:

   x = alt * cos(lat) * sin(long)
   y = alt * sin(lat)
   z = alt * cos(lat) * cos(long)
  

and you will have p_1 = (x_1, y_1, z_1) and p_2 = (x_2, y_2, z_2) points respectively.

• Finally use the Euclidean formula:

   dist = sqrt((x_2-x_1)**2 + (y_2-y_1)**2 + (z_2-z_1)**2)
  

Answer 2

Using geopy, this is the easiest and perfect solution.

https://geopy.readthedocs.io/en/stable/#geopy.distance.lonlat

>>> from geopy.distance import distance
>>> from geopy.point import Point
>>> a = Point(-71.312796, 41.49008, 0)
>>> b = Point(-81.695391, 41.499498, 0)
>>> print(distance(a, b).miles)
538.3904453677203

Answer 3

Once converted into Cartesian coordinates, you can compute the norm with numpy:

np.linalg.norm(point_1 - point_2)