найти места в пределах определенного расстояния в г

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

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, которую затем могу использовать для усреднения данных (другая переменная). Из этого метода я хочу матрицу с ИСТИНА или ЛОЖЬ для мест в пределах 500 км от показанных значений широты и долготы. Есть ли более эффективный способ сделать это?

 Emma Tebbs30 авг. 2016 г., 11:05
@tonytonov это все еще медленно, чтобы то, что мне нужно. Эта функция будет вводиться в более сложный код, поэтому должна быть быстрее
 loki30 авг. 2016 г., 10:34
Вы смотрели наrgeos пакет? ФункцияgDistance() мог решить эту проблему, вернув матрицу с расстояниями. После этого будет только матричная алгебра, чтобы узнать расстояние<=500000
 tonytonov30 авг. 2016 г., 10:41
Пытатьсяouter(lon, lat, Vectorize(function(x, y) distm(c(mylon, mylat), c(x, y), fun = distHaversine)[1, 1] < 500000))Это достаточно быстро?

Ответы на вопрос(1)

Решение Вопроса

Тайминги:

Сравнение @ Николы и моей версии дает:

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 Call в @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}
)
 Rentrop30 авг. 2016 г., 14:25
@nicola: Ваш ответ - еще одно огромное улучшение! Спасибо за дальнейшее улучшение.
 nicola30 авг. 2016 г., 14:05
Хороший ответ. Вдохновленный вашим трюком, я сделал правку, чтобы использовать как векторизацию, так и тот факт, что некоторые точки явно находятся дальше порога (как вы показали). Это должно значительно ускорить процесс.
 nicola30 авг. 2016 г., 14:33
@ Спасибо за идею. Не могли бы вы отредактировать свой ответ, чтобы сравнить новое предлагаемое решение? Tx.

Ваш ответ на вопрос