Gdzie się bujasz na rowerze?

Mam znajomego, który lubi jeździć na rowerze za każdym razem inną trasą. Skąd ma wiedzieć czy daną ulicą już jechał czy nie? Pomożemy mu robiąc heatmapę.

Idea

Przede wszystkim – trzeba mieć dane. Dzisiaj będziemy operować na danych zebranych w aplikacji Strava, zatem ruszamy w teren i kręcimy kilometry!

Osobiście trochę jeżdżę, ale nie na tyle żeby inwestować w pełne konto na Stravie. Bo pełne konto ma opcję rysowania heatmapy…

Swoją drogą co miesiąc Strava przygotowuje mapki i publikuje je na swoich stronach i tak dla przykładu jeździ się rowerem po Polsce:

a tak po Warszawie i okolicach:

Każdy chciałby mieć taką swoją mapkę, nie każdy chce tylko za to płacić. Pomożemy!

Eksport danych ze Stravy

Po zalogowaniu się w przeglądarce (w aplikacji na telefonie to nie działa), przechodzimy do swojego profilu

Tam w sekcji “My Account” na dole znajdujemy na dole “Download or Delete Your Account”

Przechodzimy do mechanizmu eksportu danych klikając w “Get Started”

W punkcie drugim – “Download Request” – wciskamy przycisk. Po jakimś czasie dostajemy maila z linkiem do archiwum

Pobieramy je, rozpakowujemy. Interesują nas tylko dwa elementy:

folder activities
plik activities.csv

Oba wypakowujemy do katalogu z naszym skryptem.

Ważne: w katalogu “activities” niektóre pliki są spakowane GZipem (mają rozszerzenie .gz) – warto je wypakować. W Linuxie wystarczy gunzip *.gz w konsoli. W Windows pewnie potrzeba jakiegoś dodatkowego programu (7-Zip dobry na wszystko).

Oczywiście można te pliki wypakować już w Pythonie, ale pominiemy ten fragment.

Przygotowanie danych

Kiedy już zgromadziliśmy nasze dane w odpowiednich strukturach możemy przystąpić do ich przetworzenia.

Zadanie będzie polegało na wyciągnięciu z każdego z plików czasu pomiaru i współrzędnych.

Struktura danych

Wszystkie trasy Strava eksportuje do plików XML w formacie GPX. Z grubsza taki plik wygląda tak:

<?xml version="1.0" encoding="UTF-8"?>
<gpx creator="StravaGPX Android" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd" version="1.1" xmlns="http://www.topografix.com/GPX/1/1">
 <metadata>
  <time>2017-06-12T15:39:43Z</time>
 </metadata>
 <trk>
  <name>Afternoon Ride</name>
  <type>1</type>
  <trkseg>
   <trkpt lat="52.2858350" lon="20.9405770">
    <ele>95.0</ele>
    <time>2017-06-12T15:39:43Z</time>
   </trkpt>
   <trkpt lat="52.2858750" lon="20.9405770">
    <ele>95.0</ele>
    <time>2017-06-12T15:40:35Z</time>
   </trkpt>
  </trkseg>
 </trk>
</gpx>

<?xml version="1.0" encoding="UTF-8"?>

</metadata>

<trk>

<name>Afternoon Ride</name>

</trkpt>

</trkpt>

</trkseg>

</trk>

</gpx>

Mamy sekcję trkseg w której znajdujemy kolekcję (listę) poszczególnych punktów pomiarowych trkpt. Każdy z tych punktów ma swoje współrzędne GPS (i to do heatmapy nam wystarczy), a także inne parametry – tutaj wysokość nad poziomem morza (ele) oraz czas pomiaru (time).

My te dane przetworzymy w zwykłego pandasowego data frame’a. Zatem potrzebujemy przejść po drzewie pliku XML i wyłowić odpowiednie elementy z poszczególnych gałęzi,

Drugą stroną jest plik activities.csv który ma całą masę metadanych o każdej z przebytych tras. Nas będą interesować tylko wybrane informacje:

Activity ID – unikalny ID przebytej trasy
Activity Name – nazwa aktywności
Activity Type – typ aktywności – pozwoli nam na wybranie jedynie tras przebytych na rowerze
Filename – ścieżka do pliku GPX z zapisem trasy

Prawdę mówiąc nie przydadzą nam się do niczego “Activity ID” oraz “Activity Name”, ale zawsze warto je zostawić gdyby np. trzeba było coś zweryfikować. Gotowe (przetworzone ze wszystkich XMLi) dane zapiszemy też do jednego zbiorczego pliku CSV – więc może do dalszej analizy ID i nazwa aktywności się przyda?

Czas na kod

Kod będzie krótki, bo algorytm jest prosty:

każdy z plików XML z zapisem trasy przerabiamy na data frame
dodajemy informacje o typie przejazdu
całość sklejamy w jeden wielki data frame (i sobie go zapisujemy, bo może zechcemy wykorzystać jakieś narzędzie BI do analizy tych danych w innych przekrojach?)
korzystając z NumPy przygotujemy sobie heatmapę do statycznej grafiki
korzystając z Folium przygotujemy interaktywną mapę (i zapiszemy do HTMLa)

Na początek importy:

import glob # lista plików
import xml.etree.cElementTree as ET # przechodzenie po XMLu
import pandas as pd # data frame, zapisywanie danych do CSV
import numpy as np # histogram 2D
from tqdm import tqdm # pasek postępu - żeby było wiadomo że coś się dzieje ;)

import matplotlib.pyplot as plt # wykresy
import matplotlib.cm as cm # palety kolorów

import glob # lista plików

import xml.etree.cElementTree as ET # przechodzenie po XMLu

import pandas as pd # data frame, zapisywanie danych do CSV

import numpy as np # histogram 2D

from tqdm import tqdm # pasek postępu - żeby było wiadomo że coś się dzieje ;)

import matplotlib.pyplot as plt # wykresy

import matplotlib.cm as cm # palety kolorów

Teraz pomocnicza funkcja czytająca jeden plik GPX (XML) i tłumacząca go na pandasowy data frame:

def parse_gpx_file(filename: str) -> pd.DataFrame:
    """Funkcja wczytuje podany plik GPX i zwraca dataframe z jego zawartością"""
    
    def find_tag(xml_node: ET.Element, tag: str) -> ET.Element:
        """Funkcja szuka w XMLu node'a o określonym tagu. Zwraca znaleziony node lub None"""
        found = False
        for node in xml_node:
            if node.tag.endswith(tag):
                found  = True
                break
        
        if found:
            return node

        return None

    # wczytujemy plik XML
    tree = ET.parse(filename)
    root = tree.getroot()

    # począwszy od roota szukamy sekcji "trk"
    trk_node = find_tag(root, "trk")
    if not trk_node:
        return pd.DataFrame()

    # mając sekcję "trk" szukamy w niej "trkseg"
    trkseg_node = find_tag(trk_node, "trkseg")
    if not trkseg_node:
        return pd.DataFrame()

    # w sekcji "trkseg" mamy kolekcję punktów na trasie
    points = []
    for node in trkseg_node:
        # z każdego punktu wyciągamy interesujące nas informacje i dokładamy do listy
        elem = node.attrib
        for nn in node:
            elem[nn.tag.split("}")[-1]] = nn.text
        points.append(elem)
        
    # z listy robimy data frame
    points_df = pd.DataFrame(points)

    # korekcja typów zebranych danych
    points_df['lat'] = points_df['lat'].astype(float)
    points_df['lon'] = points_df['lon'].astype(float)
    points_df['ele'] = points_df['ele'].astype(float)
    points_df['time'] = points_df['time'].apply(lambda s: pd.to_datetime(s))
    points_df['filename'] = filename

    return points_df[['time', 'filename', 'lon', 'lat', 'ele']]

def parse_gpx_file(filename: str) -> pd.DataFrame:

"""Funkcja wczytuje podany plik GPX i zwraca dataframe z jego zawartością"""

def find_tag(xml_node: ET.Element, tag: str) -> ET.Element:

"""Funkcja szuka w XMLu node'a o określonym tagu. Zwraca znaleziony node lub None"""

found = False

for node in xml_node:

if node.tag.endswith(tag):

found = True

break

if found:

return node

return None

# wczytujemy plik XML

tree = ET.parse(filename)

root = tree.getroot()

# począwszy od roota szukamy sekcji "trk"

trk_node = find_tag(root, "trk")

if not trk_node:

return pd.DataFrame()

# mając sekcję "trk" szukamy w niej "trkseg"

trkseg_node = find_tag(trk_node, "trkseg")

if not trkseg_node:

return pd.DataFrame()

# w sekcji "trkseg" mamy kolekcję punktów na trasie

points = []

for node in trkseg_node:

# z każdego punktu wyciągamy interesujące nas informacje i dokładamy do listy

elem = node.attrib

for nn in node:

elem[nn.tag.split("}")[-1]] = nn.text

points.append(elem)

# z listy robimy data frame

points_df = pd.DataFrame(points)

# korekcja typów zebranych danych

points_df['lat'] = points_df['lat'].astype(float)

points_df['lon'] = points_df['lon'].astype(float)

points_df['ele'] = points_df['ele'].astype(float)

points_df['time'] = points_df['time'].apply(lambda s: pd.to_datetime(s))

points_df['filename'] = filename

return points_df[['time', 'filename', 'lon', 'lat', 'ele']]

Do zmiennych przypiszmy ścieżki do poszczególnych zbiorów danych wejściowych i wyjściowych:

activities_data_gpx = "activities/*.gpx"
activities_data = "activities.csv"
activities_data_output = "activities_data.csv"

activities_data_gpx = "activities/*.gpx"

activities_data = "activities.csv"

activities_data_output = "activities_data.csv"

A teraz w jednej pętli możemy przetworzyć wszystkie pliki z trasami:

# pusty dataframe do którego będziemy doklejać dane z kolejnych plików
full_df = pd.DataFrame()

# listujemy pliki i każdy z nich parsujemy
for file in tqdm(glob.glob(activities_data_gpx)):
    tmp_df = parse_gpx_file(file)

    # wynik parsowania doklejamy do pełnego data frame
    full_df = pd.concat([full_df, tmp_df], axis=0)

# pusty dataframe do którego będziemy doklejać dane z kolejnych plików

full_df = pd.DataFrame()

# listujemy pliki i każdy z nich parsujemy

for file in tqdm(glob.glob(activities_data_gpx)):

tmp_df = parse_gpx_file(file)

# wynik parsowania doklejamy do pełnego data frame

full_df = pd.concat([full_df, tmp_df], axis=0)

Wczytajmy metadane tras i połączmy je z już otrzymanymi punktami z poszczególnych przejazdów:

activities = pd.read_csv(activities_data, usecols=['Activity ID', 'Activity Name', 'Activity Type', 'Filename'])

# join dwóch tabel
full_df = pd.merge(full_df, activities, how='left', left_on='filename', right_on='Filename')

# sortujemy po dacie i zapisujemy do CSV "na zaś"
full_df.sort_values('time').to_csv(activities_data_output, index = False)

activities = pd.read_csv(activities_data, usecols=['Activity ID', 'Activity Name', 'Activity Type', 'Filename'])

# join dwóch tabel

full_df = pd.merge(full_df, activities, how='left', left_on='filename', right_on='Filename')

# sortujemy po dacie i zapisujemy do CSV "na zaś"

full_df.sort_values('time').to_csv(activities_data_output, index = False)

Statyczna heatmapa

W zasadzie każdy z punktów na trasie ma swoje współrzędne i są to współrzędne w dwuwymiarowej przestrzeni. A tak się składa, że biblioteka NumPy ma taką funkcję jak histogram2d(). Dzięki temu możemy w prosty sposób przygotować dane do heatmapy, a później – przy pomocy zwykłego matplotlib – je pokazać:

heatmap, xedges, yedges = np.histogram2d(
    # wybieramy tylko aktywności typu "jazda na rowerze"
    x=full_df[(full_df["Activity Type"] == "Ride")]["lon"],
    y=full_df[(full_df["Activity Type"] == "Ride")]["lat"],
    bins=250 # im większa wartość tym większa dokładność
)

extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

fig = plt.figure(figsize=(9,9))
plt.imshow(heatmap.T, extent=extent, origin='lower', cmap=cm.Greys)
plt.show()

heatmap, xedges, yedges = np.histogram2d(

# wybieramy tylko aktywności typu "jazda na rowerze"

x=full_df[(full_df["Activity Type"] == "Ride")]["lon"],

y=full_df[(full_df["Activity Type"] == "Ride")]["lat"],

bins=250 # im większa wartość tym większa dokładność

)

extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]

fig = plt.figure(figsize=(9,9))

plt.imshow(heatmap.T, extent=extent, origin='lower', cmap=cm.Greys)

plt.show()

Moje przejazdy po Warszawie wyglądają mniej więcej tak:

Taki obrazek ma kilka wad. Przede wszystkim nie da się go bezpośrednio nałożyć na mapę – odwzorowanie nie jest takie do jakiego jesteśmy przyzwyczajeni, trzeba by coś porozciągać. Po drugie – bez nałożenia na mapę trudno rozpoznać konkretne miejsca. Ale tak przygotowany obrazek świetnie się sprawdzi jako prerekwizyt do przygotowania infografiki.

Heatmapa interaktywna

Dlaczego jednak nie wykorzystać gotowych rozwiązań i nie przygotować heatmapy nałożonej od razu na mapę? Użyjemy Folium:

import folium
from folium.plugins import HeatMap

1 2	import folium from folium.plugins import HeatMap

Na początek agregujemy dane, bo na tym polega heatmapa. Agregujemy w Pandasie:

# tylko przejazdy rowerem
agg_points = full_df[(full_df["Activity Type"] == "Ride")][['lon', 'lat']].copy()

# zaokrąglenie wartości współrzędnych zmniejsza liczbę punktów (i niestety dokładność)
agg_points['lon'] = agg_points['lon'].apply(lambda x: round(x, 4))
agg_points['lat'] = agg_points['lat'].apply(lambda x: round(x, 4))

agg_points_plot = agg_points.groupby(['lat', 'lon']).size().reset_index()

# tylko przejazdy rowerem

agg_points = full_df[(full_df["Activity Type"] == "Ride")][['lon', 'lat']].copy()

# zaokrąglenie wartości współrzędnych zmniejsza liczbę punktów (i niestety dokładność)

agg_points['lon'] = agg_points['lon'].apply(lambda x: round(x, 4))

agg_points['lat'] = agg_points['lat'].apply(lambda x: round(x, 4))

agg_points_plot = agg_points.groupby(['lat', 'lon']).size().reset_index()

Jedna uwaga odnośnie ostatniej linii powyższego kodu: zliczamy liczbę odczytów położenia w danym punkcie, a nie liczbę tras które przez dany punkt prowadziły. To oznacza, że chodząc w kółko dostaniemy wiele odczytów z grubsza jednego punktu. Lepszym sposobem byłoby policzenie unikalnych Activity ID w danym punkcie – to dałoby unikalną liczbę tras. Jeśli chcemy w ten sposób to robimy nieco inaczej:

agg_points = full_df[(full_df["Activity Type"] == "Ride")][['lon', 'lat', 'filename']].copy()
agg_points['lon'] = agg_points['lon'].apply(lambda x: round(x, 4))
agg_points['lat'] = agg_points['lat'].apply(lambda x: round(x, 4))

agg_points_plot = agg_points.groupby(['lat', 'lon']).nunique().reset_index()

agg_points = full_df[(full_df["Activity Type"] == "Ride")][['lon', 'lat', 'filename']].copy()

agg_points['lon'] = agg_points['lon'].apply(lambda x: round(x, 4))

agg_points['lat'] = agg_points['lat'].apply(lambda x: round(x, 4))

agg_points_plot = agg_points.groupby(['lat', 'lon']).nunique().reset_index()

Mając przygotowane agregaty wyznaczamy środek mapy (tak, żeby nie trzeba było jej jakoś przesuwać za każdym razem), samą mapę (podkład mapowy), a następnie po prostu nakładamy wyliczoną heatmapę:

centrum = (np.mean(agg_points_plot['lat']), np.mean(agg_points_plot['lon']))

# budujemy obiekt Map będący początkowo tylko podkładem z mapą
m = folium.Map(
    location =centrum,
    tiles="OpenStreetMap",
    zoom_start=10,
    zoom_control=True
)

# wyliczamy heatmapę i dodajemy ją do obiektu typu Map
HeatMap(agg_points_plot.values.tolist(), radius=10).add_to(m)

centrum = (np.mean(agg_points_plot['lat']), np.mean(agg_points_plot['lon']))

# budujemy obiekt Map będący początkowo tylko podkładem z mapą

m = folium.Map(

location =centrum,

tiles="OpenStreetMap",

zoom_start=10,

zoom_control=True

)

# wyliczamy heatmapę i dodajemy ją do obiektu typu Map

HeatMap(agg_points_plot.values.tolist(), radius=10).add_to(m)

Na koniec wygenerowaną mapę zapisujemy do pliku HTML:

m.save("heatmap.html")

1	m.save("heatmap.html")

Teraz możemy otworzyć wygenerowany plik na przykład w przeglądarce i cieszyć się wynikiem naszej pracy.

Nie jest tajemnicą gdzie mieszkam, więc powyższa heatmapa to rzeczywiste dane.

Podsumowanie

Co tu się dzisiaj wydarzyło? Nauczyliśmy się aż czterech rzeczy!

jak przeczytać plik XML
jak agregować dwuwymiarowe dane w NumPy
jak agregować je w Pandas
jak przygotować mapkę w Folium

Co można było jeszcze zrobić?

Można policzyć prędkość chwilową pomiędzy kolejnymi punktami trasy i na przykład z tych wartości (po uśrednieniu) wyrysować heatmapę.
Mając czas pomiaru w danym punkcie możemy też pokusić się o sprawdzenie czy po południu osiągamy większą prędkość jazdy niż rano. Dodając dane o kierunku wiatru (to będzie trudne – skąd wziąć tak szczegółowe w lokalizacji i czasie dane?) oraz wyliczając azymut ruchu z dwóch sąsiednich punktów pomiarowych możemy określić czy prędkość jazdy po dojęciu/dodaniu siły wiatru jest stała albo zależna od pory dnia.
W danych mamy również wartość ele czyli wysokość. Można zatem policzyć gradient wysokości w każdym z punktów i znowu zestawić to z prędkością.

W większości przypadków te pomysły są nieco absurdalne, ale zobaczcie co można uzyskać mają serię danych pomiarowych z czterema tylko cechami: czasem pomiaru, długością i szerokością geograficzną oraz wysokością nad poziomem morza. Z tych czterech parametrów wyprowadzamy kolejne: prędkość chwilowa, zmiana wysokości, azymut ruchu, cechy związane z czasem (godzina w ciągu dnia, dzień tygodnia, miesiąc, pora roku, rok, czas jaki upłynął od rozpoczęcia trasy) które możemy ze sobą dowolnie zestawiać.

Nie wspominam tutaj o samych geograficznych cechach i na przykład zastosowaniu algorytmów uczenia maszynowego, chociażby algorytmu DBSCAN który pozwoliłby na znalezienie grup blisko sąsiadujących punktów. Mogłyby to być jakieś ulubione miejsca albo coś tego typu.