div - sidebarpanel shiny

Encuentre ubicaciones dentro de cierta distancia lat/lon en r

Question

Tengo un conjunto de datos de cuadrícula, con datos disponibles en las siguientes ubicaciones:

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

Me gustaría encontrar todos los puntos de datos que se encuentran a 500 km de la ubicación:

mylat <- 47.9625
mylon <- -87.0431

Mi objetivo es usar el paquete geosphere en R, pero el método que he escrito actualmente no parece muy eficiente:

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
  }
}

Aquí, recorro cada cuadrícula en los datos y encuentro si la distancia es inferior a 500 km. Luego almaceno una variable con VERDADERO o FALSO, que luego puedo usar para promediar los datos (otra variable). A partir de este método, quiero una matriz con VERDADERO o FALSO para las ubicaciones dentro de los 500 km del lat y lon que se muestran. ¿Hay un método más eficiente para hacer esto?

Answer 1

Las funciones dist* del paquete de geosphere están vectorizadas, por lo que solo necesita preparar mejor sus entradas. Prueba esto:

#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

Como mostró la respuesta @ Floo0, hay muchos cálculos innecesarios. Podemos seguir otra estrategia: primero determinamos el rango de lon y lat que puede estar más cerca que el umbral y luego los usamos solo para calcular la distancia:

#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

De esta manera, calculas solo lg+ln+lg*ln ( lg y ln son la longitud de latgood y longood ), es decir, 531 distancias, opuestas a las 259200 con mi método anterior.

Answer 2

Solo use hutils::haversine_distance(lat, lon, mylat, mylon) < 500 directamente.

Mejora las respuestas existentes por su rapidez y robustez. En particular, no se basa en la naturaleza reticulada de los datos y funcionará con vectores largos de coordenadas. A continuación se muestran los tiempos de 100.000 puntos

# 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

Answer 3

Tiempos:

Comparando @ nicola''s y mi versión da:

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

Mi solución original: (la segunda versión de Nicola Nicola es mucho más limpia y rápida).

Puedes hacer lo siguiente (explicación abajo)

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
      }
    }
  }

Explicación:

Para el bucle aplico la siguiente lógica:

outer_loop_state se inicializa con 0. Si se encuentra una fila con al menos un punto de trama dentro del círculo, la outer_loop_state está establecida en 1. Una vez que no haya más puntos dentro del círculo para una fila dada, se romperá.

La llamada distm en la versión @nicola básicamente hace lo mismo sin este truco. Entonces calcula todas las filas.

Código de tiempos:

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}
)