找到r中某个纬度/经度距离内的位置

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 km.然后我存储一个TRUE或FALSE变量,然后我可以使用它来平均数据(其他变量).从这个方法,我想要一个TRUE或FALSE的矩阵,用于距离lat和lon 500公里范围内的位置.有没有更有效的方法来做到这一点？

Answer 1

时序:

比较@ 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中断内圆圈内没有其他点.

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

时间代码:

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 2

包的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是的长度latgood和longood),即531米的距离,反对259200我以前的方法.