EHN add strava model to compute power (#25)

doc/api.rst

@@ -72,6 +72,26 @@ Power-profile
 
    metrics.aerobic_meta_model
 
+.. _models_ref:
+
+Models
+======
+
+.. automodule:: skcycling.model
+    :no-members:
+    :no-inherited-members:
+
+.. currentmodule:: skcycling
+
+Power
+-----
+
+.. autosummary::
+   :toctree: generated/
+   :template: function.rst
+
+   model.strava_power_model
+
 .. _utils_ref:
 
 Utilities

skcycling/extraction/gradient.py

@@ -40,7 +40,7 @@ def acceleration(activity, periods=5, append=True):
                                'required. Got {} fields.'
                                .format(activity.columns))
 
-    acceleration = activity['speed'].diff(periods=periods)
+    acceleration = activity['speed'].diff(periods=periods) / periods
 
     if append:
         activity['acceleration'] = acceleration

skcycling/model/__init__.py

@@ -0,0 +1,11 @@
+"""
+The :mod:`skcycling.model` module includes algorithms to model cycling data.
+"""
+
+# Authors: Guillaume Lemaitre <g.lemaitre58@gmail.com>
+#          Cedric Lemaitre
+# License: BSD 3 clause
+
+from .power import strava_power_model
+
+__all__ = ['strava_power_model']

skcycling/model/power.py

@@ -0,0 +1,101 @@
+"""Module to estimate power from data."""
+
+# Authors: Guillaume Lemaitre <g.lemaitre58@gmail.com>
+#          Cedric Lemaitre
+# License: BSD 3 clause
+
+from __future__ import division
+
+import numpy as np
+from scipy import constants
+
+from ..extraction import gradient_elevation
+from ..extraction import acceleration
+
+
+def strava_power_model(activity, cyclist_weight, bike_weight=6.8,
+                       coef_roll_res=0.0045, pressure=101325.0,
+                       temperature=15.0, coef_drag=1, surface_rider=0.32,
+                       use_acceleration=False):
+    """Strava model used to estimate power.
+
+    Parameters
+    ----------
+    activity : DataFrame
+        The activity containing the ride information.
+
+    cyclist_weight : float
+        The cyclist weight in kg.
+
+    bike_weight : float, default=6.8
+        The bike weight in kg.
+
+    coef_roll_res : float, default=0.0045
+        Rolling resistance coefficient.
+
+    pressure : float, default=101325.0
+        Pressure in Pascal.
+
+    temperature : float, default=15.0
+        Temperature in Celsius.
+
+    coef_drag : float, default=1
+        The drag coefficient also known as Cx.
+
+    surface_rider : float, default=0.32
+        Surface area of the rider facing wind also known as S. The unit is m^2.
+
+    use_acceleration : bool, default=False
+        Either to add the power required to accelerate. This estimation can
+        become unstable if the acceleration varies for reason which are not
+        linked to power changes (i.e., braking, bends, etc.)
+
+    Returns
+    -------
+    power : Series
+        The power estimated.
+
+    References
+    ----------
+    .. [1] How Strava Calculates Power
+    https://support.strava.com/hc/en-us/articles/216917107-How-Strava-Calculates-Power
+
+    """
+    if 'gradient-elevation' not in activity.columns:
+        activity = gradient_elevation(activity)
+    if use_acceleration and 'acceleration' not in activity.columns:
+        activity = acceleration(activity)
+
+    temperature_kelvin = constants.convert_temperature(
+        temperature, 'Celsius', 'Kelvin')
+    total_weight = cyclist_weight + bike_weight  # kg
+
+    speed = activity['speed']  # m.s^-1
+    power_roll_res = coef_roll_res * constants.g * total_weight * speed
+
+    # air density at 0 degree Celsius and a standard atmosphere
+    molar_mass_dry_air = 28.97 / 1000  # kg.mol^-1
+    standard_atmosphere = constants.physical_constants[
+        'standard atmosphere'][0]  # Pa
+    zero_celsius_kelvin = constants.convert_temperature(
+        0, 'Celsius', 'Kelvin')  # 273.15 K
+    air_density_ref = (
+        (standard_atmosphere * molar_mass_dry_air) /
+        (constants.gas_constant * zero_celsius_kelvin))  # kg.m^-3
+    air_density = air_density_ref * (
+        (pressure * zero_celsius_kelvin) /
+        (standard_atmosphere * temperature_kelvin))  # kg.m^-3
+    power_wind = 0.5 * air_density * surface_rider * coef_drag * speed**3
+
+    slope = activity['gradient-elevation']  # grade
+    power_gravity = (total_weight * constants.g *
+                     np.sin(np.arctan(slope)) * speed)
+
+    power_total = power_roll_res + power_wind + power_gravity
+
+    if use_acceleration:
+        acc = activity['acceleration']  # m.s^-1
+        power_acceleration = total_weight * acc * speed
+        power_total = power_total + power_acceleration
+
+    return power_total.clip(0)

skcycling/model/tests/test_power.py

@@ -0,0 +1,109 @@
+"""Test the power package which model power using aside data."""
+
+# Authors: Guillaume Lemaitre <g.lemaitre58@gmail.com>
+#          Cedric Lemaitre
+# License: BSD 3 clause
+
+import pytest
+
+import numpy as np
+from numpy.testing import assert_array_less
+from numpy.testing import assert_allclose
+
+import pandas as pd
+from pandas.testing import assert_series_equal
+
+from skcycling.model import strava_power_model
+from skcycling.extraction import gradient_elevation
+from skcycling.extraction import acceleration
+from skcycling.exceptions import MissingDataError
+
+
+speed = np.ones(100) * 5
+distance = np.linspace(0, 99, num=100)
+elevation = np.linspace(0, 9, num=100)
+activity = pd.DataFrame({'speed': speed,
+                         'distance': distance,
+                         'elevation': elevation})
+
+
+@pytest.mark.parametrize(
+    "activity_corrupted, use_acceleration",
+    [(activity.copy().drop(columns='distance'), False),
+     (activity.copy().drop(columns='speed'), True)]
+)
+def test_strava_power_model_error(activity_corrupted, use_acceleration):
+    with pytest.raises(MissingDataError):
+        strava_power_model(activity_corrupted, cyclist_weight=70,
+                           use_acceleration=use_acceleration)
+
+
+def test_strava_power_model_auto_compute():
+    # check that the acceleration and the elevation will be auto-computed
+    power_auto = strava_power_model(activity, cyclist_weight=70)
+
+    activity_ele_acc = activity.copy()
+    activity_ele_acc = gradient_elevation(activity)
+    activity_ele_acc = acceleration(activity_ele_acc)
+    power_ele_acc = strava_power_model(activity_ele_acc, cyclist_weight=70)
+
+    assert_series_equal(power_auto, power_ele_acc)
+
+
+def test_strava_power_model():
+    # at constant speed the acceleration should not have any influence
+    power_without_acc = strava_power_model(activity, cyclist_weight=78,
+                                           use_acceleration=False)
+    power_with_acc = strava_power_model(activity, cyclist_weight=78,
+                                        use_acceleration=True)
+    assert_allclose(power_without_acc, power_with_acc)
+
+    # increase cyclist weight should increase power
+    power_initial = strava_power_model(activity, cyclist_weight=70)
+    power_increase_weight = strava_power_model(activity, cyclist_weight=78)
+    assert_array_less(power_initial, power_increase_weight)
+
+    # increase bike weight should increase power
+    power_initial = strava_power_model(activity, cyclist_weight=70,
+                                       bike_weight=7)
+    power_increase_weight = strava_power_model(activity, cyclist_weight=70,
+                                               bike_weight=8)
+    assert_array_less(power_initial, power_increase_weight)
+
+    # increase the rolling coefficient should increase power
+    power_initial = strava_power_model(activity, cyclist_weight=70,
+                                       coef_roll_res=0.0045)
+    power_increase_cr = strava_power_model(activity, cyclist_weight=70,
+                                           coef_roll_res=0.006)
+    assert_array_less(power_initial, power_increase_cr)
+
+    # increase of the pressure should increase power
+    power_initial = strava_power_model(activity, cyclist_weight=70,
+                                       pressure=101325)
+    power_increase_pressure = strava_power_model(activity, cyclist_weight=70,
+                                                 pressure=110000)
+    assert_array_less(power_initial, power_increase_pressure)
+
+    # decrease the temperature should increase the power
+    power_initial = strava_power_model(activity, cyclist_weight=70,
+                                       temperature=15.0)
+    power_decrease_temperature = strava_power_model(activity,
+                                                    cyclist_weight=70,
+                                                    temperature=10.0)
+    assert_array_less(power_initial, power_decrease_temperature)
+
+    # increase the drag coefficient should increase the power
+    power_initial = strava_power_model(activity, cyclist_weight=70,
+                                       coef_drag=0.5)
+    power_increase_coef_drag = strava_power_model(activity,
+                                                  cyclist_weight=70,
+                                                  coef_drag=1.0)
+    assert_array_less(power_initial, power_increase_coef_drag)
+
+    # increase the rider frontal surface should increase the power
+    power_initial = strava_power_model(activity, cyclist_weight=70,
+                                       surface_rider=0.32)
+    power_increase_surface_rider = strava_power_model(activity,
+                                                      cyclist_weight=70,
+                                                      surface_rider=0.5)
+    assert_array_less(power_initial, power_increase_surface_rider)