shortest path Copenhagen vs London UK
2023-10-04
Objective: 1. Using dodgr package to find shortest path.
2. Visualisation with tmap
London
Loading library
pacman::p_load(
dodgr,
sf,
tidygraph,
igraph,
dplyr,
tibble,
ggplot2,
units,
tmap,
osmdata,
leaflet,
readxl,
rio,
here
)Data
london <- read_excel("data/msoa-data.xls")## New names:
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `2005` -> `2005...13`
## • `2006` -> `2006...14`
## • `2007` -> `2007...15`
## • `2008` -> `2008...16`
## • `2009` -> `2009...17`
## • `2010` -> `2010...18`
## • `2011` -> `2011...19`
## • `2012` -> `2012...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`
## • `` -> `...30`
## • `` -> `...31`
## • `` -> `...32`
## • `` -> `...33`
## • `` -> `...34`
## • `` -> `...35`
## • `` -> `...36`
## • `` -> `...37`
## • `` -> `...38`
## • `` -> `...39`
## • `` -> `...40`
## • `` -> `...41`
## • `` -> `...42`
## • `` -> `...43`
## • `Couple household with dependent children` -> `Couple household with
## dependent children...44`
## • `Couple household without dependent children` -> `Couple household without
## dependent children...45`
## • `Lone parent household` -> `Lone parent household...46`
## • `One person household` -> `One person household...47`
## • `Other household Types` -> `Other household Types...48`
## • `Couple household with dependent children` -> `Couple household with
## dependent children...49`
## • `Couple household without dependent children` -> `Couple household without
## dependent children...50`
## • `Lone parent household` -> `Lone parent household...51`
## • `One person household` -> `One person household...52`
## • `Other household Types` -> `Other household Types...53`
## • `` -> `...54`
## • `` -> `...55`
## • `` -> `...56`
## • `` -> `...57`
## • `` -> `...58`
## • `` -> `...59`
## • `` -> `...60`
## • `` -> `...61`
## • `` -> `...62`
## • `` -> `...63`
## • `` -> `...64`
## • `` -> `...65`
## • `` -> `...66`
## • `` -> `...67`
## • `` -> `...68`
## • `` -> `...69`
## • `` -> `...70`
## • `` -> `...71`
## • `` -> `...72`
## • `` -> `...73`
## • `` -> `...74`
## • `` -> `...75`
## • `` -> `...76`
## • `` -> `...77`
## • `` -> `...78`
## • `` -> `...79`
## • `` -> `...80`
## • `` -> `...81`
## • `` -> `...82`
## • `` -> `...83`
## • `` -> `...84`
## • `` -> `...85`
## • `` -> `...86`
## • `` -> `...87`
## • `` -> `...88`
## • `` -> `...89`
## • `` -> `...90`
## • `` -> `...91`
## • `` -> `...92`
## • `` -> `...93`
## • `` -> `...94`
## • `` -> `...95`
## • `` -> `...96`
## • `` -> `...97`
## • `` -> `...98`
## • `` -> `...99`
## • `` -> `...100`
## • `` -> `...101`
## • `` -> `...102`
## • `` -> `...103`
## • `` -> `...104`
## • `` -> `...105`
## • `` -> `...106`
## • `` -> `...107`
## • `` -> `...108`
## • `` -> `...109`
## • `` -> `...110`
## • `` -> `...111`
## • `` -> `...112`
## • `` -> `...113`
## • `2005` -> `2005...114`
## • `2006` -> `2006...115`
## • `2007` -> `2007...116`
## • `2008` -> `2008...117`
## • `2009` -> `2009...118`
## • `2010` -> `2010...119`
## • `2011` -> `2011...120`
## • `2012` -> `2012...121`
## • `2005` -> `2005...123`
## • `2006` -> `2006...124`
## • `2007` -> `2007...125`
## • `2008` -> `2008...126`
## • `2009` -> `2009...127`
## • `2010` -> `2010...128`
## • `2011` -> `2011...129`
## • `2011` -> `2011...130`
## • `` -> `...132`
## • `` -> `...133`
## • `` -> `...134`
## • `` -> `...135`
## • `` -> `...136`
## • `` -> `...137`
## • `` -> `...138`
## • `` -> `...139`
## • `` -> `...140`
## • `` -> `...141`
## • `` -> `...142`
## • `` -> `...143`
## • `` -> `...144`
## • `` -> `...145`
## • `` -> `...146`
## • `` -> `...147`
## • `` -> `...148`
## • `` -> `...149`
## • `` -> `...150`
## • `` -> `...151`
## • `` -> `...152`
## • `` -> `...153`
## • `` -> `...154`
## • `` -> `...155`
## • `` -> `...156`
## • `` -> `...157`
## • `` -> `...158`
## • `` -> `...159`
## • `` -> `...160`
## • `` -> `...161`
## • `` -> `...162`
## • `` -> `...163`
## • `` -> `...164`
## • `` -> `...165`
## • `` -> `...166`
## • `` -> `...167`
## • `` -> `...168`
## • `` -> `...169`
## • `` -> `...170`
## • `` -> `...171`
## • `` -> `...172`
## • `` -> `...173`
## • `` -> `...174`
## • `` -> `...175`
## • `` -> `...176`
## • `` -> `...177`
## • `` -> `...178`
## • `` -> `...179`
## • `` -> `...180`
## • `` -> `...181`
## • `` -> `...182`
## • `` -> `...183`
## • `` -> `...184`
## • `` -> `...185`
## • `` -> `...186`
## • `` -> `...187`
## • `` -> `...188`
## • `` -> `...189`
## • `` -> `...190`
## • `` -> `...191`
## • `` -> `...192`
## • `` -> `...193`
## • `` -> `...194`
## • `` -> `...195`
## • `Fatal` -> `Fatal...196`
## • `Serious` -> `Serious...197`
## • `Slight` -> `Slight...198`
## • `Fatal` -> `Fatal...200`
## • `Serious` -> `Serious...201`
## • `Slight` -> `Slight...202`
## • `Fatal` -> `Fatal...204`
## • `Serious` -> `Serious...205`
## • `Slight` -> `Slight...206`england <- read.csv("data/MSOA_Dec_2011_PWC_in_England_and_Wales_2022_-7657754233007660732.csv")
england.2 <- read_excel("data/MLSOA_Dec_2021_PWC_in_England_and_Wales_2022_-3559472851201324412.xlsx")
colnames(england)[2] <- "MSOA Code"
colnames(england.2)[2] <- "MSOA Code"
england.3 <- left_join(england, england.2, by = "MSOA Code")
england.3 <- england.3 %>% select(c("MSOA Code", "MSOA11NM", "x.y", "y.y"))
london <- london %>% select("MSOA Code")
london.city <- left_join(london, england.3, by = "MSOA Code")london.city <- left_join(london, england, by = "MSOA Code")rac <- read.csv("data/OD_Table.csv")colnames(rac)[1] <- "MSOA Code"
tes <- left_join(rac, london.city, by = "MSOA Code")
test <- tes[1:30, ]Transform to {sf} data structure
# transform to {sf} data structure
dt <- st_as_sf(test, coords = c("longitude", "latitude"), crs = 4326)
class(dt)## [1] "sf" "data.frame"Map
# leaflet map for points
leaflet() %>% addTiles() %>% addMarkers(data = dt)# route map with osrm
osroute <- osrm::osrmRoute(loc = dt)
leaflet(data = dt) %>%
addProviderTiles("CartoDB.Positron") %>%
addPolylines(data = osroute,
label = "OSRM engine",
color = "red")Shortest paths using dodgr package
data
net <- dodgr_streetnet("city of london uk")
class(net)## [1] "sf" "data.frame"net.1 <- net[, c("osm_id", "geometry")]
head(net.1, 3)## Simple feature collection with 3 features and 1 field
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -0.0884115 ymin: 51.50617 xmax: -0.0821093 ymax: 51.51287
## Geodetic CRS: WGS 84
## osm_id geometry
## 554369 554369 LINESTRING (-0.0884017 51.5...
## 554526 554526 LINESTRING (-0.0851265 51.5...
## 742018 742018 LINESTRING (-0.0882957 51.5...map
# tmap
tm_shape(net.1) +
tm_layout(legend.outside=TRUE,bg='black') +
tm_lines("osm_id",lwd=1.5)## Warning: Number of levels of the variable "osm_id" is 9960, which is larger
## than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 9960) in the layer function to show all levels.
library(osmplotr)## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright##
## Attaching package: 'osmplotr'## The following object is masked _by_ '.GlobalEnv':
##
## londonlibrary(magrittr)
# osm base map
map <- osm_basemap(net, bg = "gray95") %>%
add_osm_objects(net, col = "gray5") %>%
add_axes() %>%
print_osm_map()shortest path
graph <- weight_streetnet(net, wt_profile = "motorcar") # convert from sf to dataframe, and weightweight it for pedestrian travel ## The following highway types are present in data yet lack corresponding weight_profile values: construction, busway, NA, elevator, road, corridor, proposed,from <- sample(graph$from_id, size = 20)
to <- sample(graph$to_id, size = 5)
# shortest path
dp <- dodgr_paths(graph, from = from, to = to)
dp [[1]] [[1]]## [1] "819157946" "9782353733" "7840826007" "1771462115" "5938625692"
## [6] "110146" "4401149224" "5927671120" "4401149225" "5927671155"
## [11] "4401149228" "4992781187" "11018999608" "11018999606" "7636278195"
## [16] "11018999607" "25477148" "11018908598" "11018908601" "11018999611"
## [21] "11018908600" "4826555771" "9030155762" "9042037853" "914822492"
## [26] "9030155760" "4992781180" "1312335310" "2387245026" "2387245027"
## [31] "2387245029" "9030155757" "9042120207" "9030155755" "9030155756"
## [36] "250327618" "9808795603" "21990630" "5565443390" "5565443397"
## [41] "5875062471" "9808795604" "5565443389" "21990629" "8923643489"
## [46] "9042053623" "5875062472" "8923643490" "21990628" "9030155754"
## [51] "9042053624" "9030155749" "6985373982" "1350430181" "1721713699"
## [56] "10677374" "1312264610" "1782718791" "6985373983" "11018990300"
## [61] "11018990299" "2902624645" "4835108996" "1144619335" "4835108997"
## [66] "7223326412" "1939558704" "2218362321" "1939558701" "1939558697"
## [71] "10660433979" "1939558667" "1939558657" "1939558651" "2218362310"
## [76] "1939558650" "2218362314" "1939558649" "1939558648" "1939558647"
## [81] "5925584204" "2961557" "24950101" "5925584203" "24950102"
## [86] "24950105" "4835108995" "2961540" "24950109" "2961549"
## [91] "24950113" "4077593104" "24950114" "8807431535" "9027328313"
## [96] "4867087535" "9027328314" "5065949596" "5049840242" "6985395258"
## [101] "5065949598" "5049840243" "1951616283" "1782714747" "1951616278"
## [106] "1951616276" "5473273376" "1951616261" "1951616231" "107627"
## [111] "4072083344" "4072083345" "259252667" "1312454045" "5925584201"
## [116] "107626" "1312454055" "8390582613" "25532464" "1547013719"
## [121] "1547013589" "8785024515" "1547013554" "1547013710" "25532494"
## [126] "10544113918" "10544113927" "1547013705" "10544113926" "4072083342"
## [131] "4072083343" "3385387607" "6985465129" "3385387604" "3385387605"
## [136] "6985465125" "3385387606" "1556324759" "24952195" "6985465124"
## [141] "6985465130" "10236420606" "6985465145" "4117882852" "4117882856"
## [146] "4117882857" "4117882860" "25278669" "25278670" "696477399"
## [151] "20625790" "4098680799" "8807410124" "1556324726" "4098680805"
## [156] "4098680798" "4098680807" "4098680806" "6967608821" "6508724312"
## [161] "4098680797" "4098680803" "4098680809" "4098680812" "1951372953"verts <- dodgr_vertices (graph)
path1 <- verts [match (dp [[1]] [[1]], verts$id), ]
head (path1)## id x y component n
## 5298 819157946 -0.0870064 51.51903 1 3511
## 16164 9782353733 -0.0872906 51.51912 1 9753
## 2420 7840826007 -0.0874469 51.51917 1 1641
## 2418 1771462115 -0.0874921 51.51919 1 1640
## 2416 5938625692 -0.0875338 51.51921 1 1639
## 1311 110146 -0.0876681 51.51926 1 887# map
dt <- st_as_sf(path1, coords = c("x", "y"), crs = 4326)
dt <- dt[1:100, ]
leaflet() %>% addTiles() %>% addMarkers(data = dt)Copenhagen
Geting bounding box
bb <- osmdata::getbb("Copenhagen")
bb## min max
## x 12.41007 12.73007
## y 55.52672 55.84672# drawing random points in a square
npts <- 10
xy <- apply(bb, 1, function(i) min(i) + runif(npts)*diff(i))head(xy, 10)## x y
## [1,] 12.61972 55.78185
## [2,] 12.48712 55.54049
## [3,] 12.46056 55.80273
## [4,] 12.63732 55.60916
## [5,] 12.70213 55.71405
## [6,] 12.45863 55.69064
## [7,] 12.56064 55.75137
## [8,] 12.61898 55.78450
## [9,] 12.61476 55.69609
## [10,] 12.61365 55.70648Shortest path with dodgr
# download street network
net <- dodgr_streetnet(bb)
trans_net <- weight_streetnet(net, wt_profile = "foot")## The following highway types are present in data yet lack corresponding weight_profile values: proposed, construction, corridor, busway, NA, platform, raceway, disused, rest_area, elevator, bus_stop,head(trans_net, 3)## geom_num edge_id from_id from_lon from_lat to_id to_lon to_lat
## 1 1 1 25981317 12.54873 55.71641 10845990327 12.54879 55.71638
## 2 1 2 10845990327 12.54879 55.71638 25981317 12.54873 55.71641
## 3 1 3 10845990327 12.54879 55.71638 422248375 12.54904 55.71623
## d d_weighted highway way_id component time time_weighted
## 1 5.612738 7.483651 residential 2796 1 4.041171 5.388228
## 2 5.612738 7.483651 residential 2796 1 4.041171 5.388228
## 3 22.043362 29.391149 residential 2796 1 15.871221 21.161627d <- dodgr_dists(trans_net, from = xy, to = xy)# Creat a small net
bbox <- matrix(c(12.5926,55.6612,12.6325,55.6728),ncol=2)
colnames(bbox) <- c("min","max")
rownames(bbox) <- c("x","y")
small_net <- dodgr_streetnet(bbox)temp <- small_net[, c("osm_id", "geometry")]tm_shape(temp) +
tm_layout(legend.outside = TRUE, bg = "black") +
tm_lines("osm_id", lwd = 1.5)## Warning: Number of levels of the variable "osm_id" is 1775, which is larger
## than max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 1775) in the layer function to show all levels.