#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
пакетом, если у вас еще нет собственных данных.