21

次の場所でデータを利用できるグリッド データセットがあります。

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

その場所から 500 km 以内にあるすべてのデータ ポイントを検索したいと思います。

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 のいずれかで変数を保存します。これを使用して、データ (他の変数) を平均化できます。このメソッドから、表示されている緯度と経度から 500 km 以内の場所の TRUE または FALSE を含むマトリックスが必要です。これを行うためのより効率的な方法はありますか?

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

私の元の解決策: (IMHO nicola の 2 番目のバージョンは、はるかにクリーンで高速です。)

次のことができます（以下の説明）

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 で初期化されます。円の内側に少なくとも 1 つのラスター ポイントを持つ行が見つかった場合は、1 に設定されます。指定された行ブレークouter_loop_stateの円内にポイントがなくなると、 .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}
)
于 2016-08-30T09:27:32.967 に答える
7

The dist* functions of the geosphere package are vectorized, so you only need to prepare better your inputs. Try this:

#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

As the @Floo0 answer showed, there is a lot of unnecessary calculations. We can follow another strategy: we first determine the lon and lat range that can be closer than the threshold and then we use only them to calculate the distance:

#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

In this way, you calculate just lg+ln+lg*ln (lg and ln are the length of latgood and longood), i.e. 531 distances, opposed to the 259200 with my previous method.

于 2016-08-30T08:56:20.143 に答える