Найти местоположения в пределах определенного расстояния lat/lon в r

У меня есть набор данных с сеткой с сеткой, данные доступны в следующих местах:

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, но метод, который я сейчас написал, кажется не очень эффективным:

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 км от показанных лат и lon. Есть ли более эффективный метод для этого?

Ответ 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

Мое оригинальное решение: (вторая версия IMHO 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 сломаться.

Вызов 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}
)

Ответ 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 с моим предыдущим методом.