如何计算R中低于某个阈值的2个坐标之间的距离？答案

【问题标题】：How to calculate distance between 2 coordinates below a certain threshold in R?如何计算R中低于某个阈值的2个坐标之间的距离？
【发布时间】：2016-08-09 17:43:19
【问题描述】：

我有 44,000 个美国邮政编码，它对应于 R 中的质心纬度/经度。这是来自 R 中的“邮政编码”包。我需要计算每个邮政编码之间的距离并保持小于 5 英里的距离。问题是要计算邮政编码之间的所有距离，我必须创建一个大小为 44,000x44,0000 的向量，由于空间问题，我不能。

我检查了 R 中的帖子，最接近我的要求是吐出 2 个具有 lat/long 的数据集之间的最小距离

DB1 <- data.frame(location_id=1:7000,LATITUDE=runif(7000,min = -90,max = 90),LONGITUDE=runif(7000,min = -180,max = 180))
DB2 <- data.frame(location_id=7001:12000,LATITUDE=runif(5000,min = -90,max = 90),LONGITUDE=runif(5000,min = -180,max = 180))

DistFun <- function(ID){
  TMP <- DB1[DB1$location_id==ID,]
  TMP1 <- distGeo(TMP[,3:2],DB2[,3:2])
  TMP2 <- data.frame(DB1ID=ID,DB2ID=DB2[which.min(TMP1),1],DistanceBetween=min(TMP1)      ) 
  print(ID)
  return(TMP2)
}

DistanceMatrix <- rbind_all(lapply(DB1$location_id, DistFun))

即使我们可以修改上面的代码以包含所有

是否有一种有效的方法可以得出所有邮政编码组合，这些组合的质心彼此之间

【问题讨论】：

我怀疑您可以跳过对数值大于大约 500 的某个数字的组合的任何距离检查。这将使组合负载更易于管理。
如何过滤掉不在以点 1 为中心的长度为 10 英里的正方形内的对，然后仅计算位于其中的对。 data.tables 应该可以帮助您非常有效地进行过滤。
您想如何存储距离？

标签： r algorithm performance geolocation distance

【解决方案1】：

一次生成整个距离矩阵将非常消耗 RAM，循环遍历唯一邮政编码的每个组合 - 非常耗时。让我们找到一些妥协。

我建议将 zipcode data.frame 分块（例如）100 行（借助 bit 包中的 chunk 函数），然后计算 44336 和 100 点之间的距离，根据目标距离阈值，然后移动到下一个数据块。在我的示例中，我将 zipcode 数据转换为 data.table 以提高速度并节省 RAM。

library(zipcode)
library(data.table)
library(magrittr)
library(geosphere)

data(zipcode)

setDT(zipcode)
zipcode[, dum := NA] # we'll need it for full outer join

仅供参考 - 这是 RAM 中每条数据的大致大小。

merge(zipcode, zipcode[1:100], by = "dum", allow.cartesian = T) %>% 
  object.size() %>% print(unit = "Mb")
# 358.2 Mb

代码本身。

lapply(bit::chunk(1, nrow(zipcode), 1e2), function(ridx) {
  merge(zipcode, zipcode[ridx[1]:ridx[2]], by = "dum", allow.cartesian = T)[
    , dist := distGeo(matrix(c(longitude.x, latitude.x), ncol = 2), 
                      matrix(c(longitude.y, latitude.y), ncol = 2))/1609.34 # meters to miles
    ][dist <= 5 # necessary distance treshold
      ][, dum := NULL]
  }) %>% rbindlist -> zip_nearby_dt

zip_nearby_dt # not the whole! for first 10 chunks only

       zip.x          city.x state.x latitude.x longitude.x zip.y     city.y state.y latitude.y longitude.y     dist
    1: 00210      Portsmouth      NH   43.00590   -71.01320 00210 Portsmouth      NH   43.00590   -71.01320 0.000000
    2: 00210      Portsmouth      NH   43.00590   -71.01320 00211 Portsmouth      NH   43.00590   -71.01320 0.000000
    3: 00210      Portsmouth      NH   43.00590   -71.01320 00212 Portsmouth      NH   43.00590   -71.01320 0.000000
    4: 00210      Portsmouth      NH   43.00590   -71.01320 00213 Portsmouth      NH   43.00590   -71.01320 0.000000
    5: 00210      Portsmouth      NH   43.00590   -71.01320 00214 Portsmouth      NH   43.00590   -71.01320 0.000000
---                                                                                                              
15252: 02906      Providence      RI   41.83635   -71.39427 02771    Seekonk      MA   41.84345   -71.32343 3.688747
15253: 02912      Providence      RI   41.82674   -71.39770 02771    Seekonk      MA   41.84345   -71.32343 4.003095
15254: 02914 East Providence      RI   41.81240   -71.36834 02771    Seekonk      MA   41.84345   -71.32343 3.156966
15255: 02916         Rumford      RI   41.84325   -71.35391 02769   Rehoboth      MA   41.83507   -71.26115 4.820599
15256: 02916         Rumford      RI   41.84325   -71.35391 02771    Seekonk      MA   41.84345   -71.32343 1.573050

在我的机器上处理 10 个块需要 1.7 分钟，所以整个处理可能需要 70-80 分钟，不算快，但可能令人满意。我们可以根据可用的 RAM 容量将块大小增加到 200 或 300 行，这将分别缩短处理时间 2 或 3 倍。

此解决方案的缺点是生成的 data.table 包含“重复”行 - 我的意思是从 A 点到 B 点以及从 B 到 A 点都有距离。这可能需要一些额外的过滤。

【讨论】：

【解决方案2】：

我猜最有效的算法会首先将空间位置转换为树状数据结构。但是，您不需要明确地执行此操作，如果您有一个算法可以 1）将纬度/经度分箱到空间索引，2）告诉您该索引的邻居，那么您可以使用它来过滤您的方形数据。（这会比构建树效率低，但可能更容易实现。）

geohash就是这样一个算法。它将连续的纬度/经度变成二维箱。有一个（相当新的）包providing geohash in R。下面是一个关于如何使用它来解决这个问题的想法：

首先，用 geohash 做一些初步校准：

将 lat/long 转换为 bin 精度 p 的哈希（比如）
评估哈希是否以与您感兴趣的距离相似的精度进行校准（例如，相邻质心之间 3-7 英里），如果没有，则返回 1 并进行调整精度p

这会产生 zipcode-哈希值 关系。

然后，为每个（唯一）哈希值计算距离

确定其（8，bc 哈希形成二维网格）最近邻，因此选择 9 个哈希值
计算 9 个散列中所有 zip 之间的成对距离（例如问题中使用 distGeo）
返回哈希值的所有 zip-zip 成对距离（例如，在矩阵中）

这会产生一个散列值-zip-zip 距离对象 关系

（在步骤 2 中，显然最好只计算每个最近邻对一次。但这可能不是必需的。）

最后，对于每个 zip

使用以上两步（通过hash值作为key）获取zip-zip
拉链的距离对象
将对象过滤到距焦点 zip 的距离（回想一下，它是与焦点 zip 相邻的一组哈希内的所有成对距离）
只保持距离< 5 miles

这会产生一个 zip-5 英里内的拉链 对象。（焦点拉链 5 英里范围内的拉链可以作为列表列（每个元素都是一个列表）存储在数据框中，紧邻焦点拉链列，或者作为单独的列表，其中焦点拉链作为名称）。

【讨论】：

【解决方案3】：

以下是使用spatialrisk 的解决方案。这些函数是用 C++ 编写的，因此速度非常快。在我的机器上大约需要 25 秒。

library(zipcodeR)
library(spatialrisk)
library(dplyr)

# Zip code data
zipcode <- zipcodeR::zip_code_db

# Radius in meters
radius_meters <- 5000

# Find zipcodes within 5000 meters
sel <- tibble(zipcode) %>%
  select(zipcode, lat, lon = lng) %>%
  filter(!is.na(lat), !is.na(lon)) %>%
  mutate(zipcode_within_radius = purrr::map2(lon, lat, ~points_in_circle(zipcode_sel, .x, .y, radius = radius_meters)[-1,])) %>%
  unnest(cols = c(zipcode_within_radius), names_repair = "unique")

【讨论】：