Kinetic and potential energy derived from Strava data for bicycle rides

bicycleScience.r

# Bicycle Science
# Analysis of data collected with Strava app 


# Load required libraries
library(xlsx)
library(ggplot2)
library(geosphere)
library(zoo)
library(pracma)
library(Dict)


# Loads data from spreadsheet, computes distances, speeds
#
# Parameters
# ----------
#   bikeFile: File containing data collected by Strava. Some pre-processing 
#             is required to get columns:
#             Time - Date/time string in format: 2022-03-15T19:41:07Z
#             Latitude - Recorded GPS latitude location 
#             Longitude - Recorded GPS longitude location 
#             Hours/minutes/seconds - Time data extracted from date/time strings
#             Running time - Hours/minutes/seconds converted to seconds
#             Corrected time - Running time minus first entry
# 
# Returns
# -------
#   bike: Structure with input fields and calculated fields
#         ptDist - Distance between successive points on the route
#         delta.time - Time between successive points
#         Speed - Speed during each time step
#         Distance - Cumulative distance along route

dataPrep <- function(bikeFile){
  # Read raw data
  bike <- read.xlsx(bikeFile,1)
  # Distances between points
  npts <- length(bike$Latitude)
  pts <- matrix(c(bike$Latitude, bike$Longitude), npts, 2)
  bike$ptDist <- c(0,distGeo(pts[1:npts-1,],pts[2:npts,]))
  # Time between recorded data points
  bike$delta.time <- c(1,bike$Corrected.time[2:npts] - bike$Corrected.time[1:npts-1])
  # Speed
  bike$Speed <- bike$ptDist / bike$delta.time
  # Cumulative distance
  bike$Distance <- cumsum(bike$ptDist)
  # Return bike structure
  return(bike)
}

# Simple rolling average smoothing function for a vector of data
# Data is extended by k-1 means of the first k-1 entries to give output vec the 
# same length as input v.
#
# Parameters
# ----------
#   v: Input vector, n x 1
#   k: Length of averaging window, k < n
# 
# Returns
# -------
#   vec: Smoothed vector

vecSmooth <- function(v,k){
  # Extend v by the mean of the first k values
  mean_v <- mean(v[1:k-1])
  # Replicate the first k-1 averages and attach to the beginning
  v <- c(repmat(mean_v,1,k-1), v)
  # Apply rolling average smoother over k values
  vec <- rollmean(v,k)
  return(vec)
}

# Calculates the kinetic, potential and total energies at each time step
#
# Parameters
# ----------
#   bike: Structure of bike data
#   bike_type: Selecte either "mbike" or "ebike"
#   mass: Structure of rider mass and masses of each bike
# 
# Returns
# -------
#   bike: Bike structure with additional fields KE, PE and Total_energy

calcEnergy <- function(bike,bike_type,mass){
  # Acceleration due to gravity (m/s^2)
  g <- 9.81
  
  # Combined mass of rider and bike
  m <- mass["rider"] + mass[bike_type]
  
  # Kinetic energy
  bike$KE <- 1/2 * m * bike$Speed.SG^2
  
  # Subtract minimum elevation from all elevations
  h <- bike$Elevation - min(bike$Elevation)
  
  # Potential energy
  bike$PE <- m * g * h
  
  # Total energy
  bike$Total_energy <- bike$KE + bike$PE
  
  return(bike)
  
}

# Plots a histogram of bike speeds grouped by meters/second
#
# Parameters
# ----------
#   bike: Structure of bike data
#   bike_type: Selecte either "mbike" or "ebike"
#
# Returns
# -------
#   Histogram of bin counts by speed

histSpeeds <- function(bike,bike_type){
  if( bike_type == "ebike"){ 
    plot_title = "Electric Bike"
    bar_color = "gray"
    bar_fill = "red"
  } else {
    plot_title = "Motorless Bike"
    bar_color = "black"
    bar_fill = "gray"    
  }  

  # Plot histograms of speeds
  ggplot(bike, aes(x = Speed.smooth)) + 
    geom_histogram(color = bar_color, fill = bar_fill,binwidth = 1) + 
    xlab("Speed (m/s)") +
    ylab("Counts") + 
    labs(title = plot_title)  
  
}

# Line plot of bike speeds at each time step
#
# Parameters
# ----------
#   bike: Structure of bike data
#   bike_type: Selecte either "mbike" or "ebike"

plotSpeeds <- function(bike,bike_type){
  if( bike_type == "ebike"){ 
    title = "Electric Bike"
  } else {
    title = "Motorless Bike"
  }
  
  plot(ebike$Corrected.time,ebike$Speed.SG,type = "l", 
       xlab = "Time (sec)", 
       ylab = "Speed (m/s)", 
       main = title) 
  
}

# Line plot of kinetic, potential and total energy at each time step
#
# Parameters
# ----------
#   bike: Structure of bike data
#   bike_type: Selecte either "mbike" or "ebike"

plotEnergy <- function(bike,bike_type){
  if( bike_type == "ebike"){ 
    plot_title = "Electric Bike"
  } else {
    plot_title = "Motorless Bike"
  }
  
  # Data frame of total, kinetic and potential energy
  df <- data.frame(bike$Corrected.time,bike$Total_energy,bike$KE,bike$PE)
  
  # Line colors
  cols <- c("#D43F3A", "#5CB85C", "#46B8DA")
  
  # Plot energies
  ggplot(data = df, aes(x = bike.Corrected.time)) + 
    geom_line(aes(y = bike.Total_energy), color = cols[1], size = 1) + 
    geom_line(aes(y = bike.KE), color = cols[2], size = 1) +
    geom_line(aes(y = bike.PE), color = cols[3], size = 1) + 
    xlab("Time (s)") +
    ylab("Energy (J)") + 
    labs(title = plot_title) + 
    scale_shape_discrete(name = "Energy",
                         labels = c("Total","KE","PE"))
}

# Line plot of bike elevation at each time step
#
# Parameters
# ----------
#   bike: Structure of bike data
#   bike_type: Selecte either "mbike" or "ebike"

plotElevation <- function(bike,bike_type){
  if( bike_type == "ebike"){ 
    plot_title = "Electric Bike"
  } else {
    plot_title = "Motorless Bike"
  }  
  plot(bike$Corrected.time,bike$Elevation, 
       xlab = "Time (secs)", 
       ylab = "Elevation (m)",
       main = plot_title)
}

# Time between any two points on the route
# Convert the GPX data from strava into KML using GPX2KML:
# https://gpx2kml.com/?results
# Load the KML into Google Earth and then create pins at the start and end 
# points. Copy the longituded and latitude data into startPoint and endPoint:
# startPoint <- cbind(<start longitude>, <start latitude>)
# endPoint <- cbind(<end longitude>, <end latitude>)
# Start and end points should be close to the route shown in red by the KML path,
# but do not need to be exact. Points used by the function are the closest route
# points to the chosen start/end points.
#
# Parameters
# ----------
#   bike: Structure of bike data
#   startPoint: Longitude and latitude at starting point
#   endPoint: Longitude and latitude at ending point
#
# Returns
# -------
#   travelTime: Time spent between starting and ending points (seconds)
#
# Example
# -------
# Houston_start <- cbind(-79.029666, 35.921946)
# Houston_end <- cbind(-79.028441, 35.920813)
# timeBtwnPoints(mbike,Houston_start,Houston_end)
# [1] 152

timeBtwnPoints <- function(bike,startPoint,endPoint){
  
  # Route points array
  LLpoints <- cbind(bike$Longitude, bike$Latitude)
  
  # Indices of start and end points
  idx_start <- which.min(distm(LLpoints,startPoint))
  idx_end <- which.min(distm(LLpoints,endPoint))
  
  # Time between points (seconds)
  travelTime <- bike$Corrected.time[idx_end] - bike$Corrected.time[idx_start]
  
  return(travelTime)
  
}


# Read data from spreadsheets
ebike <- dataPrep("ebike_LLH.xlsx")
mbike <- dataPrep("mbike_LLH.xlsx")

# Smooth speeds with simple moving average
ebike$Speed.smooth <- vecSmooth(ebike$Speed,5)
mbike$Speed.smooth <- vecSmooth(mbike$Speed,5)

# Smooth speeds with Savitzky-Golay polynomial filtering
ebike$Speed.SG <- savgol(ebike$Speed,101)
mbike$Speed.SG <- savgol(mbike$Speed,101)

# Masses of rider, motorless and electric bikes
mass <- Dict$new(
  rider = 80, 
  mbike = 15, 
  ebike = 28)

# Calculate potential, kinetic and total energies
ebike <- calcEnergy(ebike,"ebike",mass)
mbike <- calcEnergy(mbike,"mbike",mass)

# Plot energies
plotEnergy(ebike,"ebike")
plotEnergy(mbike,"mbike")