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:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
<?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> |
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 iż 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:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 |
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" |
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) |
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) |
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, iż 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() |
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 |
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() |
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, iż 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żeli 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() |
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) |
Na koniec wygenerowaną mapę zapisujemy do pliku 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.