函数返回关闭0.001的结果

Question

TL；DR：返回的浮点值被0.0001关闭。是我的函数实现导致了这一点，还是仅仅是浮动值和计算机科学的奇迹，因为不是所有的值都可以相应地表示出来？我不认为这是后者，所以我想知道是否有人可以发现问题，给一个建议，或让我试着调试的东西。你知道吗

我写了两个函数来确定地球表面多边形的质心。你知道吗

第一个函数计算两个顶点之间的方向角。它会这样做两次，一次是UR到LL，然后一次是LR到UL。你知道吗

第二个函数计算给定的两个先前计算的方向角的交点以及每个计算的起点。所以你和它的方位和左和它的方位。你知道吗

我的结果在经度上有点偏差。例如，我向interest()提供：

x1 = 120.6272
y1 = 36.9647
x2 = 120.6128
y2 = 36.9508
b1 = 265.02976409662136
b2 = 354.4744218802173

我的结果是[36.96359071593996, 120.61125138605212]。纬度已更正，但经度应为120.61135138605212。你知道吗

方程式由here转录而来。有一个计算器，您可以在“两条路径的交点”部分下提供十进制度数和方向角，以查看预期结果。当从DMS转换为DD时，该结果与ArcGIS中由同一多边形执行的质心计算结果一致。我知道一个gisstackexchange，但是我没有把重点放在这个问题的GIS部分上。我的问题在于：

问题：

在我当前的实现中，是存储函数[即cos = math.cos]，还是本机使用平方[即x ** 2与math.pow()]还是其他一些我不知道的事情导致浮点的小数位数有这么小的偏差（0.0001对我来说很重要）？

承载功能

def bearing(x1, y1, x2, y2):

    cos = math.cos
    sin = math.sin
    atan2 = math.atan2
    rad = math.radians
    deg = math.degrees

    rx1 = rad(x1)
    ry1 = rad(y1)
    rx2 = rad(x2)
    ry2 = rad(y2)

    b = atan2(sin(rx2-rx1) * cos(ry2),
              cos(ry1) * sin(ry2) - (sin(ry1) * cos(ry2) * cos(rx2-rx1)))

    return ((deg(b) + 360) % 360) # normalize to -180...+180

相交函数

def intersect(x1, y1, b1, x2, y2, b2):

    cos = math.cos
    sin = math.sin
    atan2 = math.atan2
    asin = math.asin
    acos = math.acos

    sqrt = math.sqrt

    rad = math.radians
    deg = math.degrees
    pi = math.pi

    rx1 = rad(x1)
    ry1 = rad(y1)
    rx2 = rad(x2)
    ry2 = rad(y2)
    rb1 = rad(b1)
    rb2 = rad(b2)

    dx = rx2 - rx1
    dy = ry2 - ry1

    # Angular dist x1,y1 to x2,y2
    a12 = 2 * asin(sqrt((sin(dy/2)**2) +
                        (cos(ry1) * cos(ry2) * (sin(dx/2)**2))))

    # Initial bearing from x1,y1 to x2,y2
    bi = acos((sin(ry2) - (sin(ry1) * cos(a12))) / (sin(a12) * cos(ry1)))

    if math.isnan(bi):
        bi = 0

    # Final bearing
    bf = acos((sin(ry1) - (sin(ry2) * cos(a12))) / (sin(a12) * cos(ry2)))

    if (dx) > 0:
        # Bearing from x1,y1 to x2,y2
        b12 = bi
        # Bearing from x2,y2 to x1,y1
        b21 = (2 * pi) - bf
    else:
        b12 = (2 * pi) - bi
        b21 = bf

    # Angle x2,y2 -> x1,y1 -> x3,y3
    n1 = rb1 - b12

    # Angle x1,y1 -> x2,y2 -> x3,y3
    n2 = b21 - rb2

    n3 = acos(-(cos(n1) * cos(n2)) +
              (sin(n1) * sin(n2) * cos(a12)))

    # Angular dist x1,y1 to x3,y3
    a13 = atan2(sin(a12) * sin(n1) * sin(n2),
                cos(n2) + (cos(n1) * cos(n3)))

    # Latitude
    ry3 = asin((sin(ry1) * cos(a13)) +
               (cos(ry1) * sin(a13) * cos(rb1)))

    # Longitude Delta from x1 to x3
    dx13 = atan2(sin(rb1) * sin(a13) * cos(ry1),
                 cos(a13) - (sin(ry1) * sin(ry3)))

    # Longitude
    rx3 = rx1 + dx13

    return [deg(ry3), (deg(rx3) + 540) % 360 - 180]

Answer 1

在每个函数中，最后要做的就是以有损的方式将弧度转换为度。 请尝试相同的代码，但不要在每个返回行中进行舍入：例如

return deg(b)

而不是

return ((deg(b) + 360) % 360)