left_join на основе ближайшего LAT_LON в R

Question

Я пытаюсь найти идентификатор ближайшего LAT_LON в data.frame со ссылкой на мой оригинальный data.frame. Я уже понял это, объединив оба data.frames с уникальным идентификатором и вычислив расстояние на основе функции distHaverSine из geosphere. Теперь я хочу сделать еще один шаг и присоединиться к data.frames без уникального идентификатора и найти идентификатор ближайшего LAT-LON. Я использовал следующий код после слияния:

v3 <-v2 %>% mutate(CTD = distHaversine(cbind(LON.x, LAT.x), cbind(LON.y, LAT.y)))

ДАННЫЕ:

loc <- data.frame(station = c('Baker Street','Bank'),
     lat = c(51.522236,51.5134047),
     lng = c(-0.157080, -0.08905843),
               postcode = c('NW1','EC3V'))

stop <- data.frame(station = c('Angel','Barbican','Barons Court','Bayswater'),
                lat = c(51.53253,51.520865,51.490281,51.51224),
                lng = c(-0.10579,-0.097758,-0.214340,-0.187569),
                postcode = c('EC1V','EC1A', 'W14', 'W2'))

В качестве конечного результата я хотел бы что-то вроде этого:

df <- data.frame(loc = c('Baker Street','Bank','Baker Street','Bank','Baker Street','Bank','Baker 
        Street','Bank'), 
              stop = c('Angel','Barbican','Barons Court','Bayswater','Angel','Barbican','Barons Court','Bayswater'), 
              dist = c('x','x','x','x','x','x','x','x'), 
              lat = c(51.53253,51.520865,51.490281,51.51224,51.53253,51.520865,51.490281,51.51224), 
              lng = c(-0.10579,-0.097758,-0.214340,-0.187569,-0.10579,-0.097758,-0.214340,-0.187569),
              postcode = c('EC1V','EC1A', 'W14', 'W2','EC1V','EC1A', 'W14', 'W2')
              )

Любая помощь приветствуется. Спасибо.

Answer 1

Поскольку расстояния между объектами малы, мы можем ускорить вычисления, используя евклидово расстояние между координатами. Поскольку мы не вокруг экватора, lng координаты немного сжимаются; мы можем сделать сравнение немного лучше, немного расширив lng.

cor_stop <- stop[, c("lat", "lng")]
cor_stop$lng <- cor_stop$lng * sin(mean(cor_stop$lat, na.rm = TRUE)/180*pi)
cor_loc <- loc[, c("lat", "lng")]
cor_loc$lng <- cor_loc$lng * sin(mean(cor_loc$lat, na.rm = TRUE)/180*pi)

Затем мы можем вычислить ближайшую остановку для каждого местоположения, используя пакет FNN, который использует поиск по дереву для быстрого поиска ближайших соседей K. Это должно масштабироваться до больших наборов данных (я использовал это для наборов данных с миллионами записей):

library(FNN)
matches <- knnx.index(cor_stop, cor_loc, k = 1)
matches

##      [,1]
## [1,]    4
## [2,]    2

Затем мы можем построить конечный результат:

res <- loc
res$stop_station  <- stop$station[matches[,1]]
res$stop_lat      <- stop$lat[matches[,1]]
res$stop_lng      <- stop$lng[matches[,1]]
res$stop_postcode <- stop$postcode[matches[,1]]

И вычислите фактическое расстояние:

library(geosphere)
res$dist <- distHaversine(res[, c("lng", "lat")], res[, c("stop_lng", "stop_lat")])
res

##          station      lat         lng postcode stop_station stop_lat  stop_lng
## 1 Baker Street 51.52224 -0.15708000      NW1    Bayswater 51.51224 -0.187569
## 2         Bank 51.51340 -0.08905843     EC3V     Barbican 51.52087 -0.097758
##   stop_postcode     dist
## 1            W2 2387.231
## 2          EC1A 1026.091

Если вы не уверены, что самая близкая точка в широте-долготе также является самой близкой точкой «как птица летит», вы можете использовать этот метод, чтобы сначала выбрать K ближайших точек в широте; затем рассчитайте расстояния для этих точек и затем выберите ближайшую точку.

Answer 2

Все соединения, расчеты расстояний и черчения могут быть выполнены с использованием доступных пакетов R.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3
library(nngeo)
library(mapview)

## Original data
loc <- data.frame(station = c('Baker Street','Bank'),
                  lat = c(51.522236,51.5134047),
                  lng = c(-0.157080, -0.08905843),
                  postcode = c('NW1','EC3V'))

stop <- data.frame(station = c('Angel','Barbican','Barons Court','Bayswater'),
                   lat = c(51.53253,51.520865,51.490281,51.51224),
                   lng = c(-0.10579,-0.097758,-0.214340,-0.187569),
                   postcode = c('EC1V','EC1A', 'W14', 'W2'))

df <- data.frame(loc = c('Baker Street','Bank','Baker Street','Bank','Baker Street','Bank','Baker 
        Street','Bank'), 
                 stop = c('Angel','Barbican','Barons Court','Bayswater','Angel','Barbican','Barons Court','Bayswater'), 
                 dist = c('x','x','x','x','x','x','x','x'), 
                 lat = c(51.53253,51.520865,51.490281,51.51224,51.53253,51.520865,51.490281,51.51224), 
                 lng = c(-0.10579,-0.097758,-0.214340,-0.187569,-0.10579,-0.097758,-0.214340,-0.187569),
                 postcode = c('EC1V','EC1A', 'W14', 'W2','EC1V','EC1A', 'W14', 'W2')
)



## Create sf objects from lat/lon points
loc_sf <- loc %>% st_as_sf(coords = c('lng', 'lat'), remove = T) %>%
  st_set_crs(4326) 

stop_sf <- stop %>% st_as_sf(coords = c('lng', 'lat'), remove = T) %>%
  st_set_crs(4326) 


# Use st_nearest_feature to cbind loc to stop by nearest points
joined_sf <- stop_sf %>% 
  cbind(
    loc_sf[st_nearest_feature(stop_sf, loc_sf),])


## mutate to add column showing distance between geometries
joined_sf %>%
  mutate(dist = st_distance(geometry, geometry.1, by_element = T))
#> Simple feature collection with 4 features and 5 fields
#> Active geometry column: geometry
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: -0.21434 ymin: 51.49028 xmax: -0.097758 ymax: 51.53253
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#>        station postcode    station.1 postcode.1                   geometry
#> 1        Angel     EC1V         Bank       EC3V  POINT (-0.10579 51.53253)
#> 2     Barbican     EC1A         Bank       EC3V POINT (-0.097758 51.52087)
#> 3 Barons Court      W14 Baker Street        NW1  POINT (-0.21434 51.49028)
#> 4    Bayswater       W2 Baker Street        NW1 POINT (-0.187569 51.51224)
#>                    geometry.1         dist
#> 1 POINT (-0.08905843 51.5134) 2424.102 [m]
#> 2 POINT (-0.08905843 51.5134) 1026.449 [m]
#> 3   POINT (-0.15708 51.52224) 5333.417 [m]
#> 4   POINT (-0.15708 51.52224) 2390.791 [m]



## Use nngeo and mapview to plot lines on a map
# NOT run for reprex, output image attached 
#connected <- st_connect(stop_sf, loc_sf)
# mapview(connected) + 
#   mapview(loc_sf, color = 'red') +
#   mapview(stop_sf, color = 'black')

^{Создано в 2020-01-21 с помощью пакета Представить (v0.3.0)}

Answer 3

Вы можете полностью избежать поиска ближайших соседей, если способны использовать проекционную систему координат. Если вы можете, то вы можете дешево построить полигоны Вороного вокруг каждого местоположения - эти полигоны определяют области, которые являются ближайшими к каждой из входных точек.

Затем вы можете просто использовать пересечения ГИС, чтобы найти точки l ie в каких полигонах, а затем вычислить расстояния для известных пар ближайших точек. Я думаю, что это должно быть намного быстрее. Однако вы не можете использовать полигоны Вороного с географическими c координатами.

loc <- data.frame(station = c('Baker Street','Bank'),
     lat = c(51.522236,51.5134047),
     lng = c(-0.157080, -0.08905843),
               postcode = c('NW1','EC3V'))

stop <- data.frame(station = c('Angel','Barbican','Barons Court','Bayswater'),
                lat = c(51.53253,51.520865,51.490281,51.51224),
                lng = c(-0.10579,-0.097758,-0.214340,-0.187569),
                postcode = c('EC1V','EC1A', 'W14', 'W2'))

# Convert to a suitable PCS (in this case OSGB)
stop <- st_as_sf(stop, coords=c('lng','lat'), crs=4326)
stop <- st_transform(stop, crs=27700)
loc <- st_as_sf(loc, coords=c('lng','lat'), crs=4326)
loc <- st_transform(loc, crs=27700)

# Extract Voronoi polygons around locations and convert to an sf object
loc_voronoi <- st_collection_extract(st_voronoi(do.call(c, st_geometry(loc))))
loc_voronoi <- st_sf(loc_voronoi, crs=crs(loc))

# Match Voronoi polygons to locations and select that geometry
loc$voronoi <- loc_voronoi$loc_voronoi[unlist(st_intersects(loc, loc_voronoi))]
st_geometry(loc) <- 'voronoi'

# Find which stop is closest to each location
stop$loc <- loc$station[unlist(st_intersects(stop, loc))]

# Reset locs to use the point geometry and get distances
st_geometry(loc) <- 'geometry'
stop$loc_dist <- st_distance(stop, loc[stop$loc,], by_element=TRUE)

Это дает следующий вывод:

Simple feature collection with 4 features and 4 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 524069.7 ymin: 178326.3 xmax: 532074.6 ymax: 183213.9
epsg (SRID):    27700
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
       station postcode                  geometry          loc     loc_dist
1        Angel     EC1V POINT (531483.8 183213.9)         Bank 2423.722 [m]
2     Barbican     EC1A POINT (532074.6 181931.2)         Bank 1026.289 [m]
3 Barons Court      W14 POINT (524069.7 178326.3) Baker Street 5332.478 [m]
4    Bayswater       W2 POINT (525867.7 180813.9) Baker Street 2390.377 [m]