Data preparation

Let’s first extract the relevant data, adding columns with the time-during-the-day and day-of-the-week for each entry involving the two stations under examination.
We also add flags for whether a given activity was in/out, or during the week (M-F) or week end.

sid1 <- "fe2a5f"
sid2 <- "fec8ff"

# subset on station 1
st1 <- subset(data, station1 == sid1 | station2 == sid1)
st1$t1hc <- (hour(st1$time1)*60 + minute(st1$time1))/60.0
st1$t2hc <- (hour(st1$time2)*60 + minute(st1$time2))/60.0
# the default is with Sunday = 1, but I prefer Sunday = 7.
temp <- wday(st1$time1) - 1 ; st1$wday1 <- ifelse(temp == 0, 7, temp)
temp <- wday(st1$time2) - 1 ; st1$wday2 <- ifelse(temp == 0, 7, temp)
st1$flag_day <- ifelse(st1$wday1 <= 5, "M-F", "WeekEnd")

st1$flag_out <- st1$station1 == sid1
st1$flag_in  <- st1$station2 == sid1
st1$flag_both <- (st1$station1 == sid1) & (st1$station2 == sid1)
st1$type <- NA
st1$type[st1$flag_in] <- 'in'
st1$type[st1$flag_out] <- 'out'

st1$ID <- sid1

# subset on station 2
st2 <- subset(data, station1 == sid2 | station2 == sid2)
st2$t1hc <- (hour(st2$time1)*60 + minute(st2$time1))/60.0
st2$t2hc <- (hour(st2$time2)*60 + minute(st2$time2))/60.0
temp <- wday(st2$time1) - 1  ; st2$wday1 <- ifelse(temp == 0, 7, temp)
temp <- wday(st2$time2) - 1  ; st2$wday2 <- ifelse(temp == 0, 7, temp)
st2$flag_day <- ifelse(st2$wday1 <= 5, "M-F", "WeekEnd")

st2$flag_out <- st2$station1 == sid2
st2$flag_in  <- st2$station2 == sid2
st2$flag_both <- (st2$station1 == sid2) & (st2$station2 == sid2)
st2$type <- NA
st2$type[st2$flag_in] <- 'in'
st2$type[st2$flag_out] <- 'out'

st2$ID <- sid2

# add a 'factor' column that would be handy for analysis
st1$ID <- sid1
st2$ID <- sid2

# combine the two data frames
both <- rbind(st1, st2)
both$ID <- as.factor(both$ID)
both$type <- as.factor(both$type)
both$flag_day <- as.factor(both$flag_day)

The usage patterns of these two stations, by hour of the day and by day of the week are summarized in the following simple plots, which show two clear qualitative differences.

Plots : by hour of the day

• Station fe2a5f

  • Usage slowly increases during the day peaking peaking in the late afternoon and a slow-ish decline with significant usage after 8pm and until midnight.

  • Most of the activity is pretty even between incoming and outgoing, with only modest peaks in the early morning (outgoing) and at around 5-6pm (incoming).

• Station fec8ff

  • This station exhibits a very different hourly pattern, that seems to be related to what we may call working hours, with steep rise and decline early in the morning and late afternoon, and very little activity after 8pm.

  • The tentative association with working hours is corroborated by the striking asymmetry of the split between incoming and outgoing trips: morning usage is almost exclusively arriving, while the big peak of usage in the afternoon is almost entirely leaving.

p1 <- ggplot(data = both, aes(x = t1hc, fill = type)) + theme_bw() + xlab("hour of the day")

# density plot
p1 + geom_density(alpha = 0.5) + facet_wrap(~ ID, ncol = 1) + 
        scale_fill_manual(values = c("navyblue", "gold"))

Plots : by day of the week

The by-hour usage patterns seem to suggest that a more leisure-related kind of activity characterizes station fe2acf, while the usage of bikes in/out of station fec8ff is fairly closely associated to office hours and the idea that it is driven by people using the bikes to “commute” to work.

The usage by day of the week adds weight to this idea by showing a “significant” drop in traffic for station fec8ff during the weekend.
There is instead a hint of increase of activity for station fe2a5f during the weekend, that lends further support to the idea that its location may be more residential.

p2 <- ggplot(data = both, aes(x = as.factor(wday1), fill = ID)) + theme_bw() + 
        theme(legend.position="none") +
        xlab("day of the week") 

# histogram w/ free Y-scale
p2 +  geom_bar(stat = "bin", alpha = 0.75, width = 0.8) + 
        # scale_fill_brewer(palette = "Set2") + 
        scale_fill_manual(values = c("navyblue", "gold")) + 
        facet_wrap(~ ID, ncol = 1, scale = "free_y")

Data preparation

# station IDs
sid1 <- "8f0f64" 
sid2 <- "4a4b61"

# define valid time interval
bound1 <- force_tz(ymd_hms("2013-10-30 00:00:00"), "US/Eastern")
bound2 <- force_tz(ymd_hms("2013-10-30 23:59:59"), "US/Eastern")
interval <- as.interval(86399, bound1)

We want to work with a new data frames, for lighter handling.
For further convenience, we shift columns to align time and stations of the potential arrive-depart intervals to be on the same row.

T_at_St <- data[, c(1, 3, 5, 2, 4)]
colnames(T_at_St) <- c("bike", "St_in", "T_in", "St_out", "T_out")

nl <- nrow(T_at_St)
jj <- c( seq(2, nl), 1 )

# shifting columns to align time and stations of the potential arrive-depart intervals
T_at_St$St_out <- T_at_St$St_out[jj]
T_at_St$T_out  <- T_at_St$T_out[jj]
T_at_St <- T_at_St[-nl, ]

# duration of stay
T_at_St$T_park <- difftime(T_at_St$T_out, T_at_St$T_in)
T_at_St$flag.station <- T_at_St$St_in == T_at_St$St_out

oct30 <- T_at_St[T_at_St$T_in %within% interval | T_at_St$T_out %within% interval, ]
st1 <- subset(oct30, (St_in == sid1 | St_out == sid1) )
st2 <- subset(oct30, (St_in == sid2 | St_out == sid2) )

Filtering

We define a few flags to summarize the “quality” of the data.

• flag.in, flag.out, flag.within : mark if the station was the arrival or departure station of the interval, or both.
• flag.check : summarizes the above flags in one number.
• flag.T.in : TRUE/FALSE depending on whether or not the begin time of the interval (the arrival time) is contained in the valid time window.
• flag.T.out : TRUE/FALSE depending on whether or not the end time of the interval (the departure time) is contained in the valid time window.
• flag.T.within : TRUE is the stay is fully contained within the valid time window.

st1$flag.in  <- st1$St_in  == sid1
st1$flag.out <- st1$St_out == sid1
st1$flag.within <- st1$flag.in & st1$flag.out
st1$check <- 1*st1$flag.in + 2*st1$flag.out + 4*st1$flag.within

st1$flag.T.in  <- st1$T_in  >= bound1
st1$flag.T.out <- st1$T_out <= bound2
st1$flag.T.within <- st1$flag.T.in & st1$flag.T.out 
st1$check.T <- 1*st1$flag.T.in + 2*st1$flag.T.out + 4*st1$flag.T.within

st1.condition1 <- st1$check == 7 & st1$check.T == 7 
st1.cond.xbad  <- (st1$check + st1$check.T) == 3
st1.cond.clean <- st1$check == 7 & st1$check.T == 7

st2$flag.in  <- st2$St_in  == sid2
st2$flag.out <- st2$St_out == sid2
st2$flag.within <- st2$flag.in & st2$flag.out
st2$check <- 1*st2$flag.in + 2*st2$flag.out + 4*st2$flag.within

st2$flag.T.in  <- st2$T_in  >= bound1
st2$flag.T.out <- st2$T_out <= bound2
st2$flag.T.within <- st2$flag.T.in & st2$flag.T.out 
st2$check.T <- 1*st2$flag.T.in + 2*st2$flag.T.out + 4*st2$flag.T.within

st2.condition1 <- st2$check == 7 & st2$check.T == 7 
st2.cond.xbad  <- (st2$check + st2$check.T) == 3
st2.cond.clean <- st2$check == 7 & st2$check.T == 7

best1 <- st1[st1.cond.clean, ]
best2 <- st2[st2.cond.clean, ]

Given the incompleteness of the dataset (or the “manual” re-organization of bikes among stations) there are a variety of combinations of the above conditions, some of which make some intervals unusable.

Some combinations of In/Out stations and times (e.g. out of the bounds of 10/30) that would give open intervals, or intervals with undefined left/right boundaries could be handled by making some assumptions about how long a bicycle could have been there or could have stayed there.

For simplicity I kept only intervals with the same In/Out station and with In/Out times fully within the chosen date.

Statistics

Having selected only the cleanest possible set of stays, for each of them we mark with hour-marks it encompasses, encoding it into a vector.

We then use these indices to mark times with a 0/1 in a matrix comprising all stays for each station.

# * The +1 is to shift the number to match a valid index value.
# * Hours are 0-23, but indices are 1-24.
st1.x1 <- hour(best1$T_in) + 1 + 1
st1.x2 <- hour(best1$T_out) + 1
st1.x1 <- ifelse(st1.x1 > 24, 23, st1.x1)
st1.x2 <- ifelse(st1.x2 > 24, 23, st1.x2)

# initialize 0/1 matrix
presence1 <- matrix(0, nrow = nrow(best1), ncol = 24)
for(j in 1:nrow(best1)) {
    presence1[j, st1.x1[j]:st1.x2[j]] <- 1
}
byHour1 <- colSums(presence1)

st2.x1 <- hour(best2$T_in) + 1 + 1
st2.x2 <- hour(best2$T_out) + 1
st2.x1 <- ifelse(st2.x1 > 24, 23, st2.x1)
st2.x2 <- ifelse(st2.x2 > 24, 23, st2.x2)

# initialize 0/1 matrix
presence2 <- matrix(0, nrow = nrow(best2), ncol = 24)
for(j in 1:nrow(best2)) {
    presence2[j, st2.x1[j]:st2.x2[j]] <- 1
}
byHour2 <- colSums(presence2)

df1 <- data.frame(ID = sid1, hour = 0:23, number = byHour1)
df2 <- data.frame(ID = sid2, hour = 0:23, number = byHour2)
df <- rbind(df1, df2)
df$ID <- as.factor(df$ID)

df.wide <- data.frame(hour = 0:23, station1 = byHour1, station2 = byHour2)
colnames(df.wide) <- c("hour", sid1, sid2)
# write.csv(df.wide, file = "presence_by_hour.csv", row.names = FALSE, quote = FALSE)

print(df.wide, row.names = FALSE)
#  hour 8f0f64 4a4b61
#     0      0      0
#     1      0      2
#     2      0      2
#     3      0      1
#     4      0      1
#     5      0      1
#     6      0      2
#     7      1      7
#     8      4     24
#     9     28     34
#    10     27     44
#    11     27     43
#    12     27     43
#    13     28     40
#    14     29     36
#    15     28     38
#    16     26     34
#    17     16     29
#    18      2      6
#    19      1      6
#    20      4      7
#    21      4      4
#    22      1      1
#    23      0      1

Plot

p3 <- ggplot(data = df, aes(x = hour, y = number, fill = ID)) + theme_bw() + theme(legend.position="none") +
        xlab("hour")

# histogram
p3 + geom_bar(stat = "identity", alpha = 0.5, width = 0.8) + 
        scale_fill_brewer(palette = "Set1") + 
        facet_grid(ID ~ .)