如何测量单独数据框中点之间的距离？

Posted 2023-03-23

【中文标题】如何测量单独数据框中点之间的距离？

【英文标题】：How to measure distance between points in separate data frames?

【发布时间】：2021-04-06 05:57:48

【问题描述】：

我创建了 2 个带有 geom 列（POINT 类型）的数据框。现在我想计算每对点之间的距离，例如从第一个 df 中的第一行指向第二个 df 中的第一行等。这是我的数据框：

df1 <- table %>%
  st_as_sf(coords = c("lonCust","latCust"), crs = 4326)

df2 <- table %>%
  st_as_sf(coords = c("lonApp","latApp"), crs = 4326)

我用st_distance:

distance <- st_distance(df1$geometry,df2$geometry)

但我得到了一个矩阵，其中计算了两个 geom 列中每一对的距离：

           [,1]      [,2]        [,3]         [,4]        [,5]  ...
[1,]   139.7924 7735.5718 15225.02995   558.104089  1016.58121
[2,]  8503.0544  755.2915  8764.75396  7957.289600  8788.02800
[3,] 15306.5855 9336.9008    18.96914 14876.589918 15929.51643
[4,]   548.3045 7232.0164 14898.70637     8.094068  1078.38236
[5,]   911.5635 8084.3086 15993.36365  1127.730022    46.97799
.
.

我希望在一列中计算距离，仅在相应的几何行之间：

           [,1]     
[1,]   139.7924 
[2,]  8503.0544
[3,] 15306.5855 
[4,]   548.3045
[5,]   911.5635
.
.

我读过geosphere 包，但sf 有非常好的st_distance 功能来测量距离，我想用它。最重要的是，我需要先加入这些数据框吗？来自dplyr 的简单inner_join 不允许加入两个空间数据帧，另一方面st_join 对我来说不是一个选项，因为我不想通过几何加入（两个数据帧中的几何是完全不同）

【问题讨论】：

你已经接近了。 st_distance(df1$geometry, df2$geometry, by_element = TRUE)

谢谢！顺便说一句，我有大约 25000 行，计算速度很慢，我认为最好在服务器端计算距离

library(geodist) 对于这类事情可能更快，它可以直接对data.frames 进行操作。查看geodist() 和/或geodist_vec()，并设置paired = TRUE。

【参考方案1】：

正如@mrhellmann 提到的，您可以添加by_element=T，这应该可以。如果速度仍然是个问题，我建议使用geosphere 包中的DistGeo()。但请务必查看文档以了解您的数据是否适合此功能。

library(geosphere)
library(tidyverse)
library(sf)

df1 <- table %>%
  st_as_sf(coords = c("lonCust","latCust"), crs = 4326)

doParallel::registerDoParallel()
df_crs4326 <- df1 %>%
  group_by(your_id_here) %>% 
  mutate(
    lonCust = map(geometry, 2) %>% unlist(),
    latCust= map(geometry, 1) %>% unlist(),
    # geometry_2 = st_as_sfc(coords = c("lonApp","latApp"), crs = 4326)
    ) %>%
  mutate(
    distance_to_next = distGeo(c(lonCust, latCust), c(lonApp, latApp)) %>% set_units(m),
    # distance_2 = st_distance(geometry, geometry_2, by_element = TRUE)
    ) %>%
    ungroup()

请注意，如果没有对可重现数据进行测试，我不确定被注释掉的部分是否有效。

【参考方案2】：

超快速向量化计算

此方法的工作原理：

将（经度、纬度）坐标投影到与您感兴趣的区域等距的相关坐标系。 （等距坐标系保留点之间的距离测量值，因此您可以只使用基本几何来计算距离）。 将几何图形转换为具有 X 和 Y 列的 Base R 矩阵。 最后，只需使用Pythagoras's theorem 计算点对之间的距离。

首先获取坐标参考系 (CRS)

为此，您需要一个等距的 CRS。这意味着，在感兴趣的区域内，任何距离计算都会被保留。

假设您有兴趣计算横跨美国的距离，您可以使用EPSG:102005。有关模式详细信息，请参阅此 GIS 答案。这里 CRS 的选择很关键，所以请确保你做对了，否则答案将是无稽之谈。

应用于您的示例

crs.source = 4326
crs.dest = st_crs("+proj=eqdc +lat_0=39 +lon_0=-96 +lat_1=33 +lat_2=45 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")

# coords1 and coords2 are matrixes with columns X and Y and rows of points in the `crs.dest` coordinate system.
coords1 <- table %>%
  st_as_sf(coords = c("lonCust","latCust"), crs = crs.source) %>%
  st_transform(crs.dest) %>%
  st_coordinates()
  
coords2 <- table %>%
  st_as_sf(coords = c("lonApp","latApp"), crs = crs.source) %>%
  st_transform(crs.dest) %>%
  st_coordinates()

# This is a vectorised computation, and so should be instant for a mere 25,000 rows :-)
table$distances = local(
  x_diff = coords1[, 'X'] - coords2[, 'X']
  y_diff = coords1[, 'Y'] - coords2[, 'Y']
  return(sqrt(x^2 + y^2))
)