Implementation of an analysis tool for travel mode detection in R

a rule-based approach

Abstract

With the introduction of smartphones that capture GPS data of their users’ locations, the volume of mobility data is continuously growing. One of the challenges in handling mobility data is determining the mode of transportation. Therefore, it is essential to develop methods to efficiently and accurately classify movements. In this study, a rule-based method is applied to assign the travel mode (Walk, Bike, Car, Bus, Tram, Train, and Boat) to a spatial dataset collected with the POSMO app. Various thresholds based on background layers and movement parameters are applied and the method is evaluated using a dataset which contains manually corrected travel modes. Overall, an accuracy of 60.9% was achieved. Using background information, the classification of train and boat led to convincing results, whereas it did not work well for bus and tram. Here, other criteria need to be added. Further, use of speed and acceleration have proven to be useful for the distinction of walk, bike and car segments. Overall, this study established a base for a travel mode detection procedure and tough accuracy did not yet fit the standards, the procedure implemented is a promising starting point for further development.

1 Introduction

Mobility data play a crucial role for traffic planners, public transport providers, and authorities as they serves as the basis for transportation modeling and optimization (Nitsche et al., 2014). In the past, these data were primarily collected through surveys, but today the widespread use of smartphones with GPS capabilities enables data collection on a much larger scale (Wu et al.,2016). However, GPS observations alone provide only geometric and temporal data, and additional information must be extracted if needed (Zhang et al., 2012). One particularly important attribute that requires further analysis is travel mode (Sadeghian et al., 2021; Wu et al., 2016).

In this study, we aim to develop a data science method in R for travel mode detection using GPS data from the POSMO app. Although the POSMO app has built-in in-time detection, it is prone to errors, most likely due to the lack of pre-processing. Pre-processing is a crucial step as unfiltered data often contains various inaccuracies (Wu et al., 2016). Therefore, our goal is to implement a travel mode detection method that can be applied to the data after the acquisition phase, thus enabling pre-processing.

Numerous methods have been explored for travel mode detection from GPS data, each with its own advantages and disadvantages (Sadeghian et al., 2021; Wu et al., 2016). Sadeghian et al. (2021) provided a comprehensive review of these methods, noting a growing interest in the use of machine learning algorithms, particularly unsupervised or deep learning algorithms. Consequently, several recent studies have implemented such techniques (Li et al., 2020; Markos & Yu, 2020; Yu, 2021). These algorithms are capable of handling large amounts of data and achieving highly accurate clusters. (Sadeghian et al., 2021; Wu et al., 2016). However, rule-based methods offer the advantage of including additional GIS layers into the analysis, allowing for the distinction of up to 12 transport modes, which is currently difficult to achieve with machine learning alone (Sadeghian et al., 2021). Hence, although rule-based methods may have lower overall accuracy, they offer higher precision, as many machine learning studies only categorize limited number of modes (Sadeghian et al., 2021; Wu et al., 2016; Nitsche et al., 2014). However, rule-based methods rely on a prior understanding and manually defined rules, making them more time-consuming and less transferable (Sadeghian et al., 2021). Consequently, the choice of method strongly depends on the specific objectives of the study.

Considering POSMO’s real-life applications, the ability to distinguish as many travel modes as possible is our primary focus. Therefore, this research project aims to implement a rule-based data science method in R to identify travel modes from mobile GPS data. We will assess all travel modes used in our training data, including walk, bike, car, bus, train, tram, and boat. Ultimately, we aim to answer the following research questions:

1. How can a basic analysis tool for travel mode detection from GPS data be implemented in R and what accuracy is achieved with it?

2. What are the most effective criteria, features and thresholds for the detection of different travel modes?

2 Material and Methods

2.1 Datasets & Conceptual Model

##loading our own dataset:

data <- read_delim("posmo_data/posmo_v3.csv")
head(data) #time is 2h behind
data$datetime <- as.POSIXct(data$datetime, tz = "UTC-2")  #setting time right
data <- arrange(data, datetime) #arrange, so that datetime is in the right order

data <- select(data, datetime, transport_mode, lon_x, lat_y) # keep only the  attributes needed

data <- st_as_sf(data, coords = c("lon_x","lat_y"), crs = 4326) |>  #cs is transformed to 2056
  st_transform(2056)

data_coord <- st_coordinates(data) #coordinates are extracted
data <- cbind(data, data_coord) #coordinates are binded in separate columns
head(data) #first look
summary(data)

##loading swissTLM3D:

roads <- read_sf("background_data/swisstlm3d_2023-03_2056_5728.shp/TLM_STRASSEN/swissTLM3D_TLM_STRASSE.shp") |> #loading roads
  select(OBJEKTART, geometry) |> #choosing attributes needed
  filter(OBJEKTART != "Verbindung" & OBJEKTART != "Raststaette" & OBJEKTART != "Dienstzufahrt" & OBJEKTART != "Verbindung" & OBJEKTART != "Zufahrt " & OBJEKTART != "Klettersteig") |> #deleting factorlevels which are not of interest
  st_transform(2056)#set crs to 2056 (get rid of LN02)

rails <- read_sf("background_data/swisstlm3d_2023-03_2056_5728.shp/TLM_OEV/swissTLM3D_TLM_EISENBAHN.shp") |> #loading rails
  select(VERKEHRSMI, geometry) |> #choosing attributes needed
  rename(OBJEKTART = VERKEHRSMI) |> #rename for merging
    st_transform(2056)#set crs to 2056 (get rid of LN02)

boats <- read_sf("background_data/swisstlm3d_2023-03_2056_5728.shp/TLM_OEV/swissTLM3D_TLM_SCHIFFFAHRT.shp") |> #loading boats
   select(OBJEKTART, geometry) |> #choosing attributes needed
     st_transform(2056)#set crs to 2056 (get rid of LN02)

stops <- read_sf("background_data/swisstlm3d_2023-03_2056_5728.shp/TLM_OEV/swissTLM3D_TLM_HALTESTELLE.shp")|> #loading stops
   select(OBJEKTART, geometry) |> #choosing attributes needed
    st_transform(2056) #set crs to 2056 (get rid of LN02)

##loading POSMO-dataset for validation:

val_data <- read_delim("posmo_data/posmo_validation.csv")
head(val_data) #time is probably behind as well...
val_data$datetime <- as.POSIXct(val_data$datetime, tz = "UTC-2") 
val_data <- arrange(val_data, datetime) 

val_data <- select(val_data, datetime, transport_mode, lon_x, lat_y) # keep only the  attributes needed

val_data <- st_as_sf(val_data, coords = c("lon_x","lat_y"), crs = 4326) |>  #cs is transformed to 2056
  st_transform(2056)

val_data_coord <- st_coordinates(val_data) #coordinates are extracted
val_data <- cbind(val_data, val_data_coord) #coordinates are binded in separate columns
head(val_data) #first look
summary(val_data)

To build our method, we used movement data from one of our team members, collected with the POSMO app (Genossenschaft Posmo, 2022) over a 54-day period (18.04.-10.06.2023). During this period, a total of 61’279 data points were gathered with a sampling rate set to 5 seconds. We proceeded with the following attributes:

• datetime
• geometry (X- and Y-Coordinates in CH1903+ LV95)
• transport_mode (manually corrected and validated in POSMO)

Additionally, the swissTLM3D dataset was obtained from swisstopo (2023), whereof the following feature classes and attributes were used:

• TLM_STRASSE (OBJEKTART, geometry)
• TLM_EISENBAHN (VERKEHRSMITTEL, geometry)
• TLM_SCHIFFFAHRT (OBJEKTART, geometry)
• TLM_HALTESTELLE (OBJEKTART, geometry)

A second POSMO-dataset provided by one of the other students was used for validation (see chapter 2.6). These data were collected over a period of 67 days (11.04-16.06.2023) with a sampling rate of 10s, resulting in a total of 46’305 fixes. It was eventually corrected and validated for travel mode in POSMO as well and the same attributes were analyzed further.

For preprocessing, analyses and visualizing our data we used R (v. 4.2.3; R Core Team, 2023) as well as the package ggplot2 (v. 3.4.2; Wickham, 2016).

Following Laube (2014), the movement space was conceptualized as continuous, 2D, and entity-based. Therefore, all datasets were structured as vector data. We modeled the movement as a series of unconstrained, intermittent, time-stamped fixes. The app collects GPS data using the Lagrangian perspective with event-based, active tracking.

2.2 Segmentation & Filtering

##segmentation and filtering:

p1 <- data |> #visualize raw data for example day
  filter(as.Date(datetime) == "2023-04-18") |>
  ggplot(aes(X, Y))+ 
  geom_point()+
  geom_path()+
  coord_equal() +
  theme_classic() +
  ggtitle("(A) raw data")          

data <- data |> #initial segmentation for time gaps > 10s (double the sampling rate)
  mutate(
    timelag = as.numeric(difftime(lead(datetime), datetime, units = "secs")), 
    gap = timelag > 10,
    segment_id = cumsum(gap)
  )

data <- data |> 
  group_by(segment_id) |> #calculating the stepmean with the moving temporal window within each segment
    mutate(
    stepMean = rowMeans(
      cbind(
        sqrt((lag(X, 3) - X)^2 + (lag(Y, 3) - Y)^2),
        sqrt((lag(X, 2) - X)^2 + (lag(Y, 2) - Y)^2),
        sqrt((lag(X, 1) - X)^2 + (lag(Y, 1) - Y)^2),
        sqrt((X - lead(X, 1))^2 + (Y - lead(Y, 1))^2),
        sqrt((X - lead(X, 2))^2 + (Y - lead(Y, 2))^2),
        sqrt((X - lead(X, 3))^2 + (Y - lead(Y, 3))^2)
      )
    ),
    static = if_else(is.na(stepMean) | stepMean < 2, T, F) #define every point with a stepMean of less than 2m as static
  )
  
p2 <- data |> #visualize assignment of static points
  filter(as.Date(datetime) == "2023-04-18") |>
  ggplot(aes(X, Y, color=static))+ 
  geom_point()+
  geom_path(group=1)+
  coord_equal() +
  theme_classic() +
  theme(legend.position = c(0.2,0.2))+
  labs(title = "(B) filtering static points")

rle_id <- function(vec) {
    x <- rle(vec)$lengths
    as.factor(rep(seq_along(x), times = x))
} #function to assign unique ID for each segment

data <- data |> 
   ungroup() |> 
    mutate(segment_id = rle_id(static)) |> #assign new segment ID after every break (static points)
    filter(!static) #remove static points

p3 <- data |> #visualize segmentation
  filter(as.Date(datetime) == "2023-04-18") |>
  ggplot(aes(X, Y, color=segment_id))+ 
  geom_point()+
  geom_path()+
  coord_equal() +
  theme_classic() +
  theme(legend.position = "none")+
  labs(title = "(C) segmentation")
    
    
data <- data |> #remove short segments (less than 300m or 2min)
  group_by(segment_id) |> 
  mutate(steplength = sqrt((X - lead(X, 1)) ^ 2 + (Y - lead(Y, 1)) ^ 2)) |> 
  filter(sum(steplength, na.rm=T)>300) |> 
  filter(sum(timelag, na.rm=T)>120)

p4 <- data |> #visualize final segments
  filter(as.Date(datetime) == "2023-04-18") |>
  ggplot(aes(X, Y, color=segment_id))+ 
  geom_point()+
  geom_path()+
  coord_equal() +
  theme_classic() +
  theme(legend.position = "none")+
  labs(title = "(D) filtering short segments")

In order to assign transport modes to trajectories, the raw data has to be divided into sub-segments representing individual movements first (Wu et al.,2016). An initial segmentation was performed by assigning a new segment ID whenever the time gap between consecutive fixes exceeded double the sampling rate (10s), as the event-based tracking of the POSMO app led to large gaps, which needed to be separated for the further calculations. Then, segmentation was performed following Laube & Purves (2011), were static fixes are classified as those whose average Euclidean distance to other fixes inside a temporal window v is less than some threshold d. We chose 30 seconds (8 fixes) for v and 2 meters for d. Next, sub-segments with a length less than 300m or a time-span of less than 2 minutes were removed as well. The process is visualized for a exemplary day in figure 1. Lastly, segments with a average sinuosity (according to chapter 2.3) greater than 2 were removed as well, as they mostly represented GPS errors, as visible in figure 2. Our resulting dataset contained 295 retained trajectories.

##plot for segmentation and filtering:

grid.arrange(p1,p2,p3,p4, nrow=2, ncol=2)

figure 1: segmentation and filtering of trajectories: Raw data (A) was filtered into static and dynamic points (B), then static points were removed and trajectories segmented according to these breaks (C). Last, short segments (<300m or <2min) were removed as well (D)

2.3 Calculating Movment Parameters

##calculate moving parameters:

data <- data |> 
  group_by(segment_id) |> #group by segment, so moving window starts for every segment again
  mutate(
    speed = { #distance and time passed is calculated with lag/lead of 2 -> 10s in both directions
      step_minus2 <- sqrt((lag(X, 2) - X) ^ 2 + (lag(Y, 2) - Y) ^ 2)
      time_minus2 <- abs(as.numeric(difftime(lag(datetime, 2), datetime, units = "secs")))
      step_plus2 <- sqrt((X - lead(X, 2)) ^ 2 + (Y - lead(Y, 2)) ^ 2)
      time_plus2 <- as.numeric(difftime(lead(datetime, 2), datetime, units = "secs"))
      (step_minus2/time_minus2 + step_plus2/time_plus2) / 2 #average is taken 
    },
    acc = { #change in speed and time is calculated with lag/lead of 2 -> 10s in both directions
      speed_minus2 <- speed - lag(speed, 2) 
      speed_plus2 <- lead(speed, 2) - speed
      time_minus2 <- abs(as.numeric(difftime(lag(datetime, 2), datetime, units = "secs")))
      time_plus2 <- as.numeric(difftime(lead(datetime, 2), datetime, units = "secs"))
      (speed_minus2/time_minus2 + speed_plus2/time_plus2) / 2 #average is taken
      },
    sin = { #length from every point to next point in window of 5 points is calculated
      d_minus1 <- sqrt((lag(X, 1) - X) ^ 2 + (lag(Y, 1) - Y) ^ 2) 
      d_minus2 <- sqrt((lag(X, 2) - lag(X, 1)) ^ 2 + (lag(Y, 2) - lag(Y, 1)) ^ 2)
      d_plus1  <- sqrt((X - lead(X, 1))^2 + (Y - lead(Y, 1))^2)
      d_plus2  <- sqrt((lead(X, 1) - lead(X, 2))^2 + (lead(Y,1) - lead(Y, 2))^2)
      d_straight <- sqrt((lag(X, 2) - lead(X, 2))^2 + (lag(Y,2) - lead(Y, 2))^2) #distance from first to last point of window
      (d_minus1 + d_minus2 + d_plus1 + d_plus2)/d_straight #ratio between total length and shortest distance
    }
  )

##filter with sinuosity:

p5 <- data |> 
  filter(as.Date(datetime) == "2023-04-22") |> #plot before filter
  ggplot(aes(X, Y, color=segment_id))+ 
  geom_point()+
  geom_path()+
  coord_equal() +
  theme_classic() +
  theme(legend.position = "none")+
  labs(title = "(A) before filtering sinuosity")
    
data <- data |> #filter out segments with high sinuosity
  group_by(segment_id) |> 
  filter(mean(sin, na.rm=T)<2) 

p6 <- data |> #plot after filter
  filter(as.Date(datetime) == "2023-04-22") |>
  ggplot(aes(X, Y, color=segment_id))+ 
  geom_point()+
  geom_path()+
  coord_equal() +
  theme_classic() +
  theme(legend.position = "none")+
  labs(title = "(B) after filtering sinuosity")

According to Sadeghian et al. (2011) the most common and powerful variables for transport mode detection are average speed, maximum speed and acceleration. We computed speed according to Laube & Purves (2011) using three fixes located inside a temporal window w, which was set to 20 seconds. Acceleration was calculated based on the same principle and defined as the change in velocity over the change in time. Sinuosity was not used in any of the papers in the reviews, but integrated here nonetheless (Sadeghian et al., 2021; Wu et al., 2016). The calculation was based on Laube & Purves (2011) as well, where it is defined as the ratio between a nominal track length and the line connecting the first and last points in a sampling window w consisting of 5 fixes.

##plot for sinuosity filter:

grid.arrange(p5,p6, nrow=1)

figure 2: after calculation of the moving parameters, segments with a sinuosity higher than 2 were removed as well, shown here for a examplary day before (A) and after the filtering (B)

##calculate summary:
#-> mean, max and min for all three parameters + assignment of POSMO-travel mode

data$transport_mode <- as.factor(data$transport_mode)
data_smry <- data |> 
  group_by(segment_id) |> 
  summarise(
            mean_speed = mean(speed, na.rm=T),
            max_speed = max(speed,na.rm=T),
            min_speed = min(speed,na.rm=T),
            mean_acc = mean(acc, na.rm=T),
            max_acc = max(acc,na.rm=T),
            min_acc = min(acc,na.rm=T),
            mean_sin = mean(sin, na.rm=T),
            max_sin = max(sin,na.rm=T),
            min_sin = min(sin,na.rm=T),
            transport_mode_POSMO = levels(transport_mode)[which.max(table(transport_mode))],
            percentage_tm_POSMO = max(table(transport_mode)) / length(transport_mode) * 100 #column with information what %-of fixes were of the most attributed POSMO-travel mode
)


data_smry_long <- pivot_longer(data_smry, #convert from wide to long for visualization
                           cols = -c("segment_id","transport_mode_POSMO", "percentage_tm_POSMO","geometry"),
                           names_to = "parameter") 

Next the minimum, maximum and average values were computed for the three parameters in every segment. Additionally, every segment was assigned with the adjusted travel mode from POSMO. As the segmentation done here and in POSMO may not be identical, fixes of the same segments may have several different transportation modes from POSMO. Hence, the transportation mode which was attributed to most fixes of a segment was chosen. The different features of the three parameters were visualized based on their transportation mode for exploratory data analysis (figure 3).

##plot for EDA:

data_smry_long |> 
ggplot(aes(x=transport_mode_POSMO, y=value, fill=transport_mode_POSMO))+
  geom_boxplot()+
  facet_wrap(~parameter, scales="free_y")+
  theme_classic()+
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank(),
    legend.position = "bottom")

figure 3: exploratory data analysis of the movement parameters acceleration (acc [m/s²)]), sinuosity (sin) and speed (speed [m/s]), visualized as the maximal (max), minimum (min) and average (mean) value per segment and transport mode

It can be observed, that most transportation modes show similar patterns, or at least have large overlaps. Only train trajectories differ considerably from others regarding speed and acceleration. Hence, we have made the decision to approach the transport mode detection by first embedding the trajectories within their geographic context for a better classification.

2.4 Adding Geographic Context

##spatial join: attribute nearest route

background_data <- rbind(roads, rails, boats) #merge the three datasets together
data <- st_join(data, background_data, join = st_nearest_feature) #perform spatial join with nearest feature

data$OBJEKTART <- as.factor(data$OBJEKTART) #create summary with assignment of nearest feature per segment
data_smry_context <- data |> 
  group_by(segment_id) |> 
  summarise(
            nearest_route = levels(OBJEKTART)[which.max(table(OBJEKTART))],  
            percentage_nearest_route = max(table(OBJEKTART)) / length(OBJEKTART) * 100 
) |> 
  as_tibble() |> 
  select(-c("geometry")) #spatial context is removed for join afterwards

According to Gschwend (2015) movement patterns are mostly quantified on the basis of geometric properties and the arrangement of the fixes, while the geographic environment surrounding the movement is often neglected, however would provide useful semantics insights. Hence, using external information can increase the performance of the used algorithms substantially (Sadeghian et al., 2021). As described in chapter 2.1, we used several feature classes of the swissTLM3D dataset and merged all roads, rails and boat-routes into a background layer. Next, a spatial join was performed, attributing every fix of the POSMO data to its nearest feature, as visualized for an exemplary day in figure 4. Lastly, every segment was assigned the feature, which was attributed to most of the fixes of said segment.

##visualization of spatial join

data_01 <- data |> 
  filter(as.Date(datetime) == "2023-04-18") #select example day

bbox <- st_bbox(data_01) #produce a bounding box with that day
clipped_background <- st_crop(background_data, bbox) #cut background to that bounding box

ggplot()+ #visualize
  geom_sf(data=clipped_background, aes(color=OBJEKTART))+
  geom_sf(data=data_01, alpha=0.4, size=2.5, aes(color=OBJEKTART))+
  theme_classic()+
  theme(legend.position = "bottom")

figure 4: the spatial join of the POSMO data (point) to the nearest feature of roads, rails and boats (lines)

##stop-buffer calculations:

stops <- stops |> 
  rename(stop_type = OBJEKTART)

stops_buffer <- st_buffer(stops, 75) #calculating buffer of 75m around every stop
stop_join <- st_join(data, stops_buffer, join = st_within) #spatial join with POSMO data

first_point <- stop_join |> # a vector is produced, containing the information for the first point of every segment
  group_by(segment_id) |> 
  slice_head() |> #first point is chosen
  pull(stop_type, segment_id) #information about stops is saved

last_point <- stop_join |> #same for last point of every segment
  group_by(segment_id) |> 
  slice_tail() |> 
  pull(stop_type, segment_id)

##merging all parameters:

data_smry_context <- data_smry_context |> 
  cbind(first_point, last_point) 
  
data_class <- left_join(data_smry, data_smry_context, by = "segment_id")

Further, the feature class TLM_HALTESTELLE of swissTLM3D containing bus, train and boat stops was used to compute a buffer of 75 meters around every stop. A spatial join was performed for fixes within these buffers. Next, for every segment the information, if the first and the last point were within a stop-buffer was saved additionally for further analysis.

2.5 Travel Mode Detection

##travel mode detection:

data_class$travel_mode_det <- NA

data_class <- data_class |> 
  mutate(
    travel_mode_det = case_when(
      nearest_route == "Bahn" | max_speed > 45 ~ "Train",
      nearest_route == "Tram" ~ "Tram",
      nearest_route %in% c("Autofaehre", "Personenfaehre") ~ "Boat",
      first_point == "Haltestelle Schiff" & last_point == "Haltestelle Schiff" ~ "Boat",
      (first_point %in% c("Haltestelle Bahn", "Haltestelle Bus")) & (last_point %in% c("Haltestelle Bus", "Haltestelle Bahn")) ~ "Bus",
      max_speed > 12.5 | max_acc > 0.3 ~ "Car",
      mean_speed > 2 | max_speed > 5  ~ "Bike",
      TRUE ~ "Walk"
    )
  )

Zhang et. al (2012) suggest, first separating non-motorized and motorized vehicles and then split up these sub-categories in a second classification step. However, the moving parameters of our training data seemed not to differ considerably enough to do so. Hence, we decided to use the context information for a first classification. In that way, train (closest to train rails or max speed > 45m/s), tram (closest to tram rails), boat (closest to boat line or first and last point in boat stop buffer) and bus (first and last stop in bus or train stop buffer) trajectories were classified. This was done step by step, so that trajectories classified once were not able to be classified again. Accordingly, prediction of car, bike and walk trajectories was left. In order to estimate what thresholds should be applied, these were visualized again (see figure 5). In that way, we decided to classify trajectories with a maximum speed over 12m/s or maximum acceleration over 0.3 m/s² as car segments. Then, segments with a main speed higher than 2m/s or maximal speed over 5m/s were classified as bike rides. And lastly, the remaining segments, which thus had to have a main speed lower than 2m/s were classified as walks.

##plot for EDA filtered for only Car, Bike and Walk:

data_smry_long |> 
  filter(transport_mode_POSMO=="Car" | transport_mode_POSMO =="Bike" | transport_mode_POSMO == "Walk") |> 
ggplot(aes(x=transport_mode_POSMO, y=value, fill=transport_mode_POSMO))+
  geom_boxplot()+
  facet_wrap(~parameter, scales="free_y")+
  theme_classic()+
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    legend.title = element_blank(),
    legend.position = "bottom")

Figure 5: exploratory data analysis of the movement parameters acceleration (acc [m/s²)]), sinuosity (sin) and speed (speed [m/s]) for car, bike and walk segments, visualized as the maximal (max), minimum (min) and average (mean) value per segment and transport mode

2.6 Validation

##filter out cable car as it was not present in our training data and thus model:

val_data <- val_data |> 
  filter(transport_mode != "Cable_Car")

###whole process for validation dataset:

##filtering and segmentation:
      
val_data <- val_data |> #initial segmentation for time gaps > 20s (double the sampling rate)
  mutate(
    timelag = as.numeric(difftime(lead(datetime), datetime, units = "secs")), 
    gap = timelag > 20,
    segment_id = cumsum(gap)
  )

val_data <- val_data |> 
  group_by(segment_id) |> #calculating the stepmean with the moving temporal window within each segment
    mutate(
    stepMean = rowMeans(
      cbind(
        sqrt((lag(X, 2) - X)^2 + (lag(Y, 2) - Y)^2),
        sqrt((lag(X, 1) - X)^2 + (lag(Y, 1) - Y)^2),
        sqrt((X - lead(X, 1))^2 + (Y - lead(Y, 1))^2),
        sqrt((X - lead(X, 2))^2 + (Y - lead(Y, 2))^2)
      )
    ),
    static = if_else(is.na(stepMean) | stepMean < 2, T, F) #define every point with a stepMean of less than 2m as static
  )

val_data <- val_data |> 
   ungroup() |> 
    mutate(segment_id = rle_id(static)) |> #assign new segment ID after every break (static points)
    filter(!static) #remove static points

val_data <- val_data |> #remove short segments (less than 300m or 2min)
  group_by(segment_id) |> 
  mutate(steplength = sqrt((X - lead(X, 1)) ^ 2 + (Y - lead(Y, 1)) ^ 2)) |> 
  filter(sum(steplength, na.rm=T)>300) |> 
  filter(sum(timelag, na.rm=T)>120)


##calculate moving parameters:

val_data <- val_data |> 
  group_by(segment_id) |> 
  mutate(
    speed = { #distance and time passed is calculated with lag/lead of 1 -> 10s in both directions
      step_minus1 <- sqrt((lag(X, 1) - X) ^ 2 + (lag(Y, 1) - Y) ^ 2)
      time_minus1 <- abs(as.numeric(difftime(lag(datetime, 2), datetime, units = "secs")))
      step_plus1 <- sqrt((X - lead(X, 1)) ^ 2 + (Y - lead(Y, 1)) ^ 2)
      time_plus1 <- as.numeric(difftime(lead(datetime, 1), datetime, units = "secs"))
      (step_minus1/time_minus1 + step_plus1/time_plus1) / 2 #average is taken 
    },
    acc = { #change in speed and time is calculated with lag/lead of 1 -> 10s in both directions
      speed_minus1 <- speed - lag(speed, 1) 
      speed_plus1 <- lead(speed, 1) - speed
      time_minus1 <- abs(as.numeric(difftime(lag(datetime, 1), datetime, units = "secs")))
      time_plus1 <- as.numeric(difftime(lead(datetime, 1), datetime, units = "secs"))
      (speed_minus1/time_minus1 + speed_plus1/time_plus1) / 2 #average is taken
      },
    sin = { #length from every point to next point in window of 5 points is calculated
      d_minus1 <- sqrt((lag(X, 1) - X) ^ 2 + (lag(Y, 1) - Y) ^ 2) 
      d_minus2 <- sqrt((lag(X, 2) - lag(X, 1)) ^ 2 + (lag(Y, 2) - lag(Y, 1)) ^ 2)
      d_plus1  <- sqrt((X - lead(X, 1))^2 + (Y - lead(Y, 1))^2)
      d_plus2  <- sqrt((lead(X, 1) - lead(X, 2))^2 + (lead(Y,1) - lead(Y, 2))^2)
      d_straight <- sqrt((lag(X, 2) - lead(X, 2))^2 + (lag(Y,2) - lead(Y, 2))^2) #distance from first to last point of window
      (d_minus1 + d_minus2 + d_plus1 + d_plus2)/d_straight #ratio between total length and shortest distance
    }
  )

##filter with sinuosity:

val_data <- val_data |> #filter out segments with hig sinuosity
  group_by(segment_id) |> 
  filter(mean(sin, na.rm=T)<2)

##calculate summary of moving parameters:

val_data$transport_mode <- as.factor(val_data$transport_mode)
val_data_smry <- val_data |> 
  group_by(segment_id) |> 
  summarise(
            mean_speed = mean(speed, na.rm=T),
            max_speed = max(speed,na.rm=T),
            min_speed = min(speed,na.rm=T),
            mean_acc = mean(acc, na.rm=T),
            max_acc = max(acc,na.rm=T),
            min_acc = min(acc,na.rm=T),
            mean_sin = mean(sin, na.rm=T),
            max_sin = max(sin,na.rm=T),
            min_sin = min(sin,na.rm=T),
            transport_mode_POSMO = levels(transport_mode)[which.max(table(transport_mode))],
            percentage_tm_POSMO = max(table(transport_mode)) / length(transport_mode) * 100
)

##spatial join with nearest feature:

val_data <- st_join(val_data, background_data, join = st_nearest_feature)

val_data$OBJEKTART <- as.factor(val_data$OBJEKTART) #save for every segment
val_data_smry_context <- val_data |> 
  group_by(segment_id) |> 
  summarise(
            nearest_route = levels(OBJEKTART)[which.max(table(OBJEKTART))],  
            percentage_nearest_route = max(table(OBJEKTART)) / length(OBJEKTART) * 100
) |> 
  as_tibble() |> 
  select(-c("geometry"))

##stop-buffer calculations:

val_stop_join <- st_join(val_data, stops_buffer, join = st_within) #spatial join with POSMO data

first_point <- val_stop_join |> # a vector is produced, containing the information for the first point of every segment
  group_by(segment_id) |> 
  slice_head() |> 
  pull(stop_type, segment_id) 

last_point <- val_stop_join |> #same for last point of every segment
  group_by(segment_id) |> 
  slice_tail() |> 
  pull(stop_type, segment_id)

##merging all parameters:

val_data_smry_context <- val_data_smry_context |> 
  cbind(first_point, last_point) 

val_data_class <- left_join(val_data_smry, val_data_smry_context, by = "segment_id")

##travel mode classification:

val_data_class$travel_mode_det <- NA

val_data_class <- val_data_class |> 
  mutate(
    travel_mode_det = case_when(
      nearest_route == "Bahn" | max_speed > 45 ~ "Train",
      nearest_route == "Tram" ~ "Tram",
      nearest_route %in% c("Autofaehre", "Personenfaehre") ~ "Boat",
      first_point == "Haltestelle Schiff" & last_point == "Haltestelle Schiff" ~ "Boat",
      (first_point %in% c("Haltestelle Bahn", "Haltestelle Bus")) & (last_point %in% c("Haltestelle Bus", "Haltestelle Bahn")) ~ "Bus",
      max_speed > 12.5 | max_acc > 0.3 ~ "Car",
      mean_speed > 1.94 | max_speed > 5  ~ "Bike",
      TRUE ~ "Walk"
    )
  )

First, the output of our travel mode detection was compared to the validated POSMO classification and a confusion matrix was rendered. Of course, as our procedure was built up on this data, this was was rather a comparison with the ground truth than a validation. So for a valid predictive validation according to Rykiel (1996), the second dataset was used. As the sampling rate was set at 10s, wherever a temporal window was used, number of fixes were adjusted in order to work with the same time-span. As we worked with a window of 30s during the segmentation process, we had to adjust this to 40s here, in order to get a symmetrical window. Ultimately, a confusion matrix was produced again and accuracy of our method was calculated.

3 Results

##rendering confusion matrix for training data:

conf_mat <- confusion_matrix(targets = data_class$transport_mode_POSMO, predictions = data_class$travel_mode_det)
conf_mat

##rendering confusion matrix for validation data:


val_conf_mat <- confusion_matrix(targets = val_data_class$transport_mode_POSMO, predictions = val_data_class$travel_mode_det)
val_conf_mat

On our own training dataset, we reached a overall accuracy of 91.9% (figure 6). In the predictive validation with the second dataset, overall accuracy was 60.9% and hence considerably lower. Looking at the confusion matrix in figure 7, it can be seen, that the most problematic allocations were between tram and bus, as our method frequently classified bus as tram drives (n=50). The same happened often for cars which tended to be classified as tram (n=21) too. For bus drives our method repeatedly classified car (n=26). On the other hand many walks were treated as bus drives (n=17).

##plot confusion matrix for training data:

plot_confusion_matrix(conf_mat$`Confusion Matrix`[[1]], 
                      add_normalized = FALSE,
                      add_sums = TRUE,
                      add_col_percentages = FALSE,
                      add_row_percentages = FALSE,
                      rm_zero_text = F)

figure 6: confusion matrix for the training data

##plot confusion matrix for validation data:

plot_confusion_matrix(val_conf_mat$`Confusion Matrix`[[1]], 
                      add_normalized = FALSE,
                      add_sums = TRUE,
                      add_col_percentages = FALSE,
                      add_row_percentages = FALSE,
                      rm_zero_text = F)

figure 7: confusion matrix for the validation data

4 Discussion

The accuracy of our travel mode detection procedure in the predictive validation was not as high as those in similar studies, where values above 90% are typically achieved (Sadeghian et al., 2021; Wu et al., 2016). We identified two main problems with our approach.

First, many car or bus journeys were mistakenly classified as tram rides. This may be explained by the fact that tram rails are located on streets also used by cars and buses. Consequently, these trips were allocated to the tram rail as well, resulting in incorrect classification. However, figure 8 shows that the incorrectly classified tram rides had remarkably fewer fixes per segment attributed with tram rails as the nearest feature. Hence, to address this issue, we suggest implementing a threshold. For instance, only segments where 50% of all fixes are closest to a tram rail could be classified as tram rides.

val_data_class <- val_data_class |> #add column with TRUE/FALSE for correct classification
  mutate(test = ifelse(transport_mode_POSMO == travel_mode_det, TRUE, FALSE ))

val_data_class |> #visualizie tram classification by number of fixes attributed
  filter(travel_mode_det == "Tram") |> 
  ggplot(aes(x=test, y=percentage_nearest_route, fill=test))+
  geom_boxplot()+  
  theme_classic()+
  theme(legend.position = "none")+
  labs(x="correct classification", y="percentage of fixes allocated to tram lines")

figure 8: segments of the validation dataset correctly or incorrectly classified as tram by number of fixes attributed to tram rails as nearest feature

Second, the classification of bus rides was poorly. The method of creating a buffer around bus stops thus did not work yet. Alternative approaches could involve using the actual bus network data and similarity measures, but these would require additional computational power again. On the other hand, the current method could be improved by adjusting the buffer size or enhancing the segmentation process. It is crucial to highlight the importance of pre-processing in travel mode detection. Accurate detection relies on segmentation producing trajectories that represent individual movements with the same travel mode and precise defined start and ending points (Wu et al., 2016). This becomes evident in figure 3, were many walking segments have velocities higher than what would be possible in reality. That can be explained by movements transition directly from one to another mode, without static phase, leading to incorrect segmentation (see figure 9)

In addition it has to be mentioned, that the temporal granularity plays a crucial role when calculating moving parameters, as described in Laube & Purves (2011). Conducting sensitivity analyses and simulations would be necessary to determine the appropriate granularity for the purpose here, while we only made first assumptions in this study.

data |> 
  filter(segment_id == "352") |> #filter for a segment with a high max. speed classified as walking
  ggplot(aes(X, Y))+ #plot it by transport_mode
  geom_path()+
  geom_point(aes(color=transport_mode))+
  coord_equal() +
  theme_classic() +
  theme(
    legend.position = c(0.2,0.2),
    legend.title = element_blank()
    )

figure 9: segment 352 (classified as walk) of the training data by manually validated transport mode: the transition from walk to bus did not lead to a segmentation during pre-processing

Regarding our first research question, a rule-based approach for travel mode detection appears applicable for the POSMO data and can be implemented in R as demonstrated in this study. However, the overall accuracy achieved did not meet the state-of-the-art standards. Nevertheless, by incorporating additional criteria, we would expect further improvements.

Regarding the criteria for distinguishing travel modes, as mentioned in research question two, our results heavily relied on the background GIS layers used. In particular, for train and boat segments, utilizing the background data produced convincing results, whereas it did not work well for bus and tram. Here, other criteria need to be added. Overall, our results support the assumption, that contextual information is crucial (Sadeghian et al., 2021). Furthermore, the use of maximal speed and maximal acceleration proved effective in differentiating between bike and car segments. However, some car segments were misclassified as bike segments, while none of the actual bike segments were labeled as cars. Hence, fine-tuning the thresholds by lowering them a bit could improve results. Lastly, maximal and average speed were suitable features for separating bike from walking segments, with only a few misclassifications.

In conclusion, this study established a base for a travel mode detection procedure from mobile GPS data, using both, moving parameters and further contextual information. Even tough, overall accuracy of 60.9% does not yet met the standards, the procedure implemented is a promising starting point for further developments.

5 References

• Genossenschaft Posmo (2022). POSMO Project. Version 22.01.16. https://posmo.coop/
• Gschwend, C. (2015). Relating movement to geographic context: Effects of preprocessing, relation methods and scale [Dissertation, University of Zurich]. 2015, University of Zurich, Mathematisch-naturwissenschaftliche Fakultät. https://doi.org/10.5167/uzh-135565
• Laube, P. (2014). Movement Spaces and Movement Traces. In P. Laube (Hrsg.), Computational Movement Analysis (S. 9–27). Springer International Publishing. https://doi.org/10.1007/978-3-319-10268-9_2
• Laube, P., & Purves, R. S. (2011). How fast is a cow? Cross-Scale Analysis of Movement Data. Transactions in GIS, 15(3), 401–418. https://doi.org/10.1111/j.1467-9671.2011.01256.x
• Li, L., Zhu, J., Zhang, H., Tan, H., Du, B., & Ran, B. (2020). Coupled application of generative adversarial networks and conventional neural networks for travel mode detection using GPS data. Transportation Research Part A: Policy and Practice, 136, 282–292. https://doi.org/10.1016/j.tra.2020.04.005
• Markos, C., & Yu, J. J. Q. (2020). Unsupervised Deep Learning for GPS-Based Transportation Mode Identification. 2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC), 1–6. https://doi.org/10.1109/ITSC45102.2020.9294673
• Nitsche, P., Widhalm, P., Breuss, S., Brändle, N., & Maurer, P. (2014). Supporting large-scale travel surveys with smartphones – A practical approach. Transportation Research Part C: Emerging Technologies, 43, 212–221. https://doi.org/10.1016/j.trc.2013.11.005
• R Core Team (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
• Rykiel, E. (1996). Testing ecological models: The meaning of validation. Ecological Modelling, 90, 229–244. https://doi.org/10.1016/0304-3800(95)00152-2
• Sadeghian, P., Håkansson, J., & Zhao, X. (2021). Review and evaluation of methods in transport mode detection based on GPS tracking data. Journal of Traffic and Transportation Engineering (English Edition), 8(4), 467–482. https://doi.org/10.1016/j.jtte.2021.04.004
• swisstopo (2023). swissTLM3D. Version 2.1. Bundesamt für Landestopografie. https://www.swisstopo.admin.ch/de/geodata/landscape/tlm3d.html
• Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. 216 pp. [https://ggplot2.tidyverse.org]
• Wu, L., Yang, B., & Jing, P. (2016). Travel Mode Detection Based on GPS Raw Data Collected by Smartphones: A Systematic Review of the Existing Methodologies. Information, 7(4). https://doi.org/10.3390/info7040067
• Yu, J. J. Q. (2021). Travel Mode Identification With GPS Trajectories Using Wavelet Transform and Deep Learning. IEEE Transactions on Intelligent Transportation Systems, 22(2), 1093–1103. https://doi.org/10.1109/TITS.2019.2962741
• Zhang, L., Dalyot, S., Eggert, D., & Sester, M. (2012). Multi-stage approach to travel-mode segmentation and classification of GPS traces. ISPRS - International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XXXVIII-4/W25. https://doi.org/10.5194/isprsarchives-XXXVIII-4-W25-87-2011