R代码计算多个经纬度对之间的距离并提取最近的对

Posted 2023-03-05

【中文标题】R代码计算多个经纬度对之间的距离并提取最近的对

【英文标题】：R code to calculate distance between multiple latitude-longitude pairs and extract the closest pairs

【发布时间】：2020-09-04 11:22:44

【问题描述】：

我有两个格式如下的 csv 文件

文件1.csv

Sr      Lat,Long
1       52.361176,4.899779
2       52.34061,4.871195
3       52.374749,4.893847
4       52.356624,4.912281
5       52.374026,4.883685
6       52.369956,4.919778
7       52.370895,4.8703
8       52.390454,4.915024
9       52.378576,4.900253
10      52.378372,4.896219
11      52.380056,4.899697
12      52.383744,4.875805
13      52.369981,4.881528
14      52.375954,4.904786
15      52.344417,4.891211
......1000 columns

文件2.csv

neighbourhood   LAT,LONG
Bijlmer-Centrum 52.3135175, 4.9547795
Bijlmer-Oost    52.3179787, 4.9754974
Bos en Lommer   52.3807577, 4.8545966
Buitenveldert - Zuidas  52.3382516, 4.872921499999999
Centrum-Oost    51.208107, 4.4249047
Centrum-West    52.0607927, 4.4832451
De Aker - Nieuw Sloten  52.3447535, 4.811520799999999
De Baarsjes - Oud-West  52.367746, 4.854258
De Pijp - Rivierenbuurt 52.3560276, 4.9021384
........500columns

我想计算两个文件之间的最短距离对（可能按降序排列）。此外，File1 中的每一对都应对应于 File2 中最近的位置，即不应遗漏 File1 中的任何条目。例如，考虑 file1 52.361176,4.899779 中的第一个 lat-long 对，我需要这对与 File2 中的所有其他对的距离，同样对 File1 中的所有其他条目执行此操作。这是我需要使用的公式（在 python 中）

def distance(lat1, lon1, lat2, lon2):
    p = pi/180
    a = 0.5 - cos((lat2-lat1)*p)/2 + cos(lat1*p) * cos(lat2*p) * (1-cos((lon2-lon1)*p))/2
    return 12742 * asin(sqrt(a))

我是 R 新手，因此请本论坛上的专家提供帮助。

编辑：File1 和 File2 包含的条目比这里提到的要多，这只是一个 sn-p。原始文件分别包含超过 1000 和 500 列。

【问题讨论】：

你可以看看geosphere::distHaversine

在过去的几年里，我在 SO 上多次看到过这样的问题。你看过他们的答案吗？

【参考方案1】：

这是一个使用sf 的空间连接... data.table::fread() 用于创建示例数据。

#make spatial objects
sf1 <- file1 %>% sf::st_as_sf( coords = c("Long", "Lat"), crs = 4326 )
sf2 <- file2 %>% sf::st_as_sf( coords = c("LONG", "LAT"), crs = 4326 )

st_join( sf1, sf2, join = st_nearest_feature )
# 
# Simple feature collection with 15 features and 3 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 4.8703 ymin: 52.34061 xmax: 4.919778 ymax: 52.39045
# geographic CRS: WGS 84
# First 10 features:
#    Sr         neighbourhood                  geometry
# 1   1  De Pijp - Rivierenbuurt POINT (4.899779 52.36118)
# 2   2  Buitenveldert-Zuidas POINT (4.871195 52.34061)
# 3   3  De Pijp - Rivierenbuurt POINT (4.893847 52.37475)
# 4   4  De Pijp - Rivierenbuurt POINT (4.912281 52.35662)
# 5   5  De Pijp - Rivierenbuurt POINT (4.883685 52.37403)
# 6   6  De Pijp - Rivierenbuurt POINT (4.919778 52.36996)
# 7   7  De Baarsjes - Oud-West   POINT (4.8703 52.37089)
# 8   8  De Pijp - Rivierenbuurt POINT (4.915024 52.39045)
# 9   9  De Pijp - Rivierenbuurt POINT (4.900253 52.37858)
# 10 10  De Pijp - Rivierenbuurt POINT (4.896219 52.37837)

使用的样本数据

library(sf)
library(data.table)
file1 <- data.table::fread("
Sr   Lat Long
1    52.361176 4.899779
2    52.34061 4.871195
3    52.374749 4.893847
4    52.356624 4.912281
5    52.374026 4.883685
6    52.369956 4.919778
7    52.370895 4.8703
8    52.390454 4.915024
9    52.378576 4.900253
10   52.378372 4.896219
11   52.380056 4.899697
12   52.383744 4.875805
13   52.369981 4.881528
14   52.375954 4.904786
15   52.344417 4.891211")

file2 <- data.table::fread('                            neighbourhood   LAT LONG
"Bijlmer-Centrum" 52.3135175 4.9547795
"Bijlmer-Oost"    52.3179787 4.9754974
"Bos en Lommer"   52.3807577 4.8545966
"Buitenveldert-Zuidas"  52.3382516 4.872921499999999
"Centrum-Oost"    51.208107 4.4249047
"Centrum-West"    52.0607927 4.4832451
"De Aker - Nieuw Sloten"  52.3447535 4.811520799999999
"De Baarsjes - Oud-West"  52.367746 4.854258
"De Pijp - Rivierenbuurt" 52.3560276 4.9021384')

【讨论】：

感谢您的帮助，不过我有几个问题。 1. 除了显示前 10 个特征，我可以对所有配对进行总结吗？ 2. 也可以计算距离吗？ 3. 与半正弦公式相比，这种方法的准确度如何？

你是否运行过所用函数和包的代码和帮助文件？

是的，我做到了，但不知道如何显示所有而不是前 10 个功能

你应该将结果存储到一个对象中...mything <- .....