如何计算 GeoDjango 中两点之间的 3D 距离（包括高度）

Question

序幕：

这是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)

我想编写一个关于 SO 文档的示例，但 geodjango 章节从未起飞，并且由于文档于 2017 年 8 月 8 日关闭，我将遵循 this widely upvoted and discussed meta answer 的建议并将我的示例编写为自我- 已回答的帖子。

当然，我也很乐意看到任何不同的方法！

---

问题：

假设模型：

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

我将coordinate 变量中的点存储为lan, lng, alt 点：

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

我正在尝试计算两个这样的点之间的 3D 距离：

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

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

现在d=X 以米为单位。

如果我只更改相关点之一的高度：

例如：

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

从之前的150，计算：

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

再次返回d=X，就像忽略高度一样。
如何计算我的点之间的 3D 距离？

Answer 1

阅读GEOSGeometry.distance 方法的文档：

返回此几何图形上最近点与给定几何图形（另一个 GEOSGeometry 对象）之间的距离。

注意

GEOS 距离计算是线性的——换句话说，即使 SRID 指定了地理坐标系，GEOS 也不执行球面计算。

因此，我们需要实现一种方法来计算 2 点之间更准确的 2D 距离，然后我们可以尝试应用这些点之间的高度 (Z) 差。

1.大圆二维距离计算：

计算球体表面上两点之间距离的最常用方法（因为地球很简单，但通常是建模的）是Haversine formula：

harsine 公式根据它们的经度和纬度确定球体上两点之间的great-circle distance。

虽然我们从great-circle distance wiki page 读到：

虽然这个公式对于球体上的大多数距离都是准确的，但对于对映点的特殊（并且有些不寻常）情况（在球体的相对端），它也存在舍入误差。一个对所有距离都准确的公式是 Vincenty formula 的以下特例，适用于长轴和短轴相等的椭圆体。

我们可以创建我们自己的Haversine 或Vincenty 公式的实现（如此处所示的Haversine：Haversine Formula in Python (Bearing and Distance between two GPS points)），或者我们可以使用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（文森蒂）：

      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
  

2。添加海拔高度：

如前所述，上述每个计算都会在 2 个点之间产生很大的圆距离。该距离也称为“乌鸦飞”，假设“乌鸦”在不改变高度的情况下从 A 点飞到 B 点尽可能笔直。

通过将先前方法之一的结果与点 A 和点之间的高度差 (delta) 相结合，我们可以更好地估计“步行/驾驶”（“当乌鸦走”？？）距离B、在Euclidean Formula里面进行距离计算：

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

---

以前的解决方案容易出错，尤其是点之间的实际距离越长。
出于继续评论的原因，我将其留在这里。

GeoDjango Distance 计算两点之间的二维距离，不考虑海拔差异。
为了进行 3D 计算，我们需要创建一个距离函数，该函数将在计算中考虑海拔差异：

理论：

latitude、longitude 和 altitude 是 Polar coordinates，我们需要将它们转换为 Cartesian coordinates (x, y, z)，以便在它们上应用 Euclidean Formula 并计算它们的 3D 距离。

• 假设：
  polar_point_1 = (long_1, lat_1, alt_1)
  和
  polar_point_2 = (long_2, lat_2, alt_2)

• 利用这个公式将每个点转换成它的笛卡尔等价物：

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

您将分别获得p_1 = (x_1, y_1, z_1) 和p_2 = (x_2, y_2, z_2) 积分。

• 最后用欧几里得公式：

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

Answer 2

转换成笛卡尔坐标后，就可以用 numpy 计算范数了：

np.linalg.norm(point_1 - point_2)

Answer 3

使用geopy，这是最简单完美的解决方案。

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