Есть ли способ геокодировать пересечение на основе двух названий улиц в R?

#r #coordinates #sf

Вопрос:

У меня есть набор данных, который выглядит так

 ID|Cross Street 1|Cross Street 2
 1 Beaver St      Hanover St
 2 Pearl St       Wall St
 

и я хочу, чтобы в итоге получилось 2263 значения x и y

 ID|Cross Street 1|Cross Street 2 x              |y
 1 Beaver St      Hanover St     981740.53187633 196247.34676349
 2 Pearl St       Wall St        982049.05259918 196320.89988222
 

Обратите внимание, что широта и долгота указаны здесь ниже, чтобы подтвердить местоположение.

 lat              |lon
40.705330         -74.009051
40.7055319680708  -74.00793826989502
 

Комментарии:

1. У вас есть шейп-файл дорог в этом районе, или вы ищете подход к общедоступному api (или в стиле карт Google)?

2. У меня нет шейп-файла данных. Я хочу работать с csv-файлом.

3. Координаты lat/lon должны быть откуда-то получены. Просто двух названий пересекающихся улиц без пространственных данных и информации о городе/штате/провинции/и т. Д. Будет недостаточно.

4. почтового индекса или района недостаточно, чтобы получить какую-то справку о местоположении?

Ответ №1:

Вот пример использования данных, загруженных с помощью osmdata пакета (OpenStreetMap) для небольшой части города Канберра, Австралия.

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

 library(sf)
library(tidyverse)
library(osmdata)
library(tmap)
tmap_mode('view')

# sample data with our focal area
tribble(
  ~point, ~lat, ~lon, 
  1, -35.29, 149.12, 
  2, -35.29, 149.14, 
  3, -35.27, 149.14, 
  4, -35.27, 149.12, 
) %>% 
  st_as_sf(
    coords = c('lon', 'lat'), 
    crs = 4326
  ) %>% 
  {. ->> my_points}

my_points

# Simple feature collection with 4 features and 1 field
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 149.12 ymin: -35.29 xmax: 149.14 ymax: -35.27
# Geodetic CRS:  WGS 84
# # A tibble: 4 x 2
# point        geometry
# * <dbl>     <POINT [°]>
# 1     1 (149.12 -35.29)
# 2     2 (149.14 -35.29)
# 3     3 (149.14 -35.27)
# 4     4 (149.12 -35.27)
 

И нарисуйте его, чтобы взглянуть на местность:

 tm_shape(my_points) 
  tm_dots()
 

введите описание изображения здесь

Затем мы импортируем данные OSM. Я здесь не эксперт, я просто следую примеру отсюда.

 # import OSM data
st_bbox(my_points) %>% 
  opq %>% 
  add_osm_feature(key = 'highway') %>% 
  osmdata_sf %>% 
  `[[`('osm_lines') %>% 
  {. ->> my_streets}

my_streets

# Simple feature collection with 2143 features and 87 fields
# Geometry type: LINESTRING
# Dimension:     XY
# Bounding box:  xmin: 149.1153 ymin: -35.29224 xmax: 149.1447 ymax: -35.26598
# Geodetic CRS:  WGS 84
# First 10 features:
#          osm_id              name access amenity area bicycle bridge bridge_number building  bus bus.lanes busway
# 4018757 4018757 Coranderrk Street   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4189987 4189987    London Circuit   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4189988 4189988    London Circuit   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4204003 4204003        Parkes Way   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4204005 4204005        Parkes Way   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4204006 4204006        Parkes Way   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4368141 4368141        Parkes Way   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4414782 4414782     Watson Street   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4414873 4414873      Gould Street   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
# 4574783 4574783    Macleay Street   <NA>    <NA> <NA>    <NA>   <NA>          <NA>     <NA> <NA>      <NA>   <NA>
 

Для примера я сохраняю только residential улицы ( highway столбец), а затем использую group_by() и summarise() для объединения всех сегментов улицы с одинаковым названием, прежде чем удалять улицы без названия.

 # keep only residential streets
my_streets %>% 
  filter(
    highway == 'residential'
  ) %>% 
  select(name) %>% 
  group_by(name) %>% 
  summarise %>%
  filter(
    !is.na(name)
  ) %>% 
  {. ->> my_streets_2}

my_streets_2

# Simple feature collection with 31 features and 1 field
# Geometry type: GEOMETRY
# Dimension:     XY
# Bounding box:  xmin: 149.1201 ymin: -35.28725 xmax: 149.1437 ymax: -35.26933
# Geodetic CRS:  WGS 84
# # A tibble: 31 x 2
# name                                                                                                          geometry
#   * <chr>                                                                                                   <GEOMETRY [°]>
#   1 Ainslie Avenue     MULTILINESTRING ((149.1355 -35.27924, 149.1356 -35.2792, 149.1357 -35.27915, 149.1358 -35.27909, .~
#   2 Akuna Street       LINESTRING (149.1361 -35.28036, 149.136 -35.2804, 149.1356 -35.28057, 149.1353 -35.28068, 149.134.~
#   3 Allambee Street    LINESTRING (149.1378 -35.27985, 149.1379 -35.27967, 149.138 -35.27962, 149.1391 -35.27918, 149.13.~
#   4 Allara Street      LINESTRING (149.131 -35.28574, 149.131 -35.28568, 149.1311 -35.28563, 149.1312 -35.2855, 149.1316.~
#   5 Amaroo Street      MULTILINESTRING ((149.142 -35.28725, 149.1419 -35.28721, 149.1417 -35.2871, 149.1415 -35.28698, 1.~
#   6 Batman Street      MULTILINESTRING ((149.1366 -35.27775, 149.1367 -35.2777, 149.138 -35.27721, 149.1381 -35.27717, 1.~
#   7 Binara Street      LINESTRING (149.1338 -35.28256, 149.1339 -35.2825, 149.134 -35.28245, 149.1341 -35.28239, 149.134.~
#   8 Boolee Street      MULTILINESTRING ((149.1368 -35.2813, 149.1369 -35.28125, 149.1376 -35.281, 149.1381 -35.28079, 14.~
#   9 Booroondara Street LINESTRING (149.1427 -35.28638, 149.1427 -35.28634, 149.1402 -35.28487, 149.14 -35.28479, 149.139.~
#   10 Bunda Street       LINESTRING (149.1348 -35.28205, 149.1349 -35.28193, 149.1348 -35.28181, 149.1347 -35.2815, 149.13.~
#   # ... with 21 more rows
 

Опять же, составьте график, чтобы увидеть, с чем мы работаем:

 tm_shape(my_streets_2) 
  tm_lines(lwd = 2, col = 'red')
 

введите описание изображения здесь

Теперь мы используем map() итерацию по каждой улице в коллекции и находим, где она пересекается с другими улицами (кроме самой себя), используя st_intersection() .

 # now, find which ones intersect
1:nrow(my_streets_2) %>% 
  map(function(i){
    
    my_streets_2 %>% 
      filter(
        row_number() == i
      ) %>% 
      st_intersection(
        my_streets_2 %>% 
          filter(
            row_number() != i
          )
      )
      
    }) %>% 
  bind_rows %>% 
  {. ->> my_intersections}

my_intersections

# Simple feature collection with 50 features and 2 fields
# Geometry type: GEOMETRY
# Dimension:     XY
# Bounding box:  xmin: 149.1261 ymin: -35.28565 xmax: 149.1411 ymax: -35.27228
# Geodetic CRS:  WGS 84
# # A tibble: 50 x 3
# name            name.1                                                                                        geometry
# * <chr>           <chr>                                                                                   <GEOMETRY [°]>
# 1 Ainslie Avenue  Currong Street North                            MULTIPOINT ((149.1371 -35.27858), (149.1372 -35.2789))
# 2 Ainslie Avenue  Doonkuna Street      MULTIPOINT ((149.1386 -35.27798), (149.1385 -35.27801), (149.1387 -35.27832), (1~
# 3 Ainslie Avenue  Elimatta Street      MULTIPOINT ((149.1403 -35.2777), (149.1402 -35.27773), (149.14 -35.27742), (149.~
# 4 Ainslie Avenue  Kogarah Lane                                                                POINT (149.1365 -35.27917)
# 5 Allambee Street Doonkuna Street                                                             POINT (149.1392 -35.27914)
# 6 Amaroo Street   Euree Street                                                                POINT (149.1393 -35.28565)
# 7 Batman Street   Currong Street North                                                        POINT (149.1366 -35.27775)
# 8 Batman Street   Doonkuna Street                                                             POINT (149.1381 -35.27717)
# 9 Batman Street   Elimatta Street                                                             POINT (149.1396 -35.27657)
# 10 Batman Street   Gooreen Street                                                              POINT (149.1411 -35.27599)
# # ... with 40 more rows
 

Теперь некоторые улицы пересекаются с другими улицами несколько раз (например, пересечение двух двухполосных дорог), и они представлены в виде MULTIPOINT геометрий, в то время как другие являются одиночными POINT геометриями.

Чтобы извлечь координаты (позже), я разделил все это на POINT геометрии, используя st_cast() . Обратите внимание, что предупреждающее сообщение указывает, что при этом сохраняется только первая геометрия каждой группы. Ради этого упражнения я продолжу, но если вам нужно каждое пересечение, вам, возможно, потребуется изменить эту команду.

Я также сохраняю только уникальные комбинации улиц, так как в настоящее время появляются повторяющиеся перекрестки с переменными name и name.1 значениями. Это вдвое сокращает количество функций. В качестве альтернативы вы могли бы использовать distinct(geometry) , однако это привело к уменьшению количества функций на 4 (я предполагаю , что это связано с альтернативными названиями одной и той же улицы или, возможно, перекрестков с 3 или 4 путями с одинаковыми координатами).

 my_intersections %>% 
  st_cast('POINT') %>%
  rowwise %>%
  mutate(
    streets = list(sort(c(name, name.1)))
  ) %>%
  distinct(streets, .keep_all = TRUE) %>%
  select(-streets) %>% 
  {. ->> my_int_2}

my_int_2

# Simple feature collection with 25 features and 2 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 149.1261 ymin: -35.28565 xmax: 149.1411 ymax: -35.27228
# Geodetic CRS:  WGS 84
# # A tibble: 25 x 3
# # Rowwise: 
#   name            name.1                           geometry
#   <chr>           <chr>                         <POINT [°]>
# 1 Ainslie Avenue  Currong Street North (149.1371 -35.27858)
# 2 Ainslie Avenue  Doonkuna Street      (149.1386 -35.27798)
# 3 Ainslie Avenue  Elimatta Street       (149.1403 -35.2777)
# 4 Ainslie Avenue  Kogarah Lane         (149.1365 -35.27917)
# 5 Allambee Street Doonkuna Street      (149.1392 -35.27914)
# 6 Amaroo Street   Euree Street         (149.1393 -35.28565)
# 7 Batman Street   Currong Street North (149.1366 -35.27775)
# 8 Batman Street   Doonkuna Street      (149.1381 -35.27717)
# 9 Batman Street   Elimatta Street      (149.1396 -35.27657)
# 10 Batman Street   Gooreen Street       (149.1411 -35.27599)
# # ... with 15 more rows
 

Теперь, как вы просили, мы заменяем столбец geometry столбцами X и Y (хотя, если вы намерены выполнять какие-либо дополнительные пространственные операции с этими данными, было бы разумно сохранить geometry ). Это включает в себя извлечение координат с помощью st_coordinates() и добавление их в sf коллекцию с помощью bind_cols() .

 my_int_2 %>% 
  bind_cols(my_int_2 %>% st_coordinates %>% as_tibble) %>% 
  st_drop_geometry %>% 
  {. ->> my_int_3}

my_int_3

# # A tibble: 25 x 4
# name            name.1                   X     Y
# * <chr>           <chr>                <dbl> <dbl>
# 1 Ainslie Avenue  Currong Street North  149. -35.3
# 2 Ainslie Avenue  Doonkuna Street       149. -35.3
# 3 Ainslie Avenue  Elimatta Street       149. -35.3
# 4 Ainslie Avenue  Kogarah Lane          149. -35.3
# 5 Allambee Street Doonkuna Street       149. -35.3
# 6 Amaroo Street   Euree Street          149. -35.3
# 7 Batman Street   Currong Street North  149. -35.3
# 8 Batman Street   Doonkuna Street       149. -35.3
# 9 Batman Street   Elimatta Street       149. -35.3
# 10 Batman Street   Gooreen Street        149. -35.3
# # ... with 15 more rows
 

И мы, конечно, делаем сюжет, чтобы посмотреть, как мы прошли.

 tm_shape(my_streets_2) 
  tm_lines(lwd = 2, col = 'red') 
  tm_shape(my_int_2) 
  tm_dots()
 

введите описание изображения здесь

Похоже, там есть пара немаркированных перекрестков, это могут быть:

  • пересечения двух улиц в а MULTIPOINT , где POINT сохранилась только первая
  • улицы с тем же названием
  • геометрии, которые на самом деле не пересекаются (вы могли бы использовать st_buffer() небольшое расстояние, чтобы увеличить количество пересечений, если бы вам пришлось)

Надеюсь, это приведет вас в правильном направлении, ознакомьтесь с osmdata пакетом, если у вас еще нет собственных данных.