Sample script and data for simulating a bicyclist traversing a bicycle route. Companion slides are here: https://tinyurl.com/squiggly-cosmos2018

bicycle_route_time_energy.R

# This script loads in a data file that contains latitude, longitude, and
# altitude of bicycle route, converts the route to a Cartesian path, simulates a
# bicyclist traversing the using constant power, and computes the total travel
# time and energy cost.
#
# MIT Licensed. Copyright 2018-2019 Jason K. Moore

# step 1: download and load the data
# Bicycle route collected from bikeroll.net and converted from GPX to CSV by
# www.gpsvisualizer.com/convert_input.
# this dataframe has columns "latitude", "longitude", and "altitude..ft."
df <- read.csv('first-bridge-shimanami-kaido-bicycle.csv', sep=";")

# step 2: convert latitude and longitude to cartesian x,y coordinate on surface of earth
# convert the altitude to meters and store in new variable
elevation <- df$altitude..ft. * 0.3048  # m
rad_earth <- 6371000  # m
# convert lat and lon from degrees to radians
lat <- df$latitude * pi / 180
lon <- df$longitude * pi / 180
# use equirectangular approximation
# x is east-west and y is north-south
x <- rad_earth * lon * cos(mean(lat))  # [m]
y <- rad_earth * lat  # [m]
# plot the cartesian path on the surface of the earth
# if you subtract the first values of x and y, the numbers won't be so large
plot(x - x[1], y - y[1], type="l")

# step 3: calculate the distance traveled along the path on the surface of the earth
dist <- cumsum(sqrt(diff(x)^2 + diff(y)^2))
# [-1] removes the first item in the sequence, required because diff() causes
# dist to be 1 short
plot(dist, elevation[-1], type="l")

# step 4: calculate the slope and angle at each point
slope <- diff(elevation[-1]) / diff(dist)
plot(dist[-1], slope, type="l")

angle <- atan(slope)
# plot in degrees
plot(dist[-1], angle * 180 / pi, type="l")

# step 5: create functions that interpolate the angle and elevation wrt distance
interp_angle <- approxfun(dist[-1], angle)
interp_elevation <- approxfun(dist, elevation[-1])

# step 6: define some constants
air.density <- 1.2  # kg/m^3
drag.coeff <- 1.15
frontal.area <- 0.55  # meters^2
mass <- 80  # kg
grav.acc <- 9.81  # m/s^2
roll.coeff <- 0.007
max.power <- 100  # watts

# step 7: create functions to calculate the forces
calc.drag.force <- function(velocity) {
  1 / 2 * air.density * drag.coeff * frontal.area * velocity^2
}

calc.roll.force <- function(angle) {
  roll.coeff * mass * grav.acc * cos(angle)
}

calc.grav.force <- function(angle) {
  mass * grav.acc * sin(angle)
}

calc.prop.force <- function(velocity) {
    max.power / velocity
}

# step 8: create functions that calculate the next position and velocity given
# the current position and velocity
calc.new.dist <- function(current.dist, current.vel, del.time, current.angle) {
  current.dist + current.vel * cos(current.angle) * del.time
}

calc.new.vel <- function(current.dist, current.vel, del.time, current.angle) {
  sum.of.forces <- (calc.prop.force(current.vel) - calc.drag.force(current.vel) -
                    calc.roll.force(current.angle) - calc.grav.force(current.angle))
  current.vel + sum.of.forces / mass * del.time
}

# step 9: simulate the system
current.time <- 0  # seconds
# start from the third point because of diff() truncations
current.dist <- dist[3]  # meter
current.vel <- 5  # meter/second
time.interval <- 0.1  # second

# create empty lists to store the simulation results in
times <- list()
distances <- list()
velocities <- list()
powers <- list()

# initialize a counter for storing values
i <- 1

total_distance <- dist[length(dist)]

# loop until bicyclist reaches the end of the route
while (current.dist < total_distance) {

  # print progress of simulation to the screen
  cat('--------------------------------------------\n')
  cat('Time: ', current.time, ' [s]\n')
  cat('Distance along path: ', current.dist, ' [m]\n')
  cat('Velocity: ', current.vel, '[m/s]\n')

  current.angle <- interp_angle(current.dist)

  # store the things we want to keep track of
  times[[i]] <- current.time
  distances[[i]] <- current.dist
  velocities[[i]] <- current.vel
  powers[[i]] <- current.vel * calc.prop.force(current.vel)

  # update the position and velocity
  current.time <- current.time + time.interval
  current.dist <- calc.new.dist(current.dist, current.vel, time.interval, current.angle)
  current.vel <- calc.new.vel(current.dist, current.vel, time.interval, current.angle)

  # update the counter
  i <- i + 1

}

# convert filled lists to vectors
times <- do.call(rbind, times)
distances <- do.call(rbind, distances)
velocities <- do.call(rbind, velocities)
powers <- do.call(rbind, powers)

# step 10: make a subplot showing elevation, velocity, and propulsion power
par(mfrow=c(3,1))
plot(distances, interp_elevation(distances), type="l")
plot(distances, velocities, type="l")
plot(distances, powers, type="l")

# step 11: calculate the final results (time and energy)
cat('-------\nResults\n-------\n')
cat('The total travel time is:', times[length(times)] / 60, ' [min]\n')
cat('The average speed is:', mean(velocities) * 2.23694, ' [mph]\n')

# This is only valid if power is constant, which it is here.
total.energy <- max.power * times[length(times)]

cat('The total energy consumed is:', total.energy / 1000, '[kJ]\n')
cat('The total energy consumed is:', total.energy * 0.000239006, ' [(kilo)calories]\n')

first-bridge-shimanami-kaido-bicycle.csv

type;time;latitude;longitude;altitude (ft);name;desc
T;2018-07-26 05:23:02.174;34.359617000;133.193789000;29.5;BikeRoll_Untitled route;
T;2018-07-26 05:24:02.174;34.359638000;133.193838000;26.2;;
T;2018-07-26 05:25:02.174;34.360502000;133.193728000;31.2;;
T;2018-07-26 05:26:02.174;34.360630000;133.193677000;31.8;;
T;2018-07-26 05:27:02.174;34.360717000;133.193572000;32.2;;
T;2018-07-26 05:28:02.174;34.360787000;133.193407000;32.8;;
T;2018-07-26 05:29:02.174;34.360752000;133.192910000;35.8;;
T;2018-07-26 05:30:02.174;34.360773000;133.192822000;36.4;;
T;2018-07-26 05:31:02.174;34.360998000;133.192225000;57.1;;
T;2018-07-26 05:32:02.174;34.361029000;133.192103000;62.0;;
T;2018-07-26 05:33:02.174;34.361018000;133.191950000;66.9;;
T;2018-07-26 05:34:02.174;34.360931000;133.191800000;71.9;;
T;2018-07-26 05:35:02.174;34.360845000;133.191756000;76.8;;
T;2018-07-26 05:36:02.174;34.360773000;133.191743000;78.7;;
T;2018-07-26 05:37:02.174;34.360645000;133.191757000;80.7;;
T;2018-07-26 05:38:02.174;34.360196000;133.192071000;90.2;;
T;2018-07-26 05:39:02.174;34.360008000;133.192097000;93.2;;
T;2018-07-26 05:40:02.174;34.359804000;133.192041000;96.1;;
T;2018-07-26 05:41:02.174;34.359715000;133.191946000;96.8;;
T;2018-07-26 05:42:02.174;34.359606000;133.191752000;96.1;;
T;2018-07-26 05:43:02.174;34.359560000;133.191624000;94.8;;
T;2018-07-26 05:44:02.174;34.359552000;133.191283000;88.9;;
T;2018-07-26 05:45:02.174;34.359598000;133.191151000;86.9;;
T;2018-07-26 05:46:02.174;34.359648000;133.191082000;85.0;;
T;2018-07-26 05:47:02.174;34.359932000;133.190921000;86.9;;
T;2018-07-26 05:48:02.174;34.359989000;133.190873000;87.6;;
T;2018-07-26 05:49:02.174;34.360041000;133.190814000;93.8;;
T;2018-07-26 05:50:02.174;34.360085000;133.190716000;100.1;;
T;2018-07-26 05:51:02.174;34.360086000;133.190344000;118.8;;
T;2018-07-26 05:52:02.174;34.359958000;133.190122000;124.0;;
T;2018-07-26 05:53:02.174;34.359955000;133.189980000;127.6;;
T;2018-07-26 05:54:02.174;34.359987000;133.189894000;131.2;;
T;2018-07-26 05:55:02.174;34.359995000;133.189816000;134.8;;
T;2018-07-26 05:56:02.174;34.359948000;133.189677000;138.5;;
T;2018-07-26 05:57:02.174;34.359874000;133.189612000;141.7;;
T;2018-07-26 05:58:02.174;34.359323000;133.189612000;162.4;;
T;2018-07-26 05:59:02.174;34.359239000;133.189598000;166.0;;
T;2018-07-26 06:00:02.174;34.359092000;133.189500000;166.7;;
T;2018-07-26 06:01:02.174;34.358997000;133.189481000;161.7;;
T;2018-07-26 06:02:02.174;34.358642000;133.189542000;147.0;;
T;2018-07-26 06:03:02.174;34.358482000;133.189645000;139.1;;
T;2018-07-26 06:04:02.174;34.358439000;133.189651000;138.5;;
T;2018-07-26 06:05:02.174;34.358379000;133.189635000;135.8;;
T;2018-07-26 06:06:02.174;34.358340000;133.189593000;132.9;;
T;2018-07-26 06:07:02.174;34.358312000;133.189526000;128.0;;
T;2018-07-26 06:08:02.174;34.358308000;133.189469000;123.0;;
T;2018-07-26 06:09:02.174;34.358332000;133.189389000;118.1;;
T;2018-07-26 06:10:02.174;34.358408000;133.189327000;115.5;;
T;2018-07-26 06:11:02.174;34.358478000;133.189324000;112.5;;
T;2018-07-26 06:12:02.174;34.358588000;133.189404000;109.9;;
T;2018-07-26 06:13:02.174;34.358673000;133.189428000;107.0;;
T;2018-07-26 06:14:02.174;34.358772000;133.189413000;104.3;;
T;2018-07-26 06:15:02.174;34.358958000;133.189337000;98.8;;
T;2018-07-26 06:16:02.174;34.359051000;133.189248000;98.8;;
T;2018-07-26 06:17:02.174;34.359123000;133.189052000;96.1;;
T;2018-07-26 06:18:02.174;34.359179000;133.188972000;96.8;;
T;2018-07-26 06:19:02.174;34.359399000;133.188858000;98.4;;
T;2018-07-26 06:20:02.174;34.359551000;133.188712000;100.1;;
T;2018-07-26 06:21:02.174;34.359639000;133.188690000;101.0;;
T;2018-07-26 06:22:02.174;34.359900000;133.188688000;102.7;;
T;2018-07-26 06:23:02.174;34.359993000;133.188660000;105.6;;
T;2018-07-26 06:24:02.174;34.360041000;133.188606000;112.9;;
T;2018-07-26 06:25:02.174;34.360126000;133.188391000;125.0;;
T;2018-07-26 06:26:02.174;34.360211000;133.188291000;128.6;;
T;2018-07-26 06:27:02.174;34.360637000;133.188031000;146.7;;
T;2018-07-26 06:28:02.174;34.360752000;133.187846000;156.5;;
T;2018-07-26 06:29:02.174;34.360802000;133.187719000;159.4;;
T;2018-07-26 06:30:02.174;34.360817000;133.187627000;162.1;;
T;2018-07-26 06:31:02.174;34.360725000;133.187464000;164.7;;
T;2018-07-26 06:32:02.174;34.360701000;133.187387000;167.0;;
T;2018-07-26 06:33:02.174;34.360688000;133.187276000;169.6;;
T;2018-07-26 06:34:02.174;34.360715000;133.187117000;171.9;;
T;2018-07-26 06:35:02.174;34.360725000;133.186933000;174.5;;
T;2018-07-26 06:36:02.174;34.360668000;133.186540000;172.2;;
T;2018-07-26 06:37:02.174;34.360602000;133.186435000;168.0;;
T;2018-07-26 06:38:02.174;34.360506000;133.186436000;163.4;;
T;2018-07-26 06:39:02.174;34.360084000;133.186550000;151.6;;
T;2018-07-26 06:40:02.174;34.359907000;133.186593000;139.8;;
T;2018-07-26 06:41:02.174;34.359713000;133.186134000;112.2;;
T;2018-07-26 06:42:02.174;34.358001000;133.182073000;9.2;;
T;2018-07-26 06:43:02.174;34.357271000;133.180271000;0.0;;
T;2018-07-26 06:44:02.174;34.355841000;133.176637000;0.0;;
T;2018-07-26 06:45:02.174;34.355295000;133.175342000;18.4;;
T;2018-07-26 06:46:02.174;34.354769000;133.174033000;57.4;;
T;2018-07-26 06:47:02.174;34.354736000;133.173939000;60.0;;
T;2018-07-26 06:48:02.174;34.354653000;133.173757000;62.3;;
T;2018-07-26 06:49:02.174;34.354320000;133.173638000;67.9;;
T;2018-07-26 06:50:02.174;34.354077000;133.173476000;81.4;;
T;2018-07-26 06:51:02.174;34.354053000;133.173408000;86.9;;
T;2018-07-26 06:52:02.174;34.354032000;133.173307000;93.2;;
T;2018-07-26 06:53:02.174;34.354018000;133.173238000;99.1;;
T;2018-07-26 06:54:02.174;34.353999000;133.173174000;105.0;;
T;2018-07-26 06:55:02.174;34.353931000;133.173060000;110.9;;
T;2018-07-26 06:56:02.174;34.353586000;133.172885000;132.9;;
T;2018-07-26 06:57:02.174;34.353513000;133.172807000;138.1;;
T;2018-07-26 06:58:02.174;34.353470000;133.172729000;143.4;;
T;2018-07-26 06:59:02.174;34.353457000;133.172618000;148.6;;
T;2018-07-26 07:00:02.174;34.353467000;133.172484000;150.9;;
T;2018-07-26 07:01:02.174;34.353529000;133.172380000;153.5;;
T;2018-07-26 07:02:02.174;34.353725000;133.172224000;158.5;;
T;2018-07-26 07:03:02.174;34.353916000;133.172101000;161.7;;
T;2018-07-26 07:04:02.174;34.354044000;133.171999000;161.1;;
T;2018-07-26 07:05:02.174;34.354290000;133.171836000;153.9;;
T;2018-07-26 07:06:02.174;34.354408000;133.171841000;151.2;;
T;2018-07-26 07:07:02.174;34.354476000;133.171870000;148.3;;
T;2018-07-26 07:08:02.174;34.354533000;133.171924000;145.7;;
T;2018-07-26 07:09:02.174;34.354671000;133.172137000;140.1;;
T;2018-07-26 07:10:02.174;34.354806000;133.172172000;137.1;;
T;2018-07-26 07:11:02.174;34.354883000;133.172152000;134.5;;
T;2018-07-26 07:12:02.174;34.354952000;133.172099000;131.6;;
T;2018-07-26 07:13:02.174;34.355016000;133.171986000;127.6;;
T;2018-07-26 07:14:02.174;34.355011000;133.171701000;122.4;;
T;2018-07-26 07:15:02.174;34.355046000;133.171620000;119.4;;
T;2018-07-26 07:16:02.174;34.355162000;133.171546000;116.1;;
T;2018-07-26 07:17:02.174;34.355311000;133.171591000;113.2;;
T;2018-07-26 07:18:02.174;34.355412000;133.171594000;109.9;;
T;2018-07-26 07:19:02.174;34.355497000;133.171565000;107.0;;
T;2018-07-26 07:20:02.174;34.355532000;133.171540000;103.7;;
T;2018-07-26 07:21:02.174;34.355672000;133.171324000;97.8;;
T;2018-07-26 07:22:02.174;34.355835000;133.171181000;93.2;;
T;2018-07-26 07:23:02.174;34.355846000;133.171147000;91.2;;
T;2018-07-26 07:24:02.174;34.355845000;133.171039000;89.2;;
T;2018-07-26 07:25:02.174;34.355796000;133.170973000;86.3;;
T;2018-07-26 07:26:02.174;34.355745000;133.170934000;83.3;;
T;2018-07-26 07:27:02.174;34.355667000;133.170949000;80.4;;
T;2018-07-26 07:28:02.174;34.355612000;133.171003000;77.4;;
T;2018-07-26 07:29:02.174;34.355577000;133.171182000;74.5;;
T;2018-07-26 07:30:02.174;34.355516000;133.171318000;72.8;;
T;2018-07-26 07:31:02.174;34.355424000;133.171377000;70.9;;
T;2018-07-26 07:32:02.174;34.355326000;133.171377000;69.2;;
T;2018-07-26 07:33:02.174;34.355190000;133.171332000;69.2;;
T;2018-07-26 07:34:02.174;34.355071000;133.171354000;69.2;;
T;2018-07-26 07:35:02.174;34.355006000;133.171414000;69.2;;
T;2018-07-26 07:36:02.174;34.354924000;133.171527000;69.6;;
T;2018-07-26 07:37:02.174;34.354813000;133.171537000;69.9;;
T;2018-07-26 07:38:02.174;34.354723000;133.171480000;70.2;;
T;2018-07-26 07:39:02.174;34.354660000;133.171413000;70.5;;
T;2018-07-26 07:40:02.174;34.354540000;133.171384000;70.9;;
T;2018-07-26 07:41:02.174;34.354447000;133.171423000;70.5;;
T;2018-07-26 07:42:02.174;34.354219000;133.171598000;68.6;;
T;2018-07-26 07:43:02.174;34.354141000;133.171595000;69.9;;
T;2018-07-26 07:44:02.174;34.354059000;133.171490000;71.2;;
T;2018-07-26 07:45:02.174;34.353818000;133.171664000;75.1;;
T;2018-07-26 07:46:02.174;34.353361000;133.171991000;87.3;;
T;2018-07-26 07:47:02.174;34.353054000;133.172164000;97.1;;
T;2018-07-26 07:48:02.174;34.352815000;133.172258000;104.0;;
T;2018-07-26 07:49:02.174;34.352598000;133.172334000;110.9;;
T;2018-07-26 07:50:02.174;34.352095000;133.172425000;123.4;;
T;2018-07-26 07:51:02.174;34.351889000;133.172454000;126.3;;
T;2018-07-26 07:52:02.174;34.351717000;133.172460000;127.6;;
T;2018-07-26 07:53:02.174;34.351276000;133.172489000;128.3;;
T;2018-07-26 07:54:02.174;34.351116000;133.172531000;125.7;;
T;2018-07-26 07:55:02.174;34.350888000;133.172620000;119.8;;
T;2018-07-26 07:56:02.174;34.350744000;133.172703000;114.5;;
T;2018-07-26 07:57:02.174;34.350034000;133.173218000;68.6;;
T;2018-07-26 07:58:02.174;34.349844000;133.173348000;58.7;;
T;2018-07-26 07:59:02.174;34.349717000;133.173438000;53.8;;
T;2018-07-26 08:00:02.174;34.349687000;133.173465000;48.9;;
T;2018-07-26 08:01:02.174;34.349622000;133.173526000;44.0;;
T;2018-07-26 08:02:02.174;34.349770000;133.173747000;35.1;;
T;2018-07-26 08:03:02.174;34.349933000;133.173987000;29.5;;
T;2018-07-26 08:04:02.174;34.350125000;133.174226000;28.9;;
T;2018-07-26 08:05:02.174;34.350442000;133.174571000;33.5;;
T;2018-07-26 08:06:02.174;34.350536000;133.174622000;34.8;;
T;2018-07-26 08:07:02.174;34.350678000;133.174665000;35.8;;
T;2018-07-26 08:08:02.174;34.350704000;133.174671000;38.4;;