packages = c('igraph', 'tidygraph', 'ggraph', 'visNetwork', 'lubridate', 'tidyverse', 'ggplot2', 'ggmap', 'leaflet', 'ggiraph', 'sf', 'tmap')
for(p in packages){library
if(!require(p, character.only = T)){
install.packages(p)
}
library(p, character.only = T)
}
#top_bus <- c('51','61','961','66','67','14','30','143','10','147','178','197','154','167','166','174','196','171','2','12','21','170','970','13')
busstops <- read_csv('data/busstops_with_planning_area_XY.csv')
busstops$BusStopCode <- as.character(busstops$BusStopCode)
busstops <- busstops[c('BusStopCode', 'Description', 'RoadName', 'planning_area', 'Latitude', 'Longitude')] %>%
filter(planning_area == 'QUEENSTOWN')
#busstops <- busstops[busstops$BusStopCode %in% as.list(unique(busroute['BusStopCode']))[['BusStopCode']], ]
busstops
## # A tibble: 226 x 6
## BusStopCode Description RoadName planning_area Latitude Longitude
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 10289 Opp SIS Bldg Alexandra ~ QUEENSTOWN 1.29 104.
## 2 10299 Aft Dawson Rd Alexandra ~ QUEENSTOWN 1.29 104.
## 3 10321 Aft Margaret Dr Tanglin Rd QUEENSTOWN 1.29 104.
## 4 10469 Bef Margaret Dr Kay Siang ~ QUEENSTOWN 1.30 104.
## 5 10479 Aft Tanglin Rd Kay Siang ~ QUEENSTOWN 1.30 104.
## 6 11009 Ghim Moh Ter Ghim Moh Rd QUEENSTOWN 1.31 104.
## 7 11021 Aft Queenstown NPC ~ Queensway QUEENSTOWN 1.29 104.
## 8 11029 Blk 166 Queensway QUEENSTOWN 1.29 104.
## 9 11031 Opp Queen's Cl Queensway QUEENSTOWN 1.29 104.
## 10 11039 Blk 19 Queensway QUEENSTOWN 1.29 104.
## # ... with 216 more rows
busroute <- read_csv('data/bus_route_overall.csv')
busroute$BusStopCode <- as.character(busroute$BusStopCode)
busroute <- busroute[c('BusStopCode', 'Direction', 'Distance', 'ServiceNo', 'StopSequence')]
# %>% dplyr::filter(ServiceNo %in% top_bus)
busroute <- busroute[busroute$BusStopCode %in% as.list(unique(busstops['BusStopCode']))[['BusStopCode']], ]
busroute
## # A tibble: 1,144 x 5
## BusStopCode Direction Distance ServiceNo StopSequence
## <chr> <dbl> <dbl> <chr> <dbl>
## 1 15141 1 26.6 10 59
## 2 15151 1 26.9 10 60
## 3 15161 1 27.1 10 61
## 4 15191 1 27.6 10 62
## 5 15201 1 27.9 10 63
## 6 15221 1 28.5 10 64
## 7 16201 1 28.8 10 65
## 8 16011 1 29.1 10 66
## 9 16021 1 29.6 10 67
## 10 16031 1 29.9 10 68
## # ... with 1,134 more rows
busroute_2 <- busroute
busroute_2['StopSequence'] = busroute_2['StopSequence']-1
busroute_2['BusStopCode_dest'] = busroute_2['BusStopCode']
busroute_2 <- busroute_2[c('BusStopCode_dest', 'Direction', 'ServiceNo', 'StopSequence')]
busstops_from_to <- dplyr::inner_join(busroute, busroute_2, by =c('StopSequence', 'ServiceNo', 'Direction'))
busstops_from_to
## # A tibble: 975 x 6
## BusStopCode Direction Distance ServiceNo StopSequence BusStopCode_dest
## <chr> <dbl> <dbl> <chr> <dbl> <chr>
## 1 15141 1 26.6 10 59 15151
## 2 15151 1 26.9 10 60 15161
## 3 15161 1 27.1 10 61 15191
## 4 15191 1 27.6 10 62 15201
## 5 15201 1 27.9 10 63 15221
## 6 15221 1 28.5 10 64 16201
## 7 16201 1 28.8 10 65 16011
## 8 16011 1 29.1 10 66 16021
## 9 16021 1 29.6 10 67 16031
## 10 16031 1 29.9 10 68 16041
## # ... with 965 more rows
#join the two tables together
busroute_busstop <- dplyr::inner_join(busstops_from_to, busstops, by ='BusStopCode')
keeps <- c('BusStopCode', 'BusStopCode_dest', 'Direction', 'Distance', 'ServiceNo', 'Latitude', 'Longitude', 'planning_area', 'Description')
busroute_busstop <- busroute_busstop[, keeps, drop = FALSE] %>%
rename(from = BusStopCode) %>%
rename(to = BusStopCode_dest)
head(busroute_busstop)
## # A tibble: 6 x 9
## from to Direction Distance ServiceNo Latitude Longitude planning_area
## <chr> <chr> <dbl> <dbl> <chr> <dbl> <dbl> <chr>
## 1 15141 15151 1 26.6 10 1.27 104. QUEENSTOWN
## 2 15151 15161 1 26.9 10 1.27 104. QUEENSTOWN
## 3 15161 15191 1 27.1 10 1.27 104. QUEENSTOWN
## 4 15191 15201 1 27.6 10 1.28 104. QUEENSTOWN
## 5 15201 15221 1 27.9 10 1.28 104. QUEENSTOWN
## 6 15221 16201 1 28.5 10 1.28 104. QUEENSTOWN
## # ... with 1 more variable: Description <chr>
busroute_busstop_aggregated <- busroute_busstop %>%
group_by(from, to, planning_area) %>%
summarise(Weight = n()) %>%
filter(from!=to) %>%
filter(Weight > 1) %>%
ungroup()
busroute_busstop_aggregated$from <- as.character(busroute_busstop_aggregated$from)
busroute_busstop_aggregated$to <- as.character(busroute_busstop_aggregated$to)
busroute_busstop_aggregated
## # A tibble: 195 x 4
## from to planning_area Weight
## <chr> <chr> <chr> <int>
## 1 10299 10289 QUEENSTOWN 5
## 2 10321 10479 QUEENSTOWN 2
## 3 10469 11499 QUEENSTOWN 2
## 4 10479 10469 QUEENSTOWN 2
## 5 11009 11359 QUEENSTOWN 6
## 6 11021 11031 QUEENSTOWN 12
## 7 11029 11519 QUEENSTOWN 4
## 8 11031 11041 QUEENSTOWN 11
## 9 11039 11029 QUEENSTOWN 12
## 10 11041 11051 QUEENSTOWN 11
## # ... with 185 more rows
#busstops_nodes_a <- unique(busroute_busstop_aggregated[c('target', 'planning_area')]) %>%
# rename(Id = target)
#busstops_nodes_b <- unique(busroute_busstop_aggregated[c('source', 'planning_area')]) %>%
# rename(Id = source)
#total <- unique(rbind(busstops_nodes_a, busstops_nodes_b))
total <- busstops %>%
rename(id = BusStopCode)
total$id <- as.character(total$id)
total
## # A tibble: 226 x 6
## id Description RoadName planning_area Latitude Longitude
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 10289 Opp SIS Bldg Alexandra Rd QUEENSTOWN 1.29 104.
## 2 10299 Aft Dawson Rd Alexandra Rd QUEENSTOWN 1.29 104.
## 3 10321 Aft Margaret Dr Tanglin Rd QUEENSTOWN 1.29 104.
## 4 10469 Bef Margaret Dr Kay Siang Rd QUEENSTOWN 1.30 104.
## 5 10479 Aft Tanglin Rd Kay Siang Rd QUEENSTOWN 1.30 104.
## 6 11009 Ghim Moh Ter Ghim Moh Rd QUEENSTOWN 1.31 104.
## 7 11021 Aft Queenstown NPC HQ Queensway QUEENSTOWN 1.29 104.
## 8 11029 Blk 166 Queensway QUEENSTOWN 1.29 104.
## 9 11031 Opp Queen's Cl Queensway QUEENSTOWN 1.29 104.
## 10 11039 Blk 19 Queensway QUEENSTOWN 1.29 104.
## # ... with 216 more rows
bus_graph <- tbl_graph(nodes = total, edges = busroute_busstop_aggregated, directed = TRUE)
# activate(edges) %>%
# arrange(desc(Weight))
#%>%
# filter(id %in% union(.E()$to, .E()$from))
bus_graph
## # A tbl_graph: 226 nodes and 195 edges
## #
## # A directed simple graph with 52 components
## #
## # Node Data: 226 x 6 (active)
## id Description RoadName planning_area Latitude Longitude
## <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 10289 Opp SIS Bldg Alexandra Rd QUEENSTOWN 1.29 104.
## 2 10299 Aft Dawson Rd Alexandra Rd QUEENSTOWN 1.29 104.
## 3 10321 Aft Margaret Dr Tanglin Rd QUEENSTOWN 1.29 104.
## 4 10469 Bef Margaret Dr Kay Siang Rd QUEENSTOWN 1.30 104.
## 5 10479 Aft Tanglin Rd Kay Siang Rd QUEENSTOWN 1.30 104.
## 6 11009 Ghim Moh Ter Ghim Moh Rd QUEENSTOWN 1.31 104.
## # ... with 220 more rows
## #
## # Edge Data: 195 x 4
## from to planning_area Weight
## <int> <int> <chr> <int>
## 1 1 2 QUEENSTOWN 5
## 2 3 4 QUEENSTOWN 2
## 3 5 6 QUEENSTOWN 2
## # ... with 192 more rows
ggraph(bus_graph) +
geom_edge_link() +
geom_node_point()
bus_graph=bus_graph%>%mutate(betweenness_centrality = centrality_betweenness(normalized = T)) %>%mutate(closeness_centrality = centrality_closeness(normalized = T))%>%mutate(centrality_degree=centrality_degree(mode='out',normalized = T))%>%mutate(centrality_eigen=centrality_eigen(weights=bus_graph$Weight,directed=T))
bus_graph
## # A tbl_graph: 226 nodes and 195 edges
## #
## # A directed simple graph with 52 components
## #
## # Node Data: 226 x 10 (active)
## id Description RoadName planning_area Latitude Longitude betweenness_cen~
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 10289 Opp SIS Bl~ Alexand~ QUEENSTOWN 1.29 104. 0.00228
## 2 10299 Aft Dawson~ Alexand~ QUEENSTOWN 1.29 104. 0
## 3 10321 Aft Margar~ Tanglin~ QUEENSTOWN 1.29 104. 0
## 4 10469 Bef Margar~ Kay Sia~ QUEENSTOWN 1.30 104. 0.0000992
## 5 10479 Aft Tangli~ Kay Sia~ QUEENSTOWN 1.30 104. 0.000159
## 6 11009 Ghim Moh T~ Ghim Mo~ QUEENSTOWN 1.31 104. 0.000179
## # ... with 220 more rows, and 3 more variables: closeness_centrality <dbl>,
## # centrality_degree <dbl>, centrality_eigen <dbl>
## #
## # Edge Data: 195 x 4
## from to planning_area Weight
## <int> <int> <chr> <int>
## 1 1 2 QUEENSTOWN 5
## 2 3 4 QUEENSTOWN 2
## 3 5 6 QUEENSTOWN 2
## # ... with 192 more rows
G=bus_graph
#Normally we need to create two vectors in order to plot both the network data and the coordinates on a map
#plot_vector takes the longitude and latitude to be plotted on a map.
#plot_vector1 takes the lon, lat and the network centrality values.
plot_vector<- as.data.frame(cbind(V(G)$Longitude,V(G)$Latitude))
plot_vector1<- as.data.frame(cbind(V(G)$Longitude,V(G)$Latitude,V(G)$betweenness_centrality,V(G)$closeness_centrality))
#we are taking the edgelist which is being used to get the origin and destination of the airport data.
#edgelist[,1]- takes the origin values and the edgelist[,2] takes the destination values.
edgelist <- get.edgelist(G)
edgelist[,1]<-as.numeric(match(edgelist[,1],V(G)))
edgelist[,2]<-as.numeric(match(edgelist[,2],V(G)))
#the edges now consists of the edge between the origin and the destination values
edges <- data.frame(plot_vector[edgelist[,1],], plot_vector[edgelist[,2],])
#naming the coloumns obtained to plot
colnames(edges) <- c("X1","Y1","X2","Y2")
#Here we are taking ggmap as our base layer. on top of ggmap we are plotting the origin and destination coordinates using the plot_vector as our 2nd layer, then we are plotting the coordinates of all the centrality from the plot_vector1 as our 3rd layer and then we are ploting the geom_segment i.e. the edges as our 4th layer.
#register_google(key = ) #create your own key
#Singapore
#map <- get_map(location = c(lon = 103.8279248, lat = 1.3438295), zoom=11, source = "google", maptype = 'roadmap')
#Queenstown
map <- get_map(location = c(lon = 103.7662684, lat = 1.2869421), zoom=14, source = "google", maptype = 'toner')
p <- ggmap(map)
z=p +
geom_segment(aes(x=X1, y=Y1, xend = X2, yend = Y2), color='black',data=edges,arrow = arrow(length = unit(0.2,"cm")))+
#geom_point(aes(V1, V2), data=plot_vector1)+
geom_point(aes(x=V1,y=V2,size=plot_vector1$V3,colour=plot_vector1$V4),plot_vector1)+
scale_colour_gradientn(colours=rainbow(5))+scale_size_continuous(range = c(1, 10))
z
#we use the node3 to get our origin and destination airport for plotting it on our map
node1=data.frame(plot_vector1[edgelist[,1],])
node2=data.frame(plot_vector1[edgelist[,2],])
node3=data.frame(cbind(node1,node2))
#Then we use a for loop in order to put the destination and origin airport on to the leaflet
map3 = leaflet(node3) %>%
addProviderTiles("CartoDB", group = "CartoDB") %>%
addTiles()
for(i in 1:nrow(node3)){
map3 <- addPolylines(map3, lat = as.numeric(node3[i, c(2, 6)]),
lng = as.numeric(node3[i, c(1, 5)]))
#origin
map3<-addCircleMarkers(map3, lat = as.numeric(node3[i, 2]),
lng = as.numeric(node3[i,1]),radius =(as.numeric(node3[i, 4]))*100,color='red')
#destination
map3<-addCircleMarkers(map3, lat = as.numeric(node3[i, 6]),
lng = as.numeric(node3[i, 5]),radius =(as.numeric(node3[i, 8]))*100,color='red' )
}
map3
map <- get_map(location = c(lon = 103.7662684, lat = 1.2869421), zoom=14, source = "google", maptype = 'watercolor')
bus_graph$layout=cbind(V(bus_graph)$Longitude,V(bus_graph)$Latitude)
#arrow formating
a <- arrow(type = "closed", length = unit(.09, "inches"))
#from the below code we are telling ggmap to take ggraph as your base layer with layout as the cordinates.
#Here we get the curvature using geom_edge_arc.
g<-ggmap(map, fullpage = TRUE, base_layer = ggraph(bus_graph)) +
geom_edge_arc(aes(edge_size=E(bus_graph)$centrality_degree),edge_colour = 'red', edge_alpha = 0.5,curvature = 0.2, arrow=a, end_cap = circle(.008, 'inches')) +
geom_node_point(aes(colour = closeness_centrality,size=betweenness_centrality))+ scale_colour_gradientn(colours=rainbow(3))+scale_size_continuous(range = c(1, 10))
plot(g)
