21

我有一个网格数据集，数据可在以下位置获得：

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 公里以内的位置。有没有更有效的方法来做到这一点？

4

4 回答 4

10

时间：

比较@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

我原来的解决方案：（恕我直言，尼古拉的第二个版本更干净，更快。）

您可以执行以下操作（以下说明）

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}
)
于 2016-08-30T09:27:32.967 回答
7

该包的dist*功能geosphere是矢量化的，因此您只需要准备更好的输入。尝试这个：

#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的长度），即531个距离，与我之前的方法的259200相反。latgoodlongood

于 2016-08-30T08:56:20.143 回答
1

我在下面添加了一个使用 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 方法的时间安排：

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 关于如何测试一个点是否在一个圆圈内的出色答案。请参阅：https ://stackoverflow.com/a/7227057/5440749

于 2019-10-24T18:35:26.700 回答
0

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

如果假设这些点是给定latand的交叉连接，则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 包（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
于 2019-02-09T15:30:09.807 回答