对角线上非零值的 R 距离矩阵 (rdist.earth )答案

【问题标题】：R distance matrix with non-zero values on diagonal (rdist.earth )对角线上非零值的 R 距离矩阵 (rdist.earth )
【发布时间】：2018-03-12 23:00:38
【问题描述】：

空间距离矩阵对角线上的元素应该为零，因为它们代表每个位置与其自身之间的距离。但是来自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。首先，为什么第二列的条目显示的数字比其他两个多？为什么第二条对角线上的条目不像其他条目那样从零到八位小数？差异似乎不足以归因于浮点舍入误差。

现在将rdist.earth() 的输出与另一个包geosphere 和函数distGeo() 的输出进行比较，后者计算两点之间的距离（不是完整的距离矩阵）。在这里，我们计算每个点与其自身之间的距离。输出向量单位为米：

> library(geosphere)
> distmat2 <- distGeo(LLdat,LLdat)
> distmat2
[1] 0 0 0

因此，对于distGeo()，所有三个距离度量都一致并且适当地为零。

我有什么遗漏吗？或者这是否表明rdist.earth() 有问题？

【问题讨论】：

标签： r geographic-distance

【解决方案1】：

不幸的是，这是一个舍入错误。

如果您查看源代码，您可以复制问题：

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

【讨论】：

【解决方案2】：

正如@thc 所解释的，这确实是一个数字问题，显然与formula choice 有关。特别要注意的是，在使用acos 之前，所有值都非常接近于 1。acos 在 x 处的导数是 -(1-x^2)^(-1/2)，随着 x 到 -Inf 发散1，所以公式敏感也就不足为奇了。

为了解决这个问题，您可以在 wikipedia 页面中实施其他提议且更稳定的解决方案之一，使用 geosphere 因为他们 seem to have 更谨慎的实施，或者当然您可以设置 @987654325 @。但是，后一种选择可能并不可取，因为当点在空间中非常接近时，这些数字问题也可能存在于非对角项中。

【讨论】：