在 r 中的特定 lat/lon 距离内查找位置

Question

我有一个网格数据集，数据位于以下位置：

lon <- seq(-179.75,179.75, by = 0.5)
lat <- seq(-89.75,89.75, by = 0.5)

我想找到位置 500 公里以内的所有数据点：

mylat <- 47.9625
mylon <- -87.0431

我打算在R中使用geosphere包，但是我目前写的方法似乎不是很有效：

require(geosphere)
dd2 <- array(dim = c(length(lon),length(lat)))
for(i in 1:length(lon)){
  for(ii in 1:length(lat)){
    clon <- lon[i]
    clat <- lat[ii]
    dd <- as.numeric(distm(c(mylon, mylat), c(clon, clat), fun = distHaversine))
    dd2[i,ii] <- dd <= 500000
  }
}

在这里，我循环遍历数据中的每个网格，并查找距离是否小于 500 公里。然后我用 TRUE 或 FALSE 存储一个变量，然后我可以用它来平均数据（其他变量）。通过这种方法，我想要一个矩阵，其中 TRUE 或 FALSE 用于距离显示的纬度和经度 500 公里以内的位置。有没有更有效的方法来做到这一点？

Answer 1

geosphere包的dist*函数是向量化的，所以你只需要更好地准备你的输入。试试这个：

#prepare a matrix with coordinates of every position
allCoords<-cbind(lon,rep(lat,each=length(lon)))
#call the dist function and put the result in a matrix
res<-matrix(distm(cbind(mylon,mylat),allCoords,fun=distHaversine)<=500000,nrow=length(lon))
#check the result
identical(res,dd2)
#[1] TRUE

正如@Floo0 的回答所示，有很多不必要的计算。我们可以采用另一种策略：我们首先确定可以接近阈值的lon和lat范围，然后我们只使用它们来计算距离：

#initialize the return
res<-matrix(FALSE,nrow=length(lon),ncol=length(lat))
#we find the possible values of longitude that can be closer than 500000
#How? We calculate the distances between us and points with our same lon 
longood<-which(distm(c(mylon,mylat),cbind(lon,mylat))<=500000)
#Same for latitude
latgood<-which(distm(c(mylon,mylat),cbind(mylon,lat))<=500000)
#we build the matrix with only those values to exploit the vectorized
#nature of distm
allCoords<-cbind(lon[longood],rep(lat[latgood],each=length(longood)))
res[longood,latgood]<-distm(c(mylon,mylat),allCoords)<=500000

这样，你只计算lg+ln+lg*ln（lg和ln是latgood和longood的长度），即531个距离，反对我以前的方法的 259200。

Answer 2

时间：

比较@nicola 和我的版本给出：

Unit: milliseconds

               min         lq      mean     median         uq       max neval
nicola1 184.217002 219.924647 297.60867 299.181854 322.635960 898.52393   100
floo01   61.341560  72.063197  97.20617  80.247810  93.292233 286.99343   100
nicola2   3.992343   4.485847   5.44909   4.870101   5.371644  27.25858   100

我最初的解决方案：（恕我直言，nicola 的第二个版本更干净、更快。）

您可以执行以下操作（下面有解释）

require(geosphere)
my_coord <- c(mylon, mylat)
dd2 <- matrix(FALSE, nrow=length(lon), ncol=length(lat))
outer_loop_state <- 0
for(i in 1:length(lon)){
    coods <- cbind(lon[i], lat)
    dd <- as.numeric(distHaversine(my_coord, coods))
    dd2[i, ] <- dd <= 500000
    if(any(dd2[i, ])){
      outer_loop_state <- 1
    } else {
      if(outer_loop_state == 1){
        break
      }
    }
  }

解释：

对于循环，我应用以下逻辑：

outer_loop_state 初始化为 0。如果找到圆内至少有一个光栅点的行，outer_loop_state 设置为 1。一旦圆内不再有点给定行 i 中断。

@nicola 版本中的 distm 调用基本上没有这个技巧。所以它计算所有行。

时间代码：

microbenchmark::microbenchmark(
  {allCoords<-cbind(lon,rep(lat,each=length(lon)))
  res<-matrix(distm(cbind(mylon,mylat),allCoords,fun=distHaversine)<=500000,nrow=length(lon))},
  {my_coord <- c(mylon, mylat)
  dd2 <- matrix(FALSE, nrow=length(lon), ncol=length(lat))
  outer_loop_state <- 0
  for(i in 1:length(lon)){
    coods <- cbind(lon[i], lat)
    dd <- as.numeric(distHaversine(my_coord, coods))
    dd2[i, ] <- dd <= 500000
    if(any(dd2[i, ])){
      outer_loop_state <- 1
    } else {
      if(outer_loop_state == 1){
        break
      }
    }
  }},
  {#intitialize the return
    res<-matrix(FALSE,nrow=length(lon),ncol=length(lat))
    #we find the possible value of longitude that can be closer than 500000
    #How? We calculate the distance between us and points with our same lat 
    longood<-which(distm(c(mylon,mylat),cbind(lon,mylat))<500000)
    #Same for latitude
    latgood<-which(distm(c(mylon,mylat),cbind(mylon,lat))<500000)
    #we build the matrix with only those values to exploit the vectorized
    #nature of distm
    allCoords<-cbind(lon[longood],rep(lat[latgood],each=length(longood)))
    res[longood,latgood]<-distm(c(mylon,mylat),allCoords)<=500000}
)

Answer 3

直接用hutils::haversine_distance(lat, lon, mylat, mylon) < 500就可以了。

如果假设这些点是给定的 lat 和 lon 的交叉连接，请先使用交叉连接来获取它们：

library(data.table)
library(hutils)

lon <- seq(-179.75,179.75, by = 0.5)
lat <- seq(-89.75,89.75, by = 0.5)

mylat <- 47.9625
mylon <- -87.0431

Points <- CJ(lon = lon,
             lat = lat)
Points[, dist := haversine_distance(lat, lon, mylat, mylon)]
Points[, sum(dist < 500)]
#> [1] 379

^{由 reprex package (v0.3.0)}

于 2019-10-24 创建

它通过速度和稳健性改进了现有答案。特别是，它不依赖于数据的网格化性质，并且可以处理长坐标向量。以下是100,000积分的时间

# A tibble: 2 x 14
  expression         min        mean      median         max `itr/sec`  mem_alloc  n_gc n_itr  total_time
  <chr>         <bch:tm>    <bch:tm>    <bch:tm>    <bch:tm>     <dbl>  <bch:byt> <dbl> <int>    <bch:tm>
1 nicola2    39891.120ms 39891.120ms 39891.120ms 39891.120ms    0.0251 8808.632MB     0     1 39891.120ms
2 hutils        15.492ms    15.591ms    15.578ms    15.728ms   64.1       5.722MB     0    33   514.497ms

Answer 4

我在下面添加了一个使用 spatialrisk 包的解决方案。这个包中的关键函数是用 C++ (Rcpp) 编写的，因此速度非常快。

首先加载数据：

mylat <- 47.9625
mylon <- -87.0431

lon <- seq(-179.75,179.75, by = 0.5)
lat <- seq(-89.75,89.75, by = 0.5)
df <- expand.grid(lon = lon, lat = lat)

函数spatialrisk::points_in_circle() 计算距离中心点半径范围内的观测值。请注意，距离是使用 Haversine 公式计算的。

与@Hugh 版本相比，空间风险方法的时间安排：

spatialrisk::points_in_circle(df, mylon, mylat, radius = 5e5)

Unit: milliseconds
       expr       min        lq      mean    median        uq       max neval cld 
spatialrisk  3.071897  3.366256  5.224479  4.068124  4.809626  17.24378   100   a 
     hutils 17.507311 20.788525 29.470707 25.061943 31.066139 268.29375   100   b

结果可以很容易地转换为矩阵。

看看@philcolbourn 关于如何测试一个点是否在圆内的出色回答。参见：

在 r 中的特定 lat/lon 距离内查找位置

find locations within certain lat/lon distance in r

r

geosphere