library(readr)
library(dplyr)
library(ggplot2)
library(leaflet)
library(tidyr)
library(lubridate)
library(stringi)
library(RColorBrewer)
Junction elements:
Comparing the tasks and performance of junctions is a great way compare and improve traffic situations. However, the roads and junctions are a network and therefore inter-connected. Hence the poor performance of a junction may be due to the poor performance of a nearby system (road, junction, etc.) and not because of its own components/configuration. However, if a junction has excellent performance we can pride its components and configuration. Statements:
This means that you can’t simply replace the components of a bad performing junction with components of a similar (in task) components of a good performing junction. Because the bad performance may be due to the system. Hence, you first have to analyse the traffic as a whole.
I would suggest that analysis is done to inspect the origin of bad traffic. We should assume that this is either at a junction, or at a place of an accident or similar unpredicted event.
If we can identify the junction that first starts clogging and slow down overal traffic than we can start by applying junction analysis on that junction. Identify its task, identify other junctions with a similar task, and do a perfomance analysis.
In the following exploration we have skipped a proper analysis of localising the origin of traffic problems and we started to analyse junction A13 N470 because we were told it is a troublesome junction.
I’ve requested and gathered a lot of data from the NDW open data. A request can be submitted here. You can receive intensity or speed data for 10 locations per request and data is averaged to hourly blocks, an average over all the days you request. I’ve focused on Thursday data from all Thursdays in April and May 2016.
Shamelessly, I preprocessed the data in Excel. Although it wasn’t really processing, more moving the data together in neat tables. Both the ‘raw’ sheets and the combined tables are available on Github, if not than remind me.
I’ve collected speed and intensity data from ten locations around the A13 N470 junction.
Hour_A13N470_1 <- read_csv("NDW Data/A13/Hourly Data A13 N470 1.0.csv")
Info_A13N470_1 <- read_csv("NDW Data/A13/Info A13 N470 1.0.csv")
data <- tbl_df(Hour_A13N470_1)
Rename columns for easier manipulation and plotting.
names(data)[c(1,6,11)] <- c("Time", "IpH", "SpH")
Create a map of the locations at which there are intensity and speed measurements.
leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addCircleMarkers(data = data,
lat = ~ Lat,
lng = ~ Lon,
popup= ~as.character(ID),
radius = 3)
Plot the intensity versus speed of each of the ten measurement locations.
data %>%
ggplot(aes(x = IpH, y = SpH, col = factor(ID))) +
geom_point(alpha = 0.4) +
geom_smooth(se = FALSE)
The plot shows a clear separation by speed of the roads.
The N470 seems equally busy on both sides but the traffic is much faster on Kruithuisweg (ID 1,2,3 West of the A13), than on the namesless N470 (ID 4,5 East of the A13). I assume that’s because they are different road types, with different speed limits.
The A13 also shows differences. South-East bound towards Rotterdam seems to drop a lot of speed when it gets really busy, while North-West bound towards Den Haag doesn’t drop in speed for the same intensity.
We are now left with six directions with ID 1, 2, 4, 5, 7 and 9. We will add four more directions on the A13, in both directions of the A13 outside of both N470 exits. An image will illustrate all directions.
rm(data)
rm(Hour_A13N470_1)
rm(Info_A13N470_1)
First we filter the data we are interested in, with ID numbers: 1, 2, 4, 5, 7, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19.
Hour_A13N470_2 <- read_csv("NDW Data/A13/Hourly Data A13 N470 1.1.csv")
Info_A13N470_2 <- read_csv("NDW Data/A13/Info A13 N470 1.1.csv")
data <- tbl_df(Hour_A13N470_2 %>% filter(ID %in% c(1, 2, 4, 5, 7, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19)))
info <- tbl_df(Info_A13N470_2 %>% filter(ID %in% c(1, 2, 4, 5, 7, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19)))
Rename columns for easier manipulation and plotting.
names(data)[c(1,6,11)] <- c("Time", "IpH", "SpH")
Create a map of the locations at which there are intensity and speed measurements.
leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addCircleMarkers(data = data,
lat = ~ Lat,
lng = ~ Lon,
popup= ~as.character(ID),
radius = 3)
The newly added data at ID 11 to 19 probably have a lot of similarity. Expect ID’s to be similar between 11 and 12, 16 and 17, 13 and 14 and 15, and 18 and 19.
data %>%
filter(ID %in% 11:19) %>%
ggplot(aes(x = IpH, y = SpH, col = factor(ID))) +
geom_point(alpha = 0.4) +
geom_smooth(se = FALSE)
The expectation was correct hence we can drop some ID numbers. Based on map locality I pick 12, 14, 17, and 18 to remain.
That results in a map:
leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addCircleMarkers(data = data,
lat = ~ Lat,
lng = ~ Lon,
popup= ~as.character(ID),
radius = 3)
With corresponding intensity-speed plot:
data <- data %>% filter(ID %in% c(1, 2, 4, 5, 7, 9, 12, 14, 17, 18))
info <- info %>% filter(ID %in% c(1, 2, 4, 5, 7, 9, 12, 14, 17, 18))
data %>%
ggplot(aes(x = IpH, y = SpH, col = factor(ID))) +
geom_point(alpha = 0.4) +
geom_smooth(se = FALSE)
We might not have the full data of what happens around the junction but we extrapolate more info from the data we have. We try to understand Kruithuisplein, the actual junction with traffic lights sitting on top of the A13.
Traffic light locations with labels.
We would love the traffic light data for this junction, but without it we’ll try to undestand as much as possible with the following breakdown:
Hence, in terms of intensity
Separate intensity data and generate the traffic light data.
TL <- data %>%
select(c(1,2,6)) %>%
mutate(ID = paste0("ID_", ID)) %>%
spread(ID, IpH)
TL <- TL %>%
mutate(TL_1 = ID_5,
TL_3 = ID_14 - ID_7,
TL_5 = ID_2,
TL_7 = ID_17 - ID_9,
TL_2p = TL_3 - ID_4,
TL_4p = TL_5 - (ID_18 - ID_9),
TL_6p = TL_7 - ID_1,
TL_8p = TL_1 - (ID_12 - ID_7))
Plot the traffic light data. Note that half the traffic lights are calculated as greater than (they are plotted with dotted lines).
TL %>%
select(c(1,12:19)) %>%
gather("TrafficLight", "I", 2:9) %>%
mutate("Time" = as.numeric(substr(Time, 1, 2)),
"Actual" = grepl("p", TrafficLight)) %>%
ggplot(aes(x = Time, y = I, col = TrafficLight)) +
geom_point(alpha = 0.5) +
geom_path(aes(linetype = Actual))
Instead of looking at Kruithuisplein we now explore the ten important flows.
intensity <- data %>%
select(c(1,2,6)) %>%
mutate(ID = paste0("ID_", ID)) %>%
spread(ID, IpH) %>%
transmute(
Time = as.numeric(substr(Time, 1, 2)),
Junction = "A13 N470",
W_O = ID_1,
W_I = ID_2,
E_O = ID_4,
E_I = ID_5,
S_T = ID_7,
N_T = ID_9,
S_I = ID_14 - ID_7,
N_O = ID_12 - ID_7,
N_I = ID_17 - ID_9,
S_O = ID_18 - ID_9)
Plot the intensity ‘signature’ of the junction.
intensity %>%
gather("Dir", "I", 3:12) %>%
mutate(Dir2 = factor(substr(Dir, 1, 1)),
Dir3 = factor(substr(Dir, 3, 3))) %>%
ggplot(aes(x = Time, y = I, group = Dir, col = Dir2)) +
geom_point(alpha = 0.5, aes(shape = Dir3), size = 2) +
geom_path()
To explain the plot:
So the red line with circles is the road on the east entering the junction (ID5). Or the green line with squares is the road on the north going underneath the junction (ID9).
My apologies for the complication. When we compare with other junctions it might makes sense and be useful.
Before we head on we may as well quickly gather the velocity data for the similar directions. We check the velocities only on the roads where people are driving. We check traffic from the south driving through (S_T is ID_7), from the north going through (N_T is ID_9), and West and East coming in and going out.
velocity <- data %>%
select(c(1, 2, 11)) %>%
mutate(ID = paste0("ID_", ID)) %>%
spread(ID, SpH) %>%
transmute(
Time = as.numeric(substr(Time, 1, 2)),
Junction = "A13 N470",
W_O = ID_1,
W_I = ID_2,
E_O = ID_4,
E_I = ID_5,
S_T = ID_7,
N_T = ID_9)
And plot.
velocity %>%
gather("Dir", "I", 3:8) %>%
mutate(Dir2 = factor(substr(Dir, 1, 1)),
Dir3 = factor(substr(Dir, 3, 3))) %>%
ggplot(aes(x = Time, y = I, group = Dir, col = Dir2)) +
geom_point(alpha = 0.5, aes(shape = Dir3), size = 2) +
geom_smooth(se = FALSE, span = 0.5)
Most importantly is the dip in Southbound direction in the afternoon. This is it’s worst performing part. We can analyse the southbound average speed around the junction and see if the performance increases, decreases or remains the same as traffic passes this junction. Passing location ID’s 16, 17, 9, 18, 19 in that order.
Select data:
SB_vel <- tbl_df(Hour_A13N470_2) %>%
select(c(1, 2, 11)) %>%
filter(ID %in% c(16, 17, 9, 18, 19)) %>%
transmute(Time = as.numeric(substr(`Uur op de dag`, 1, 2)),
ID = factor(ID, levels = c(16, 17, 9, 18, 19)),
SpH = `Average Speed per hour`)
Plot:
ggplot(data = SB_vel, aes(x = Time, y = SpH, col = ID)) +
geom_point() +
geom_smooth(se = FALSE, span = 0.5, alpha = 0.3)
Hour_A13N470_3 <- read_csv("NDW Data/A13/Hourly Data A13 N470 1.2.csv")
SB_vel <- tbl_df(Hour_A13N470_3) %>%
select(c(1:5, 11)) %>%
filter(ID %in% c(20, 21, 22, 23, 16, 17, 9, 18, 19, 24, 25, 26, 27, 28)) %>%
transmute(Time = as.numeric(substr(`Uur op de dag`, 1, 2)),
ID = factor(ID, levels = c(20, 21, 22, 23, 16, 17, 9, 18, 19, 24, 25, 26, 27, 28)),
SpH = `Average Speed per hour`,
Lat = Lat,
Lon = Lon,
Dis = as.numeric(stri_sub(Location, -4, -1)))
colours <- colorRampPalette(brewer.pal(8,"Dark2"))(8)
SB_vel %>%
filter(Time %in% c(7:10, 16:19)) %>%
arrange(ID, Time) %>%
ggplot(aes(x = Dis, y = SpH, col = factor(Time))) +
geom_point(position = position_jitter(0.1), alpha = 0.5) +
geom_line(stat = "smooth", method = "loess", span = 0.3, size = 2, alpha = 0.7) +
geom_vline(xintercept = c(7.2, 9.1, 11.1, 16.7), color = "darkgrey", linetype = 2) +
geom_text(aes(x = 7.2, label = "Brassenkade", y = 50), colour="darkgrey", angle = 90, vjust = 1.2) +
geom_text(aes(x = 9.1, label = "Oostporweg", y = 50), colour="darkgrey", angle = 90, vjust = 1.2) +
geom_text(aes(x = 11.1, label = "N470", y = 50), colour="darkgrey", angle = 90, vjust = 1.2) +
geom_text(aes(x = 16.7, label = "N209", y = 50), colour="darkgrey", angle = 90, vjust = 1.2) +
scale_color_manual("Hour of \nthe day", values = colours) +
labs(title = "Hourly averaged speed southbound along the A13 - Grey lines indicate junctions", x = "Distance by 'hectometer paaltjes' (km)", y = "Average speed for the following hour (km/hr)")
geom_vline(xintercept=200, colour=“grey”) + geom_text(aes(x=200, label=“the strong cars”, y=20), colour=“blue”, angle=90, vjust = 1.2, text=element_text(size=11))+ geom_text(aes(x=200, label=“the weak cars”, y=20), colour=“red”, angle=90, vjust = -1, text=element_text(size=11))
leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addCircleMarkers(data = SB_vel,
lat = ~ Lat,
lng = ~ Lon,
popup= ~as.character(ID),
radius = 3)
Now that we have a profile of intenstiy and speed across the A13 N470 junction we look for junctions with similar intensity. After that we can compare the speed data of both to see which junction handles the same amount of traffic more effectively.
Scrolling the map of data points on ndwopendata.kxa.nl I’ve identified about 20 junction that could match A13-N470.
The list is that has been identified.
I’ve collected speed and intensity data from twenty locations around the A8 N516 junction (Coenplein).
Hour_A8N516_1 <- read_csv("NDW Data/A8-N516/Hourly Data A8 N516 1.0.csv")
Info_A8N516_1 <- read_csv("NDW Data/A8-N516/Info A8 N516 1.0.csv")
coen <- tbl_df(Hour_A8N516_1)
Rename columns for easier manipulation and plotting.
names(coen)[c(1,6,11)] <- c("Time", "IpH", "SpH")
Create a map of the locations at which there are intensity and speed measurements.
# Jitter some geo coords because the overlap completely
coen[coen$ID %in% c(32, 34, 50),"Lat"] <- coen[coen$ID %in% c(32, 34, 50),"Lat"] + 0.0001
leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addCircleMarkers(data = coen,
lat = ~ Lat,
lng = ~ Lon,
popup= ~as.character(ID),
radius = 3)
Plot the intensity versus speed of each of the ten measurement locations.
coen %>%
ggplot(aes(x = IpH, y = SpH, col = factor(ID))) +
geom_point(alpha = 0.4) +
geom_smooth(se = FALSE)
Explore the overlap and interconnectedness of values at different locations.
Southbound A8: 36, 35, 33, 34, 41, 42, 43, 51, 52, 53
# Intensity taking the exit on Southbound A8
S1 <- (filter(coen, ID == 36) %>% select(IpH) - filter(coen, ID == 35) %>% select(IpH))[[1]]
S1
## [1] 71 39 29 32 71 237 673 767 747 566 508 532 548 587 572 548 595
## [18] 550 415 417 293 267 210 136
# 35 turns into 33 and 34
S3 <- filter(coen, ID == 33) %>% select(IpH) + filter(coen, ID == 34) %>% select(IpH)
(S3 - filter(coen, ID == 35) %>% select(IpH))[[1]]
## [1] 14 9 7 11 14 34 35 46 31 71 56 59 82 68 80 73 74 79 61 63 48 54 47
## [24] 27
S3 <- S3[[1]]
# 34 becomes 43 becomes
(filter(coen, ID == 34) %>% select(IpH) - filter(coen, ID == 43) %>% select(IpH))[[1]]
## [1] 17 25 49 130 288 711 914 694 840 264 244 119 143 233 194 3 41
## [18] 120 -63 53 87 127 121 45
# (41 + 42) becomes (51 + 52)
S4 <- (filter(coen, ID == 41) %>% select(IpH) - filter(coen, ID == 51) %>% select(IpH))[[1]]
S5 <- (filter(coen, ID == 42) %>% select(IpH) - filter(coen, ID == 52) %>% select(IpH))[[1]]
S4 + S5
## [1] 10 5 7 8 15 44 48 54 18 12 27 34 7 26 19 31 40
## [18] 39 35 -11 12 16 13 8
# 43 becomes 53
S6 <- (filter(coen, ID == 43) %>% select(IpH) - filter(coen, ID == 53) %>% select(IpH))[[1]]
S6
## [1] 1 0 0 1 1 10 14 -2 19 -30 5 -1 0 0 2 2 4
## [18] 4 1 -1 2 1 2 1
# 41 + 42 + 43 equals 51 + 52 + 53
S4 + S5 + S6
## [1] 11 5 7 9 16 54 62 52 37 -18 32 33 7 26 21 33 44
## [18] 43 36 -12 14 17 15 9
# 41 + 42 + 43 - (33 + 34) is entering Southbound
S7 <- (filter(coen, ID == 41) %>% select(IpH) + filter(coen, ID == 42) %>% select(IpH) + filter(coen, ID == 43) %>% select(IpH) - S3)[[1]]
S7
## [1] 95 57 45 46 95 292 747 1104 1065 817 713 697 792 808
## [15] 789 796 868 917 596 492 385 371 299 188
We want to keep southbound S1 (out), S3 (through) and S7 (in) to plot against A13 N470 data.
Northbound A8: 55, 54, 57, 56, 46, 44, 48
# Intenstiy entering the A8 northbound
# 55 becomes 46
(filter(coen, ID == 55) %>% select(IpH) - filter(coen, ID == 46) %>% select(IpH))[[1]]
## [1] 1 0 0 0 -1 NA NA NA NA NA NA NA NA 2 2 12 16
## [18] 13 -26 -2 NA NA NA NA
# 56 becomes 48
(filter(coen, ID == 56) %>% select(IpH) - filter(coen, ID == 48) %>% select(IpH))[[1]]
## [1] -2 -2 0 0 2 4 1 -1 -5 -5 -6 7 -6 0 4 -3 31
## [18] 3 -49 -15 -9 -3 -6 -10
# (46 + 48) becomes 47 becomes 40. (don't trust 47 with it's missing data)
(filter(coen, ID == 46) %>% select(IpH) + filter(coen, ID == 48) %>% select(IpH) - filter(coen, ID == 47) %>% select(IpH))[[1]]
## [1] 1 1 2 1 2 NA NA NA NA NA NA NA NA 26 55 124 129
## [18] 89 -4 37 NA NA NA NA
(filter(coen, ID == 40) %>% select(IpH) - filter(coen, ID == 47) %>% select(IpH))[[1]]
## [1] 0 -1 3 1 1 0 -4 -2 2 -6 4 2 -5 8 25 62 2
## [18] -33 19 25 1 -2 1 1
# 54 becomes 44 becomes 45 becomes 37
(filter(coen, ID == 54) %>% select(IpH) - filter(coen, ID == 44) %>% select(IpH))[[1]]
## [1] -5 -6 -1 -5 -3 -4 -2 3 0 -6 -3 -5 -1 -2 -5 2 2 -3 -7 -6 -5 -4 -7
## [24] -7
(filter(coen, ID == 44) %>% select(IpH) - filter(coen, ID == 45) %>% select(IpH))[[1]]
## [1] -1 0 -3 1 -1 2 10 10 8 8 16 16 15 11 18 10 13 9 0 10 6 7 6
## [24] -2
(filter(coen, ID == 45) %>% select(IpH) - filter(coen, ID == 37) %>% select(IpH))[[1]]
## [1] -6 -1 -1 -1 0 -1 -6 -3 -16 -13 -15 -15 -16 -13 -21 -17 -12
## [18] -19 -18 -11 -10 -9 -10 -7
# (37 + 40) becomes (38 exit traffic)
N2 <- filter(coen, ID == 38) %>% select(IpH)
N1 <- (filter(coen, ID == 37) %>% select(IpH) + filter(coen, ID == 40) %>% select(IpH) - N2)[[1]]
N1
## [1] 129 74 54 46 59 129 393 657 778 639 634 651 709 746
## [15] 834 986 1003 1005 794 593 468 424 412 290
N2 <- N2[[1]]
# (38 + entering traffic becomes 39)
N3 <- (filter(coen, ID == 39) %>% select(IpH) - filter(coen, ID == 38) %>% select(IpH))[[1]]
N3
## [1] 67 35 20 21 30 65 186 340 393 351 397 431 470 501 587 707 936
## [18] 981 719 399 305 277 221 132
We want to keep northbound N1 (out), N2 (through), N3 (in).
intensity_coen <- data_frame(
Time = c(0:23),
Junction = rep("A8 N516", 24),
W_O = coen[coen$ID == 32, "IpH"][[1]],
W_I = coen[coen$ID == 31, "IpH"][[1]],
E_O = coen[coen$ID == 49, "IpH"][[1]],
E_I = coen[coen$ID == 50, "IpH"][[1]],
S_T = N2,
N_T = S3,
S_I = N1,
N_O = N3,
N_I = S1,
S_O = S7)
Plot the intensity ‘signatures’ of two junctions the junction.
bind_rows(intensity, intensity_coen) %>%
gather("Dir", "I", 3:12) %>%
mutate(Dir2 = factor(substr(Dir, 1, 1)),
Dir3 = factor(substr(Dir, 3, 3)),
Junction = factor(Junction)) %>%
ggplot(aes(x = Time, y = I, group = Dir, col = Dir2)) +
geom_point(alpha = 0.5, aes(shape = Dir3), size = 2) +
geom_path() +
facet_grid(. ~ Junction) +
labs(title = "Junction Intensity Profile", x = "Time of day", y = "Intensity, cars per hour")
They have reasonably similar signatures hence we want to plot their velocities. First generate the velocity data for A8 N516.
velocity_coen <- coen %>%
select(c(1, 2, 11)) %>%
mutate(ID = paste0("ID_", ID)) %>%
spread(ID, SpH) %>%
transmute(
Time = as.numeric(substr(Time, 1, 2)),
Junction = "A8 N516",
W_O = ID_32,
W_I = ID_31,
E_O = ID_49,
E_I = ID_50,
S_T = ID_38,
N_T = ID_35)
And plot.
bind_rows(velocity, velocity_coen) %>%
gather("Dir", "I", 3:8) %>%
mutate(Dir2 = factor(substr(Dir, 1, 1)),
Dir3 = factor(substr(Dir, 3, 3)),
Junction = factor(Junction)) %>%
ggplot(aes(x = Time, y = I, group = Dir, col = Dir2)) +
geom_point(alpha = 0.5, aes(shape = Dir3), size = 2) +
geom_smooth(se = FALSE, method = "loess", span = 0.4) +
facet_grid(. ~ Junction) +
labs(title = "Junction Speed Profile", x = "Time of day", y = "Average hourly speeds (km/hr)")
You can see that the velocity dips are pretty equal, but note that the the traffic from the North going through (green) is equally busy in the afternoon in both places and the at the A8 N516 the average speed remains above 100 km per hour.
First we start with the usual. Read the data.
Hour_A67N634_1 <- read_csv("NDW Data/A67-N634/Hourly Data A67 N634 1.0.csv")
Info_A67N634_1 <- read_csv("NDW Data/A67-N634/Info A67 N634 1.0.csv")
geldrop <- tbl_df(Hour_A67N634_1)
Rename columns for easier manipulation and plotting.
names(geldrop)[c(1,6,11)] <- c("Time", "IpH", "SpH")
Create a map of the locations at which there are intensity and speed measurements.
# Jitter some geo coords because the overlap completely
geldrop[geldrop$ID %in% c(62, 67, 80, 73),"Lon"] <- geldrop[geldrop$ID %in% c(62, 67, 80, 73),"Lon"] + 0.0001
geldrop[geldrop$ID %in% c(65, 78, 81),"Lon"] <- geldrop[geldrop$ID %in% c(65, 78, 81),"Lon"] - 0.0001
geldrop[geldrop$ID %in% c(63, 77, 74, 71),"Lat"] <- geldrop[geldrop$ID %in% c(63, 77, 74, 71),"Lat"] + 0.00005
geldrop[geldrop$ID %in% c(84, 79, 64),"Lat"] <- geldrop[geldrop$ID %in% c(84, 79, 64),"Lat"] - 0.00005
leaflet() %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addCircleMarkers(data = geldrop,
lat = ~ Lat,
lng = ~ Lon,
popup= ~as.character(ID),
radius = 3)
Plot the intensity versus speed of each of the ten measurement locations.
geldrop %>%
ggplot(aes(x = IpH, y = SpH, col = factor(ID))) +
geom_point(alpha = 0.4) +
geom_smooth(se = FALSE)