空间距离矩阵对角线上的条目应为零，bc 它们表示每个位置与其自身的距离。但是fieldsR package的rdist.earth()函数有时会在对角线上给我非零：
> # Set number of decimals of output display
> options(digits=8)
> # Some longitude, latitude data
> LLdat
lon lat
1 -105.85878 43.65797
2 -105.81812 43.57009
3 -105.80796 43.57748
>
> # Create distance matrix
> library(fields)
> distmat <- rdist.earth(LLdat,LLdat)
> distmat
1 2 3
1 0.0000000 6.410948951394 6.12184338
2 6.4109490 0.000059058368 0.72150586
3 6.1218434 0.721505863563 0.00000000
在上面的距离矩阵中，对角线上的第二个条目是0.000059058368，以英里为单位(默认单位)，而其他两个条目是0.0000000。首先，为什么第二列的条目显示的数字比其他两列多？为什么第二个对角线上的条目不像其他条目那样是零到小数点后 8 位？这种差异似乎还不够小，不足以归因于浮点舍入误差。
现在将rdist.earth()的输出与不同的包、geosphere和函数distGeo()的输出进行比较，后者计算两点之间的距离(不是全距离矩阵)。在这里，我们计算每个点和自身之间的距离。输出矢量单位以米为单位：
> library(geosphere)
> distmat2 <- distGeo(LLdat,LLdat)
> distmat2
[1] 0 0 0
所以对于distGeo()，所有三个距离度量都是一致的，并且适当地为零。
我错过了什么吗？或者这是否表明rdist.earth()有问题？
不幸的是，这是一个舍入错误。
如果您查看源代码，则可以复制该问题：
x1 <- LLdat
R <- 3963.34
coslat1 <- cos((x1[, 2] * pi)/180)
sinlat1 <- sin((x1[, 2] * pi)/180)
coslon1 <- cos((x1[, 1] * pi)/180)
sinlon1 <- sin((x1[, 1] * pi)/180)
pp <- cbind(coslat1 * coslon1, coslat1 * sinlon1, sinlat1) %*%
t(cbind(coslat1 * coslon1, coslat1 * sinlon1, sinlat1))
return_val = (R * acos(ifelse(abs(pp) > 1, 1 * sign(pp), pp)))
该函数首先计算一个中间矩阵 pp：
print (pp)
[,1] [,2] [,3]
[1,] 1.0000000000 0.9999986917 0.9999988071
[2,] 0.9999986917 1.0000000000 0.9999999834
[3,] 0.9999988071 0.9999999834 1.0000000000
似乎对角线都是一样的。 然而：
print(pp, digits=22)
[,1] [,2] [,3]
[1,] 1.0000000000000000000000 0.9999986917465573110775 0.9999988070789928018556
[2,] 0.9999986917465573110775 0.9999999999999998889777 0.9999999834298258782894
[3,] 0.9999988070789928018556 0.9999999834298258782894 1.0000000000000000000000
> acos(0.9999999999999998889777) * R
[1] 5.905836821e-05
> acos(1.0000000000000000000000) * R
[1] 0
正如@thc所解释的那样，这确实是一个数字问题，显然与公式选择有关。特别要注意的是，在使用acos之前，所有值都非常接近 1。在 x 处的 acos 导数是 -(1-x^2)^(-1/2)，当 x 变为 1 时发散到 -Inf，因此公式敏感也就不足为奇了。
至于处理这个问题，您可以在维基百科页面中实现其他建议的和更稳定的解决方案之一，使用geosphere，因为它们似乎具有更谨慎的实现，或者当然您可以设置diag(M) <- 0。但是，后一种选择可能是不可取的，因为当点在空间中非常接近时，这些数值问题也可能保留在非对角线术语中。