如何在GeoDjango中计算两点之间的三维距离(包括海拔)

5

序言:

这是 SO 上经常出现的问题:

我本来想在 SO 文档中撰写一个例子,但 geodjango 章节从未成立,并且由于文档在 2017 年 8 月 8 日关闭了,因此我将遵循这篇广受赞誉和讨论的元回答的建议并将我的例子写成一个自问自答的帖子。

当然,我也很乐意看到不同的方法!!


问题:

假设有如下模型:

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

我正在尝试计算两个这样的点之间的三维距离:
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,就像高程被忽略一样。
我如何计算我的点之间的三维距离?

3个回答

7

阅读 GEOSGeometry.distance 方法的文档:

返回此几何图形和给定几何图形(另一个 GEOSGeometry 对象)之间最近点的距离。

注意

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

因此,我们需要实现一种方法来计算更精确的二维两点距离,然后再尝试应用这些点之间的高度(Z)差异。

1. 大圆 2D 距离计算 (查看下面的解释中 2022 年更新以获取使用 geopy 更好的方法)

计算球面上 2 点之间距离(如通常简化但通常用于建模的地球),最常用的方式是Haversine 公式

Haversine 公式确定了给定经纬度在球面上的两点间的大圆距离

虽然从大圆距离维基页面上我们可以读到:

尽管此公式对于球体上大多数距离是准确的,但在对蚊子的点(在球体的两端)的特殊(有些不寻常)情况下,也会出现舍入误差。一个对所有距离都准确的公式是等轴椭球体的 Vincenty 公式

我们可以创建我们自己的 Haversine 或 Vincenty 公式实现(如此处所示的 Haversine:Python 中的 Haversine 公式(两个 GPS 位置之间的方位角和距离)),或者我们可以使用包含在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更新:关于使用geopy进行2D距离计算:

从版本1.14.0开始,GeoPy不建议使用Vincenty更改日志中说明:

更改:使用Vincenty现在会发出警告。应该使用大地线代替Vincenty。计划在geopy 2.0中删除Vincenty。(#293)

因此(特别是如果我们将在WGS84椭球体上应用计算),我们应该使用geodesic距离

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. 增加高度因素:

如上所述,以上每个计算都可以得出两点之间的大圆距离。这个距离也被称为“直线距离”,假设“乌鸦”在飞行过程中不改变高度,而是尽可能地从A点飞向B点。

我们可以通过将先前方法之一的结果与点A和点B之间的高度差(delta)结合到欧几里得距离公式中来更好地估计“步行/驾驶”(“乌鸦走路”??)距离:

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

之前的解决方案容易出错,特别是实际距离越长。
我留下它是为了继续评论。

GeoDjango Distance计算两点之间的2D距离,并不考虑高度差异。
为了获得3D计算,我们需要创建一个距离函数,在计算中考虑高度差异:

理论:

纬度经度海拔极坐标系,我们需要将它们转换为笛卡尔坐标系(xyz,以便在它们上应用欧几里得公式并计算它们的3D距离。

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

现在你有两个点,分别为p_1 = (x_1, y_1, z_1)p_2 = (x_2, y_2, z_2)

  • Finally use the Euclidean formula:

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

1
这里计算的距离线可能会穿过地球内部(例如,日本的一个点和美国的另一个点)。这不会导致不准确的答案吗?使用Haversine公式计算大圆距离更加准确(参见en.wikipedia.org/wiki/Haversine_formula)。 - Alan Evangelista
@AlanEvangelista 你说得对。你的评论让我找到了一个更少错误的解决方案。看一下编辑后的答案 :) - John Moutafis
1
sqrt(great_circle(p1, p2).m**2, (p1.z - p2.z)**2) 应该改为 sqrt(great_circle(p1, p2).m**2 + (p1.z - p2.z)**2),对吧?用加号替换逗号? - Krupip
@Krupip 捕捉得好,谢谢! - John Moutafis

1
使用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

geopy的距离函数默认使用WGS-84椭球体,这是全球最准确的。https://geopy.readthedocs.io/en/stable/#module-geopy.distance 我不明白@John Moutafis为什么会对我的答案投反对票。 - ericobi
2
如果您运行此代码并使用除“0”以外的任何值作为高度,将会出现“ValueError:不支持计算具有不同高度的点之间的距离”。 - Dan Taninecz Miller

-1
一旦转换为笛卡尔坐标,你可以使用numpy计算其范数:
np.linalg.norm(point_1 - point_2)

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接