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

How to measure distance between points in separate data frames?

我创建了 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)

但我得到了一个矩阵,其中计算了来自两个几何列的每一对的距离:

           [,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 对我来说不是一个选项,因为我不想按几何连接(两个数据框中的几何图形完全不同)

正如@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()

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

超快速向量化计算

此方法适用于:

  1. 将(经度、纬度)坐标投影到与您感兴趣的区域等距的相关坐标系。 (等距坐标系保留点之间的距离测量值,因此您可以只使用基本几何来计算距离)。
  2. 将几何图形转换为具有 X 和 Y 列的基本 R 矩阵。
  3. 最后,简单地使用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))
})