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