Lift or slope classifier

The goal of this notebook is to separate ski lifts from ski slopes, using a set of features and an external dataset with the ski lifts of the world (OpenSnowMap). This shouldn’t be too difficult a task, but maybe just difficult enough to justify some feature engineering and training a classifier.

Separating ski lifts from slopes is useful, since the activity’s statistics (average speed, heart rate etc.) can be improved by removing the ski lifts from the data.

My secondary goal is trying out

import json
import os
import struct

from base64 import b64encode
from datetime import timedelta
from pprint import pprint

import iso8601
import numpy
import pandas
import untangle
import geopandas

from cachier import cachier
from shapely.geometry import Point, MultiPoint
from sklearn.linear_model import LogisticRegressionCV, LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.base import TransformerMixin, BaseEstimator
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from IPython.core.display import HTML
from shapely.ops import nearest_points
from haversine import haversine, Unit

I’m using the output of a Polar Vantage V sports watch. The output is a GPX file with all the GPS registrations, and a CSV file. This CSV file has the added boolean column Lift, which is hand-labeled data (True for ski-lift; False for snowboarding or schnapps). These files are loaded in the code block below and will be the training set.

For this, I’m using the untangle library, which converts XML into native Python objects.

Parsing the relevant parts from OpenSnowMap is quite heavy, so I serialise it into lifts.jsonl so it doesn’t need to run each time.

[Download lifts.jsonl]

The dataset consists of nodes and ways. A way can have the label “aerialway”, which seem to be the ski lifts. For each way that has this label, I collect all the nodes and save them in JSON lines format. The XML parsing is done using untangle again.

if not os.path.exists('data/lifts.jsonl'):
    opensnowmap = untangle.parse('data/planet_pistes.osm')
    nodes = {node['id']: (node['lat'], node['lon']) for node in opensnowmap.osm.node}
    with open('data/lifts.jsonl', 'w') as w:
        lifts = []
        for way in opensnowmap.osm.way:
                for tag in way.tag:
                    if tag['k'] == 'aerialway':
                        lifts.append([nodes[nd['ref']] for nd in way.nd])
            except AttributeError:
if os.path.exists('data/lifts.jsonl'):
    lifts = []
    lift_points = []
    with open('data/lifts.jsonl', 'r') as f:
        for i, line in enumerate(f.readlines()):
            skilift = json.loads(line)
            for lat, lon in skilift:
                lift_points.append({'lift_id': i,
                                    'lift_point': Point(float(lat), float(lon))})
    lift_points = geopandas.GeoDataFrame(lift_points)

In the code blocks below, I add some crafted features.

def sinuosity_index(window):
    source = list(window)
    window = []
    for latlon in source:
        lat, lon = numpy.frombuffer(bytes(latlon), dtype=numpy.float32)
        window.append((lat, lon))
    last_lat, last_lon = window[0]
    first_lat, first_lon = previous_lat, previous_lon = window.pop()
    distance = 0.
    for lat, lon in window:
        distance += haversine((previous_lat, previous_lon), (lat, lon), unit=Unit.METERS)
        previous_lat, previous_lon = lat, lon
        sinuosity = (haversine((first_lat, first_lon), (last_lat, last_lon), unit=Unit.METERS)) / distance
    except ZeroDivisionError:
        sinuosity = 0.
    return sinuosity
skilifts_multipoint = MultiPoint(lift_points['lift_point'].tolist())

def distance_to_lift(row):
    query, result = nearest_points(Point(float(row['lat']), float(row['lon'])), 
    row['Distance to ski lift (meters)'] = haversine(
        (query.x, query.y), 
        (result.x, result.y), 
    return row
class ReadSnowboardingDataset(TransformerMixin, BaseEstimator):

    def __init__(self, sinuosity_window=10, altitude_window=15):
        self.sinuosity_window = sinuosity_window
        self.altitude_window = altitude_window

    def fit(self, X):
        return self.transform(X)

    def transform(self, X):
        snowboarding_datasets = []
        for snowboarding_filename in X:
            snowboarding = pandas.read_csv(snowboarding_filename, skiprows=2)

            trackpoints = untangle.parse(snowboarding_filename.replace(".csv", ".gpx"))
            start_time = iso8601.parse_date(trackpoints.gpx.metadata.time.cdata)
            trackpoints = pandas.DataFrame(
                    {'lat': trackpoint['lat'],
                     'lon': trackpoint['lon'],
                     'latlon': numpy.frombuffer(
                         bytes(numpy.float32(trackpoint['lat'])) + 
                     'timestamp': iso8601.parse_date(trackpoint.time.cdata)}
                    for trackpoint in trackpoints.gpx.trk.trkseg.trkpt

            snowboarding['timestamp'] = snowboarding['Time'].apply(str)
            del snowboarding['Time']

            # TODO proper formatting of timedelta or timestamp using proper utils
            trackpoints['timestamp'] = trackpoints['timestamp'] - start_time
            trackpoints['timestamp'] = trackpoints['timestamp'].apply(
                lambda x: str(x).split('days ')[1][0:8]

            snowboarding = snowboarding.merge(trackpoints, on='timestamp')
            snowboarding = snowboarding.apply(distance_to_lift, axis=1)

            snowboarding['Sinuosity index'] = snowboarding['latlon'].rolling(self.sinuosity_window).apply(
                sinuosity_index, raw=True

            # TODO alignment with closest ski lift

            snowboarding['Altitude change (m)'] = snowboarding['Altitude (m)'].diff()
            snowboarding['Altitude change smoothed (m)'] = snowboarding['Altitude change (m)']\
            snowboarding['Speed smoothed (km/h)'] = snowboarding['Speed (km/h)']\

            snowboarding['Absolute speed difference between smoothed and current (km/h)'] = \
                snowboarding['Speed (km/h)'] - snowboarding['Speed smoothed (km/h)']

            snowboarding['Absolute speed difference between smoothed and current (km/h)'] = \
                snowboarding['Speed (km/h)'] - snowboarding['Speed smoothed (km/h)']

            snowboarding['Relative speed difference between smoothed and current (km/h)'] = \
                snowboarding['Absolute speed difference between smoothed and current (km/h)'] / \
                snowboarding['Speed (km/h)']

        return pandas.concat(snowboarding_datasets)

Q: What sorcery is happening with the latlon field?

A: pandas currently makes it hard to apply a function on a rolling window, for Series that are non-numeric [1]. The same goes for rolling calculations that need multiple fields [2]. Therefore, I mash two 32 bit floats into a single 64 bit float, and unpack it in the function sinuosity_index.

class SplitFeaturesClass(TransformerMixin, BaseEstimator):

    def fit(self, X):
        return self.transform(X)

    def transform(self, X):
        snowboarding_selection = X[
            ['Altitude change smoothed (m)',
            'Speed (km/h)',
            'Absolute speed difference between smoothed and current (km/h)',
            'Relative speed difference between smoothed and current (km/h)',
            'Distance to ski lift (meters)',
            'Sinuosity index',

        movement_features = snowboarding_selection[
            ['Altitude change smoothed (m)',
            'Speed (km/h)',
            'Absolute speed difference between smoothed and current (km/h)',
            'Relative speed difference between smoothed and current (km/h)',
            'Distance to ski lift (meters)',
            'Sinuosity index']]\
            .replace([numpy.inf, -numpy.inf], 0.)
        is_lift = snowboarding_selection['Lift']

        return movement_features, is_lift

For this project, I wanted the complexity to be in the feature engineering step, and then just fit a very simple model (logistic regression). In the code block below, the data is split in the features and the target classes (X and y respectively in scikit learn terms).

snowboarding_pipeline = Pipeline([
    ('read_snowboarding_dataset', ReadSnowboardingDataset()),
    ('split_features_class', SplitFeaturesClass()),
snowboarding_filenames = ['./data/Bart_Broere_2020-02-05_14-56-24.csv']
features, is_lift = snowboarding_pipeline.transform(snowboarding_filenames)
# features = PolynomialFeatures(degree=2, interaction_only=True).fit_transform(features)

model = LogisticRegressionCV(max_iter=10000), y=is_lift)

cross_validated_scores = cross_val_score(model,

[1.0, 0.993849938499385, 0.9876998769987699, 0.9408866995073891]
for column, weight in zip(features.columns, list(model.coef_[0])):
    print(f"{column}: {weight}")
Altitude change smoothed (m): 3.5423181866480538
Speed (km/h): 0.21813147529809632
Absolute speed difference between smoothed and current (km/h): -0.026170639562450405
Relative speed difference between smoothed and current (km/h): 0.22554038490129336
Distance to ski lift (meters): 0.0046837551993429835
Sinuosity index: 0.21784616373482674

Although the classifier’s performance is quite good already, its robustness could probably be improved by more labeled data. Currently it’s a classifier that assigns the most value to the altitude change (are we going up?). If we add down-hill ski lifts, this could probably be improved. This hopefully causes weights to shift to features like how constant the speed is. This does require more data labeling, which is boring. Searching for better hyperparameters (like the window for curviness) only makes sense if there’s more training data. With the limited set of training data available now, the hyperparameters can’t be searched reliably.

Ski lifts definitely go down sometimes

O, and I like, but I sometimes miss my debugger.