Haversine 公式和 Python 3 - 数学域错误

Question

我正在尝试实现 Haversine 的公式来确定给定位置的纬度和经度是否在指定半径内。我在这里使用了一个详细的公式：
Calculate distance between two latitude-longitude points? (Haversine formula)

我遇到了以下可重现输入的数学域错误，它不会一直发生，但经常足以让我认为我编写了不正确的代码：

from math import atan2, sqrt, sin, cos

# All long / lat values are in radians and of type float
centerLongitude = -0.0391412861306467
centerLatitude = 0.9334153362515779

inputLatitudeValue = -0.6096173085842176
inputLongitudeValue = 2.4190393564390438

longitudeDelta = inputLongitudeValue - centerLongitude # 2.4581806425696904
latitudeDelta = inputLatitudeValue - centerLatitude # -1.5430326448357956

a = (sin(latitudeDelta / 2) ** 2 + cos(centerLatitude) * cos(centerLongitude) 
     * sin(longitudeDelta / 2) ** 2)
# a = 1.0139858858386017

c = 2 * atan2(sqrt(a), sqrt(1 - a)) # Error occurs on this line

# Check whether distance is within our specified radius below

Answer 1

您不能在负数上使用sqrt：

>>> sqrt(-1)   
ValueError: math domain error

使用cmath.srt:

>>> import cmath
>>> cmath.sqrt(-1)
1j

在你的情况下：

>>> a = 1.0139858858386017
>>> sqrt(1-a)
ValueError: math domain error

Answer 2

笼统地说，简单地说……变量a必须受到保护。绝对不能大于 1.0，也不能小于 0.0，通常情况下，对于干净的数据，它应该在这个范围内。

问题在于浮点运算的常见流行计算机实现如何进行近似和舍入，以及这些结果有时如何超出内置数学函数的范围或域。运算的组合导致后续数学函数超出域的数字并不特别常见，但在确实发生的情况下，取决于算法，它始终是可重现的，并且需要考虑并减轻，超出了正常的理论算法直觉。我们在计算机中编码的内容取决于软件和硬件的理论概念是如何实现的。在纸上，用铅笔，我们仍然必须有关于何时以及如何舍入浮点数学结果的指南。计算机实现没有什么不同，但有时我们很高兴地没有意识到引擎盖下发生的这些事情。流行的计算机实现不仅需要知道要以什么精度进行舍入，而且还要处理在机器级别实际计算的二进制数表示的转换中如何以及何时进行近似和舍入。

关于半正弦公式，我在说什么a？

如本版本代码所示（供参考）：

import math

a = math.sin(dlat / 2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2)**2
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))

在上面的示例中，a 没有得到适当的保护。这是一个潜伏的问题，等待下一行的 c 计算崩溃，在某些情况下，特别是对于 hasrsine 公式。

如果某些 latlon 数据组合导致

a = 1.00000000000000000000000000000000000001

以下 c 计算会导致错误。

如果某些 latlon 数据组合导致

a = -0.00000000000000000000000000000000000000000000001

下面的c计算会出错。

这是您的语言/平台的浮点数学实现及其舍入近似值的方法，可能导致罕见但实际且始终可重复的稍微超出域 在未受保护的正弦编码中导致错误的条件。

几年前，我对 0 和 1 以及 179 和 180 之间的相对角距离进行了为期三天的蛮力测试，步长值非常小。球体的半径为 1，单位球体，因此半径值无关紧要。当然，以角距离以外的任何单位计算球体表面上的近似距离都需要在其单位中包括球体的半径。但是我正在测试haversine逻辑实现本身，半径为1消除了复杂性。当相对角距离为 0 -1 或 179-180 时，这些是半正弦可能遇到困难的情况，使用流行的浮点实现来实现，这些实现涉及在低系统级别转换为二进制，如果 a 不受保护。理论上，Haversine 应该适用于小角距离，但 FPA（浮点算法）的机器或软件实现并不总是与球面几何理论的理想精确配合。在我进行了 3 天的蛮力测试后，记录了数千个 latlon 组合，这些组合会使未受保护的流行发布的半正弦公式崩溃，因为 a 没有受到保护。您必须保护 a。如果它高于 1.0 或低于 0.0，即使是最轻微的一点，只需测试该条件并将其推回范围即可。很简单。

我保护了 -0.000000000000000000002 的 a对这些数据采取保护措施是必要的。

我保护了一个 1.000000000000000000002 的 a需要对该数据采取行动。

在 c 计算之前是简单的 2 行。你可以把它全部挤在一条线上。代码的保护行在 a 计算之后，c 计算之前。保护您的a。

您是否会因为这些轻微的轻推而失去精确度？不超过使用该数据的浮点数学已经引入的近似值和舍入。它崩溃的数据不应该崩溃纯理论算法，一个没有 FPA 罕见问题的算法。只需保护 a，这应该可以减轻这些错误，使用这种形式的 hasrsine。 hasrsine 有其他替代品，但如果人们了解它非常适合的地方，haversine 完全适合许多用途。我将它用于天球计算，其中地球的椭球形状与任何事物无关。 Haversine 简单快速。但请记住保护您的a。