Blog

Runkeeper GPS-Tracks mit Python analysieren (Teil 1)

Einen GPS-Track (GPX-Datei) mit Python untersuchen und auf eine interaktive Karte plotten

Fitness-Tracking-Apps wie Runkeeper sammeln eine Menge interessanter Daten. Mit ein paar Python-Kenntnissen können die Läufe noch detaillierter ausgewertet werden als in der App. Im ersten Teil dieser Serie werde ich einen einzelnen GPS-Track (GPX-Datei) mit Python untersuchen, dabei verwende ich Geopandas und Folium. Manche Ideen stammen von diesem Blogbeitrag, aber ich habe einen eleganteren Weg gefunden, ohne durch den GeoDataFrame zu iterieren. Außerdem füge ich nützliche Tooltips zur Karte hinzu. Im zweiten Teil ändere ich den Code, um einen GeoDataFrame mit allen meinen GPS-Tracks zu erhalten. Dieser eignet sich hervorragend zum Plotten und Analysieren der Läufe. Teil 2 enthält ein Jupyter-Notebook, das Sie mit Ihren eigenen GPS-Tracks verwenden können.

Download the GPS Tracks

Um die mit Runkeeper aufgenommenen GPS-Tracks herunterzuladen, muss man sich auf https://runkeeper.com/ einloggen. Dann auf das Zahnrad-Icon klicken, „Account Settings“ wählen und dann „Export Data“. Man bekommt ein ZIP, das die GPS-Tracks als einzelne GPX-Dateien enthält. (Vor ein paar Monaten musste ich 2 Tage warten, bis der endgültige Download-Link erschienen ist. Letzte Woche funktionierte es sofort.)

GPX-Dateien in Geopandas öffnen

Geopandas kann GPX-Dateien öffnen und parsen (dabei wird GDAL verwendet). Um einen GeoDataFrame mit allen Punkten entlang des Tracks zu bekommen:

import pandas as pd
import geopandas as gpd
import numpy as np
from shapely.geometry import LineString

file = "gpx/2021-01-16-110027.gpx"

gdf = gpd.read_file(file, layer='track_points')
gdf.head()


Zu Höhenangaben von Runkeeper siehe Mit Python einem GPS-Track die Höhe hinzufügen. Wir sehen viele Spalten mit „None“-Werten, die wir löschen können. Beachten Sie die Spalten „track_fid“, „track_seg_id“ und „track_seg_point_id“: Sie werden von GDAL verwendet, um verschiedene Tracks, Segmente und Trackpoints zu identifizieren, die sich in einer einzelnen GPX-Datei befinden können. Runkeeper beginnt nach einer Pause ein neues Track-Segment, daher können wir nicht einfach track_seg_point_id als Index verwenden. Ich werde den von Geopandas erstellten Index als „point_id“ benennen und die anderen Spalten löschen. Das Ergebnis enthält die Pausen. Ich wandle die Zeit zu Datetime.

# Better use the geodataframe index as point_id
gdf.index.names = ['point_id'] 
# And turn time to Datetime
gdf['time'] = pd.to_datetime(gdf['time'])

# Drop useless columns
gdf.dropna(axis=1, inplace=True)
gdf.drop(columns=['track_fid', 'track_seg_id', 'track_seg_point_id'], inplace=True)

gdf.head()

Informationen aus dem GPS-Track extrahieren

Wie lange hat der Lauf gedauert (einschließlich Pausen)?

total_time = gdf.iloc[-1]['time'] - gdf.iloc[0]['time']
total_time

Timedelta(‚0 days 00:42:26‘)

Jetzt berechnen wir die Entfernung von Punkt zu Punkt. Zunächst müssen wir die Geometrie auf UTM umstellen, um Koordinaten in Metern statt in Grad zu erhalten. Mit neueren Versionen von geopandas und pyproj ist es einfach, die beste UTM-Zone zu ermitteln:

gdf = gdf.to_crs(gdf.estimate_utm_crs())

Jetzt erstelle ich einen zweiten GeoDataFrame, dessen Index um -1 verschoben ist. Mit diesem Trick ist es einfach, die Entfernung von einem Punkt zum nächsten zu berechnen.

shifted_gdf = gdf.shift(-1)
gdf['time_delta'] = shifted_gdf['time'] - gdf['time']  
gdf['dist_delta'] = gdf.distance(shifted_gdf)

Hinweis: Die letzte Zeile von shifted_gdf besteht aus NaN-Werten. Daher erhalten wir für die letzte Zeile von gdf auch NaN als distance (vom letzten Trackpoint gibt es keinen Abstand zum nächsten Trackpoint).

Nun können wir die Geschwindigkeit in verschiedenen Formaten berechnen:

# speed in various formats
gdf['m_per_s'] = gdf['dist_delta'] / gdf.time_delta.dt.seconds 
gdf['km_per_h'] = gdf['m_per_s'] * 3.6
gdf['min_per_km'] = 60 / (gdf['km_per_h'])

Die zurückgelegte Strecke (in Metern) und die verstrichene Zeit:

gdf['distance'] = gdf['dist_delta'].cumsum()
gdf['time_passed'] = gdf['time_delta'].cumsum()

Weitere nützliche Spalten:

# Minutes instead datetime might be useful
gdf['minutes'] = gdf['time_passed'].dt.seconds / 60

# Splits (in km) might be usefull for grouping
gdf['splits'] = gdf['distance'] // 1000

# ascent is elevation delta, but only positive values
gdf['ele_delta'] = shifted_gdf['ele'] - gdf['ele']  
gdf['ascent'] = gdf['ele_delta']
gdf.loc[gdf.ascent < 0, ['ascent']] = 0

# Slope in %
# (since ele_delta is not really comparable)
gdf['slope'] = 100 * gdf['ele_delta'] / gdf['dist_delta']

# Ele normalized: Startpoint as 0
gdf['ele_normalized'] = gdf['ele'] - gdf.loc[0]['ele']

# slope and min_per_km can be infinite if 0 km/h
# Replace inf with nan for better plotting
gdf.replace(np.inf, np.nan, inplace=True)
gdf.replace(-np.inf, np.nan, inplace=True)

# Note the NaNs in the last line
gdf.tail()

Einfache Statistik:

gdf.describe()

Der (gesamte) Aufstieg in Höhenmetern:

gdf['ascent'].sum()

Wir können die Geometrie wieder zu WGS84 konvertieren. Dies ist besonders dann sinnvoll, wenn wir Tracks aus verschiedenen Regionen vergleichen wollen (Teil 2).

gdf = gdf.to_crs(epsg = 4326)
shifted_gdf = shifted_gdf.to_crs(epsg = 4326)

Jetzt sind Ihre Visualisierungsfähigkeiten gefragt. Z.B. ein einfaches Profil:

sns.relplot(x="distance", y="ele", data=gdf, kind="line")
plt.ylim(0, 50)

Anmerkung: Ich habe für die y-Achse Grenzen gesetzt, da Berlin ziemlich flach ist und ein geringfügiges Auf und Ab stark übertrieben wäre.

Die Punkte mit Linien verbinden

Wir wollen keine Punkte auf der Karte darstellen, sondern die Linien, die diese Punkte verbinden. Ich erstelle einen weiteren GeoDataFrame mit Linien anstelle von Punkten als Geometrie. Beachten Sie, dass dadurch Informationen wie Zeit und Höhe des letzten Punktes wegfallen. Die Lambda-Funktion verwendet LineString (aus dem Modul shapely), um eine Verbindungslinie zu erstellen. Dann verwenden wir die Linie als neue Geometrie des GeoDataFrame.

lines = gdf.iloc[:-1].copy() # Drop the last row
lines['next_point'] =  shifted_gdf['geometry']
lines['line_segment'] = lines.apply(lambda row: LineString([row['geometry'], row['next_point']]), axis=1) 

lines.set_geometry('line_segment', inplace=True, drop=True)
lines.drop(columns='next_point', inplace=True)
lines.index.names = ['segment_id'] 

Wie unterscheiden sich die Splits?

Unsere Spalte "splits" (Abschnitte von 1 km Länge) ist praktisch, um Geschwindigkeit, Steigung usw. entlang der Strecke zu vergleichen. Bin ich am Ende müde geworden? Der folgende Code erstellt eine Tabelle mit vielen interessanten Details für jeden Split, z. B. Mittelwert, Median und Standardabweichung der Geschwindigkeit, Steigung und die für den Split benötigte Zeit.

lines.groupby('splits').aggregate(
     {'km_per_h': ['mean', 'median', 'max', 'std'],
     'min_per_km': ['mean', 'median', 'max', 'std'],
     'ascent': 'sum',
     'time_delta' : 'sum' })

Karte mit Folium

Mit Folium können wir eine interaktive Karte unseres Tracks zeichnen, einschließlich einem praktischen Tooltip mit Informationen zu jedem Linienabschnitt (siehe Screenshot am Anfang dieses Beitrags).

import folium
import branca.colormap as cm
from folium import plugins

# The approx. center of our map
location = lines.dissolve().convex_hull.centroid

Ich möchte Markierungen für jeden Kilometer einzeichnen. Wir können den ersten Punkt eines jeden Splits verwenden.

splitpoints = gdf.groupby('splits').first()

Folium mag keine Zeitstempel, diese Spalten müssen gelöscht werden. Ich entferne auch Zeilen mit der Geschwindigkeit NaN.

lines_notime = lines.drop(columns=['time', 'time_delta', 'time_passed'] )
lines_notime.dropna(inplace=True)

Werte runden, für hübschere Tooltips:

lines_notime['distance'] = (lines_notime['distance']/1000).round(decimals=3)
lines_notime['km_per_h'] = lines_notime['km_per_h'].round(decimals=1)
lines_notime['minutes'] = lines_notime['minutes'].round(decimals=2)
lines_notime['slope'] = lines_notime['slope'].round(decimals=2)

Und schließlich können wir die Karte zeichnen, wobei die Geschwindigkeit als Farbe dargestellt wird. Die Style-Funktion weist jedem Liniensegment Farbe und Tooltip zu. Für die Marker (die "Meilensteine" der Splits) verwende ich folium.plugins.BeautifyIcon, dieses Plugin kann Zahlen auf den Marker plotten.

m = folium.Map(location=[location.y, location.x], zoom_start=15, tiles='cartodbpositron')

# Plot the track, color is speed
max_speed = lines['km_per_h'].max()
linear = cm.LinearColormap(['white', 'yellow', 'red'], vmin=0, vmax=max_speed)
route = folium.GeoJson(
     lines_notime,                         
     style_function = lambda feature: {
            'color': linear(feature['properties']['km_per_h']),
            'weight': 3},
            tooltip=folium.features.GeoJsonTooltip(
            fields=['distance', 'ele', 'km_per_h', 'minutes', 'slope'], 
            aliases=['Distance', 'Elevation', 'Speed', 'Minutes', 'Slope'])
     ) 

m.add_child(linear)
m.add_child(route) 

# Add markers of the splitpoints
for i in range(0,len(splitpoints)):
    folium.Marker(
        location=[splitpoints.iloc[i].geometry.y, splitpoints.iloc[i].geometry.x],
       # popup=str(i),
        icon=plugins.BeautifyIcon(icon="arrow-down",
                                  icon_shape="marker", 
                                  number=i, 
                                  border_color= 'white', 
                                  background_color='lightgreen',
                                  border_width=1)
    ).add_to(m)

m

... und wir bekommen die ganz oben gezeigte Karte.

Zur nächsten Stufe

Es ist schön, einen einzigen GPS-Track auszuwerten. Aber es wäre noch schöner, die Daten aller GPS-Tracks zu analysieren. Ich werde dies im zweiten Teil tun.

JupyterLab Spickzettel

Markdown Syntax

# Überschrift 1
## Überschrift 2 etc.

Hyperlink: [Link Text](Link URL)

**fette Schrift**
*kursive Schrift*
> blockquote

Geordnete Liste:
1. Bla
2. Bla

Ungeordnete Liste:
– Bla
– Bla

`Code steht in Akzenten`

Bild: ![alt text](image.jpg)

Siehe auch: Extended Synthax. Zum Beispiel:

~~Durchgestrichen~~
Tiefgestellt: H~2~O
Hochgestellt: H^2^

Code-Block mit Syntax Highlighting (das sind Akzente):

```python
s = "das ist python code"
print(s)
```

Latex in Markdown

Im Absatz steht die Latex-Formel $5-4\alpha$ in Dollarzeichen.

Freigestellte Formel:

$$f(x) = \sum\limits_1^k x^2 – \frac{1}{x-1}$$

Neues conda environment in JupyterLab anzeigen

Im environment ipykernel installieren:

(my-env)$ conda install ipykernel

Den Kernel registrieren:

(my-env)$ ipython kernel install --user --name=my-env

Im Browser den Tab mit JupyterLab neu laden, dann wird der Kernel angezeigt.

Einen Kernel aus der Liste entfernen:

jupyter kernelspec uninstall my-env

Extensions

Ich verwende folgende Plugins (per conda installieren, nicht in Jupyterlab):

Mit Python einem GPS-Track die Höhe hinzufügen

Wie kann einer GPX-Datei die Höhe hinzugefügt werden? Mit Python ist das dank elevation, geopandas, rasterio und den freien SRTM-Daten relativ leicht.

Leider sind GPS-Geräte ziemlich ungenau, was die Höhe angeht: Der vertikale Fehler ist etwa bei ± 15 m (bei einem guten GPS-Signal und freie Sicht auf den Himmel). Für ein Profil entlang eines GPS-Tracks ist die gemessene Höhe unbrauchbar. Tracking-Apps wie Runkeeper umgehen das Problem, indem sie jeden Punkt an einen Webservice schicken, der die Höhe „von der Karte“ zurückgibt, welche dann von der App in die GPX-Datei eingetragen wird (siehe auch Runkeeper GPS-Tracks mit Python analysieren Teil 1 und Teil 2). Wer ungenaue Höhendaten einer GPX-Datei ersetzen oder die fehlende Höhe ergänzen will, kann das aber auch leicht mit Python erledigen. Basierend auf den frei zugänglichen SRTM-Daten können wir leicht die Höhe von Punkten ermitteln und das Ergebnis als GPX speichern.

Die SRTM-Daten wurden von einem Space Shuttle aus per Radar gemessen. Das Modul elevation lädt die Daten unseres Gebiets aus dem Netz herunter und erstellt aus den einzelnen Kacheln eine einzige Datei, die als Geotiff abgespeichert wird. Das ist einfach ein TIF-Bild mit Georeferenz in den Metadaten. Es hat (wie ein Graustufen-Bild) nur einen Kanal und der Wert jedes Pixels ist die Höhe in Meter. Die Auflösung der Datei beträgt 1 Bogensekunde, also ungefähr 1 Pixel alle 30 m. Das Geotiff mit unserem DEM (digital elevation model) können wir mit rasterio öffnen und dann den Wert einzelner Koordinaten auslesen.

Vorbereitung

import pandas as pd
import geopandas as gpd
import numpy as np
import os
import elevation
import rasterio

track_in = 'gpx/2021-05-22-163014.gpx' 
track_out = 'test.gpx'

dem_file = '2021-05-22-163014.tif' 
demfolder = "dem" # Wo das DEM gespeichert wird

GPX-Track mit Geopandas öffnen

Wir können einfach Geopandas verwenden, um die GPX-Datei zu öffnen (siehe auch hier).

gdf = gpd.read_file(track_in, layer='track_points')

Dieser Befehl verwendet intern GDAL. Mit layer='track_points' bekommen wir alle Trackpoints des Tracks, für Wegpunkte bräuchten wir layer='waypoints'.

GDAL gibt einige nutzlose Spalten mit „None“-Werten zurück, wir können sie entfernen. Und ‚time‘ sollte in datetime sein (wichtig, wenn wir als GPX speichern wollen).

gdf.dropna(axis=1, inplace=True)
gdf['time'] = pd.to_datetime(gdf['time'])

gdf.head()

Zu den Spalten ‚track_fid‘ und ‚track_seg_id‘: Sie werden von GDAL verwendet, um verschiedene Tracks und Tracksegmente zu referenzieren, die eventuell in der GPX-Datei enthalten sind. Diese ID-Werte sind in der GPX-Datei nicht zu finden, wenn wir sie in einem Texteditor öffnen. Es ist jedoch wichtig, diese Spalten zu behalten, um den Track wieder in eine GPX-Datei speichern zu können.

Wir brauchen auch die Grenzen unserer Region. Ich füge etwas Rand hinzu, um das DEM für andere Aufgaben in der Region verwenden zu können.

# Get the bounds of the track
minx, miny, maxx, maxy = gdf.dissolve().bounds.loc[0]

# You might want to add some margin if you want to use the DEM for other tracks/tasks
minx, miny, maxx, maxy = bounds  = minx - .05, miny - .05, maxx + .05, maxy + .05

Das DEM mit elevation herunterladen

# Make sure the dem folder does exist 
# (otherwise the module elevation fails)
if not os.path.exists(demfolder):
    os.mkdir(demfolder)

# And use an absolute path, 
# the module elevation does not work with relative paths.    
dem_file = os.path.join(os.getcwd(), demfolder, dem_file)

# Download DEM 
elevation.clip(bounds=bounds, output=dem_file)

# Clean temporary files
elevation.clean()

Das DEM mit rasterio öffnen und sampeln

Die Datei mit rasterio öffnen und die Daten des ersten (und einzigen) Kanals in ein Numpy-Array laden:

dem_data = rasterio.open(dem_file)
dem_array = dem_data.read(1)

Hinweis: Unser DEM verwendet bereits dasselbe CRS wie unsere GPS-Daten (WGS84 = EPSG 4326), wir brauchen also keine Reprojektion. Wir können das CRS und die Größe des DEMs überprüfen: dem_data.crs, dem_data.width, dem_data.height, dem_data.bounds. Die Anzahl der Kanäle gibt dem_data.count.

Das Wichtigste passiert in der nächsten Zeile: dem_data.index() nimmt die Koordinaten (Längengrad, Breitengrad) und gibt die dazugehörigen Indexwerte unseres Arrays zurück. Wir können diese Indexwerte verwenden, um den Wert aus dem Array zu erhalten. Hinweis: GPX verwendet das Tag <ele> </ele> für Höhe, ‚ele‘ muss der Spaltenname sein. Und wir erhalten ein int, aber es muss ein float sein, um den GeoDataFrame als GPX-Datei zu speichern.

gdf['ele'] = dem_array[dem_data.index(gdf['geometry'].x, gdf['geometry'].y)]

# Note that we get ele values as int. 
# To save as GPX, ele must be as float.
gdf['ele'] = gdf['ele'].astype(float)

gdf.head()

Natürlich sind unsere Daten noch in 1-m-Schritten und ein Profil sieht entsprechend stufig aus (siehe unten). Ich habe mein Ergebnis mit der von Runkeeper angegebenen Höhe verglichen: Es gibt nur geringe Unterschiede. Ich vermute, dass sie auch die SRTM-Daten verwenden, aber interpolieren (und einen leichten Weichzeichner einsetzen, um das DEM zu glätten?). Wir können das auch, siehe unten. Aber erst einmal speichern wir unsere Daten als GPX-Datei.

GeoDataFrame als GPX speichern

Geopandas macht es sehr einfach, unseren GeoDataFrame als GPX-Datei zu speichern:

gdf.to_file(track_out, 'GPX', layer='track_points')

Wichtig: Das funktioniert nicht, wenn die Datei bereits existiert! Sie müssen ggf. zuerst eine bestehende Datei löschen.

Wenn wir mit Wegpunkten gearbeitet haben, müssen wir layer='waypoints' verwenden. (Nicht vergessen: ‚time‘ muss in datetime sein, und um die Punkte als Track zu speichern, benötigen wir ‚track_fid‘ und ‚track_seg_id‘.)

Extra: Ein Blick auf das DEM

Wir können auch einen kurzen Blick auf unser DEM werfen. Verwandeln wir es in eine Karte und zeichnen die Strecke ein:

import matplotlib.pyplot as plt

vmin = dem_array.min()
vmax = dem_array.max()
extent = minx, maxx, miny, maxy

fig, ax = plt.subplots()
cax = plt.imshow(dem_array, extent=extent, 
                  cmap='Spectral_r', vmin=vmin, vmax=vmax)
fig.colorbar(cax, ax=ax)
gdf.plot(ax=ax, markersize=1, zorder=11)

Ein histplot der Höhen:

import seaborn as sns
sns.displot(dem_array.ravel())

Glätten und skalieren des DEM

Um ein glatteres Profil entlang der Strecke zu erhalten, können wir unser DEM glätten und hochskalieren. Hinweis: Die hochskalierten Daten sind nicht unbedingt genauer. Zuerst müssen wir unser Array in Float umwandeln. Dann können wir einen leichten Weichzeichner anwenden, um die Stufen im Terrain zu glätten. Natürlich sollten die gerundeten Werte immer noch mit den ursprünglichen int-Werten übereinstimmen (sigma=0,33 funktioniert in meinem Beispiel gut, der Wert muss ggf. angepasst werden).

dem_array_float = dem_array.astype(float)

from scipy import ndimage

# Blur to reduce the steps in the DEM.
# With sigma 0.33 I am still very close to the original.
dem_array_float = ndimage.gaussian_filter(dem_array_float, sigma=0.33)

Jetzt können wir unser Array hochskalieren, um alle paar Meter ein Pixel zu erhalten. Hinweis: scipy.ndimage.zoom() mit order=3 (Standard) verwendet kubische Interpolation. Wir brauchen auch eine neue Transformationsmatrix. Sie wird von rasterio verwendet, um den Koordinaten (Breitengrad, Längengrad) die jeweiligen Pixelindizes zuzuweisen.

upscale_factor = 5

# Upscale the array
dem_array_float = ndimage.zoom(dem_array_float, upscale_factor).round(decimals=1)


# Scale the transform matrix
transform = dem_data.transform * dem_data.transform.scale(
        (dem_data.width / dem_array_float.shape[-1]),
        (dem_data.height / dem_array_float.shape[-2])
    )

Jetzt können wir die Höhendaten für unsere Trackpunkte ermitteln:

gdf['ele_resample'] = dem_array_float[rasterio.transform.rowcol(transform, gdf['geometry'].x, gdf['geometry'].y)]

Wir können ein Profil zeichnen (nachdem wir die Entfernung entlang der Strecke berechnet haben):

Die Abbildung zeigt einen Teil meines Tracks: blau unter Verwendung der Höhe des ursprünglichen DEM mit int-Werten, orange mit dem geglätteten und hochskalierten DEM. Da ich einen Runkeeper-Track verwendet habe, kann ich die Differenz zwischen meinen Modellen und der von Runkeeper angegebenen Höhe berechnen. Mit describe() bekomme ich:

Mit dem hochskalierten DEM liegt meine Höhe bei mehr als 50 % aller Trackpoints innerhalb von ± 1 m des von Runkeeper angebenen Wertes. Ziemlich beeindruckend! Es gibt einige Ausreißer, aber ich weiß nicht, welches Modell näher am realen Gelände ist.

Folium und Geopandas: FeatureGroup für kategorische Daten

Wie kategorische Daten auf einer mit Folium erstellten Karte einzelnen Ebenen (FeatureGroup) zugeordnet werden können.

Folium ist eine Python-Bibliothek zum Erstellen von interaktiven Karten. Die geplotteten Marker oder Polygone können mit folium.FeatureGroup() einzelnen Ebenen zugeordnet werden, die mit einem Mausklick an- und ausgeschaltet werden können.

Alle Beispiele, die ich im Netz gefunden habe, definieren jede einzelne Ebene explizit mit einem Befehl wie:

 feature_group1 = FeatureGroup(name='Foo')

Wenn ich kategorische Daten in einem Geopandas.GeoDataFrame habe, dann möchte ich sicher nicht per Hand für jede Kategorie eine Ebene erstellen. Zum Glück können wir das automatisieren, indem wir jede initialisierte FeatureGroup einfach an eine Liste anhängen. Hier ein Beispiel mit Vulkanen.

Ein kurzer Blick auf die Daten:

 volcanoes.info()


RangeIndex: 1233 entries, 0 to 1232
Data columns (total 13 columns):
 #   Column             Non-Null Count  Dtype   
---  ------             --------------  -----   
 0   number             1233 non-null   int64   
 1   name               1233 non-null   object  
 2   country            1233 non-null   object  
 3   vtype              1233 non-null   object  
 4   evidence           1233 non-null   object  
 5   last_eruption      1233 non-null   object  
 6   region             1233 non-null   object  
 7   subregion          1233 non-null   object  
 8   elevation          1233 non-null   int64   
 9   rocks              1198 non-null   object  
 10  tectonic_setting   1230 non-null   object  
 11  last_eruption_int  787 non-null    float64 
 12  geometry           1233 non-null   geometry
dtypes: float64(1), geometry(1), int64(2), object(9)
memory usage: 125.4+ KB

Ich möchte nun eine Karte, die nach Gesteinstyp gruppiert ist.

grouped = volcanoes.groupby('rocks')

Da die Farbpaletten von Seaborn besonders schön sind, erzeuge ich damit eine Liste mit Farben im HTML-Hex-Format. Die Anzahl der benötigten Farben ist len(grouped).

import seaborn as sns
pal = sns.color_palette("husl", len(grouped)).as_hex() 

Da ich in diesem Beispiel nicht möchte, dass auch die Basiskarte (das Tilelayer) im LayerControl angezeigt wird, initialisiere ich die Karte ohne tiles und füge ein TileLayer mit control=False hinzu:

m = folium.Map(tiles=None)
folium.TileLayer('cartodbpositron', control=False).add_to(m)

Nun iterieren wir durch den gruppierten GeoDataFrame. Für jede Kategorie wird eine FeatureGroup erzeugt, die an eine Liste angehängt wird. Dann werden Marker für entsprechende Vulkane der zuletzt initialisierten FeatureGroup hinzugefügt.

f_groups = []

for group_name, group_data in grouped:
    f_groups.append(folium.FeatureGroup(group_name))
    color = pal.pop()
    
    for i in range(0,len(group_data)):

        # html for popup of markers
        html=f"""
            <h2>  {group_data.iloc[i]['name']} </h2>
            <small>
            <p> Country: {group_data.iloc[i]['country']}  <br/>
            Elevation: {group_data.iloc[i]['elevation']}  <br/>
            Last Eruption: 
            {group_data.iloc[i]['last_eruption']}  <br/>
            Rocks: {group_data.iloc[i]['rocks']}  <br/>
            Tectonic Setting: 
            {group_data.iloc[i]['tectonic_setting']}  </p>
            </small>
            """

        iframe = folium.IFrame(html=html, width=300, height=200)
        popup = folium.Popup(iframe, max_width=650)

        # Add markers to last FeatureGroup    
        folium.CircleMarker(
            location=[group_data.iloc[i].geometry.y, 
                      group_data.iloc[i].geometry.x],
            radius=5,
            popup=popup,
            tooltip=group_data.iloc[i]['name'],
            fill_color=color,
            stroke = False, 
            fill_opacity = 1,
        ).add_to(f_groups[-1])

    # Add last featureGroup to Map
    f_groups[-1].add_to(m)

folium.LayerControl().add_to(m)

m

Das Alter der ozeanischen Lithosphäre

Die Bewegungen der Plattentektonik werden auf Karten, die das Alter der ozeanischen Kruste zeigen, besonders deutlich.

Alter der ozeanischen Lithosphäre

Die Untersuchung der Ozeanböden brachte entscheidende Argumente, die dazu führten, dass sich die Theorie der Plattentektonik in den 1960er-Jahren durchsetzen konnte. Tatsächlich kann man die Drift der Kontinente und das Werden und Vergehen der Ozeane auf Karten, die das Alter der ozeanischen Kruste zeigen, regelrecht ablesen. Durch die Prozesse der Plattentektonik wird ozeanische Kruste ständig recycelt; sie entsteht neu an den mittelozeanischen Rücken, an den Subduktionszonen verschwindet sie wieder (siehe auch mein Buch Bewegte Bergwelt). Aus diesem Grund ist der Boden der Ozeane im Maßstab der Erdgeschichte nirgends sonderlich alt.

Ich zeige hier zwei Karten, die ich mit PyGMT geplottet habe: Die Weltkarte ist eine Mollweide-Projektion, die Colormap ist „inferno“. Die andere Karte zeigt den Atlantik. Dabei wird besonders deutlich, wie sich dieser Ozean wie eine Schere öffnet, während Sibirien und Alaska noch immer zusammenhängen. Dies ist eine orthographische Projektion.

Alter der ozeanischen Lithosphäre im Atlantik

Seit ich PyGMT entdeckt habe, ist dies eines meiner bevorzugten Werkzeuge, um Karten zu zeichnen. Man kann damit nicht nur die Topografie (DEM) herunterladen und plotten, sondern auch das Alter der ozeanischen Lithosphäre:

age = pygmt.datasets.load_earth_age(resolution="10m")

Mithilfe von Masken habe ich das mit der Topografie kombiniert. Die Lage der Vulkane ist vom Global Volcanism Program.

PyGMT: Verschiedene colormaps für Land und Wasser

Eine korrekte Küstenlinie in der Karte dank grdlandmask, damit Holland nicht unter Wasser ist.

Mit Python und PyGMT ist es ziemlich einfach, eine einfache Topokarte zu zeichnen, die auf den Höhendaten der NASA (Public Domain!) beruht. Allerdings beruht der einfachste Ansatz darauf, dass eine colormap alles unter 0 m NN in Blautönen, alles über 0 m NN in Grün-/Brauntönen einfärbt (z.B. mit cmap=“geo“). Bei einer Weltkarte oder Bergen sieht das Ergebnis gut aus, bei Küstenlinien gibt es aber Probleme. Eine Karte von Holland macht das Problem deutlich:

import pygmt
region=[3.2, 6.8, 51.4, 53.5]
grid = pygmt.datasets.load_earth_relief(resolution="03s", 
       region=region)

fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="M15c", cmap="geo")
fig.colorbar(frame=["a10", "x+lElevation", "y+lm"])
fig.coast(shorelines=True, region=region, resolution="f")
fig.show()

Land unter! In den Tutorials und Docs habe ich keine Anleitung gefunden, wie das Problem gelöst werden kann. Mithilfe von pygmt.grdlandmask war es dann ganz einfach. Der Ansatz beruht auf der Küstenlinie der Funktion fig.coast(). Die erzeugte Maske setzt für jedes Pixel den für „Land“ bzw. „Wasser“ festgelegten Wert, z.B. 0 oder 1. Der erste Trick ist, dass diese Maske einfach mit dem grid multipliziert werden kann, dann werden alle Wasser-Pixel auf 0 gesetzt. Jetzt können wir das Land mit einer geeigneten colormap zeichnen, z.B. cmap=“dem1″. Als spacing für die Maske verwende ich den Wert der resolution des verwendeten Datasets; resolution der Maske ist „f“ (full resolution), wie bei fig.coast().

Das Wasser plotten wir im zweiten Schritt darüber. Allerdings können wir jetzt nicht einfach „maskvalues=[1, 0]“ verwenden, weil dann das Land zu Flachwasser der Tiefe 0 wird. Es ist aber zum Glück möglich, bei der Maske „NaN“ (als string) als Wert anzugeben. Und in fig.grdimage können wir mit nan_transparent=True sicherstellen, dass das Land nun transparent ist.

land = grid * pygmt.grdlandmask(region=region, 
     spacing="03s", 
     maskvalues=[0, 1], 
     resolution="f")
wet = grid * pygmt.grdlandmask(region=region, 
     spacing="03s", 
     maskvalues=[1, "NaN"], 
     resolution="f")

fig = pygmt.Figure()

# Plotte Land (mit Hillshading)
fig.grdimage(grid=land, projection="M15c", 
     cmap="dem1", shading=True) 
fig.colorbar(frame=["a10", "x+lElevation", "y+lm"])

# Plotte Wasser (NaN transparent)
fig.grdimage(grid=wet, 
     cmap="seafloor", nan_transparent=True)

fig.coast(shorelines=True, projection="M15c", 
     region=region, resolution="f")
fig.coast(rivers="1/0.5p,blue") # Flüsse

fig.show()

Trekking-Karte mit PyGMT und Python zeichnen

Einfache Karte mit GPS-Track in Python und PyGMT am Beispiel der Tour du Mont Blanc.

Karte der Tour du Mont Blanc
Karte der Tour du Mont Blanc, erstellt mit PyGMT

Mithilfe von PyGMT ist es ziemlich einfach, gut aussehende Karten zu zeichnen. Als Beispiel zeige ich hier, wie eine einfache Trekking-Karte der Tour du Mont Blanc erstellt werden kann. Im Detail ist die Syntax etwas esoterisch und geht auf die Syntax des General Mapping Tools zurück (siehe Docs von PyGMT).

import pygmt
import geopandas as gpd
import pandas as pd

Wenn der gesamte GPS-Track eine einzige Datei ist, können wir ihn mit dem folgenden Befehl in einen GeoDataFrame einlesen. Handelt es sich um viele Dateien, siehe GPS track (GPX-Datei) in Geopandas öffnen.

track = gpd.read_file('gpx/tmb.gpx', layer='tracks')

Wir können mit track.plot() einen ersten Eindruck des Tracks bekommen, und track.total_bounds gibt die maximalen/minimalen Werte von Länge und Breite.

Mit pygmt.datasets.load_earth_relief lade ich die Topo-Daten (der NASA) aus dem Web. Für diesen Maßstab ist resolution="01s" (1 Bogensekunde) geeignet, bei Übersichtskarten (ganze Regionen, Kontinente, Welt) muss unbedingt eine andere Auflösung gewählt werden. Die Topodaten werden mit grdimage gezeichnet. Die Funktion coast kann nicht nur Küsten, sondern auch große Flüsse und Grenzen zeichnen. Die Funktion plot zeichnet die Geometrie des Tracks.

lat_min, long_min, lat_max, long_max = track.total_bounds

# Add some margin
lat_min -= 0.1
lat_max += 0.1
long_max += 0.06
long_min -= 0.05

region = [lat_min, lat_max, long_min, long_max]

grid = pygmt.datasets.load_earth_relief(resolution="01s", region=region)

fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="M15c", cmap="dem2", shading=True)
#fig.colorbar(frame=["a1000", "x+lElevation", "y+lm"])

# Plot borders
fig.coast(borders=["1/0.5p,red,..-"], resolution="f")

# Plot track
fig.plot(data=track, pen="1.5p,white")
fig.plot(data=track, pen="0.8p,blue")
fig.show()

So sieht die Karte schon ganz gut aus, es fehlen aber noch die Orte und wichtige Gipfel. Ich fange mit einer einfachen Liste an und versuche, die Koordinaten per Geocoding zu bekommen.

peaks_list = [ 'Mont Blanc', 'Grandes Jorasses', 'Mont Dolent']
towns_list = ['Chamonix', 'Les Houches', 'Les Contamines', 'Courmayeur', 'Praz de Fort', 'Champex', 'Martigny']

# Get coordinates with gpd.tools.geocode
# (requires GeoPy)

peaks = gpd.tools.geocode(peaks_list)
peaks['label'] = peaks_list

towns = gpd.tools.geocode(towns_list)
towns['label'] = towns_list

Mit dem Code unten können diese Orte über die obige Karte geplottet werden. Allerdings war Les Contamines auf der Karte nicht zu sehen, weil die von geocode zurückgegebenen Koordinaten nicht korrekt waren. Ich bereite mit einer Liste einen weiteren Geodataframe vor und hänge ihn an. Außerdem definiere ich weitere „places“.

correct_towns = [['Les Contamines', 45.822746, 6.726553]]

correct_towns = pd.DataFrame(correct_towns, columns = ['label', 'latitude', 'longitude'])
correct_towns = gpd.GeoDataFrame(
    correct_towns, geometry=gpd.points_from_xy(correct_towns.longitude, correct_towns.latitude))

towns = towns.append(correct_towns)

# These places did not get usable coordinates
# with geocoding
# and I also want different markers
places = [['Col de la Forclaz', 46.059, 7.0022],
         ['Flégère', 45.9604, 6.8865]]
 
places = pd.DataFrame(places, columns = ['label', 'latitude', 'longitude'])

Jetzt können wir die Karte vervollständigen. Mit plot zeichnen wir verschiedene Markierungen für Städte, Gipfel und andere Orte; text schreibt das Label daneben (x mit kleinem Abstand).

fig.plot(data=towns, style="c0.2c", color="white", pen="black")
fig.text(text=towns.label, x=towns.geometry.x + 0.005, y=towns.geometry.y, justify="ML", font="10p")

fig.plot(x=places.longitude, y=places.latitude, style="c0.15c", color="black")
fig.text(text=places.label, x=places.longitude + 0.005, y=places.latitude, justify="ML", font="8p")

fig.plot(data=peaks, style="t0.2c", color="black")
fig.text(text=peaks.label, x=peaks.geometry.x + 0.005, y=peaks.geometry.y, justify="ML", font="10p")

fig.show()

Leider sind die Topo-Daten nicht immer ganz fehlerfrei, mit virtuellen Gipfeln und Löchern. Und die Küsten von Seen sind bei diesem Maßstab ziemlich ungenau. Die Idee, eine Karte des Torres del Paine Cirquit zu zeichnen, habe ich daher verworfen.

KML mit GoogleEarth Placemarks in Geopandas öffnen

Alle Placemarks in allen Ordnern einer KML-Datei in einen GeoDataFrame mit minidom importieren (Python)

Beim Spielen mit Python und Geopandas wollte ich meine GoogleEarth-Ortsmarken in einen Geopandas GeoDataFrame importieren. Das war komplizierter als erwartet. Geopandas kann theoretisch KML (das von GoogleEarth verwendete Format) öffnen und parsen:

geopandas.read_file() 

Es verwendet die Bibliothek fiona dafür. Diese importiert jedoch nur den ersten Ordner und meine Ortsmarken sind in vielen Ordnern organisiert. Nachdem ich mehrere auf KML spezialisierte Module ausprobiert hatte, stellte sich heraus, dass es einfacher ist, das KML mit dem Standard-Python-Modul minidom zu parsen. KML ist schließlich XML. Ich habe auch eine Spalte mit dem Pfad des Ordners erstellt (die meisten meiner Ortsmarken sind unbenannt, die Ordnernamen sind die nützlichsten Informationen). Das ist nicht wirklich schnell, aber es funktioniert.

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from xml.dom.minidom import *

# KML-Datei öffnen

dom = parse('travel.kml')

# Definiere eine Funktion, um den Pfad einer Ortsmarke zu bekommen

def subfolders(node):
    if node.parentNode == dom.documentElement:
        return ""
    else:
        foldername = node.getElementsByTagName("name")[0].firstChild.data
        path = subfolders(node.parentNode) + "/" + foldername
    return path

# Parse das DOM der KML-Datei
# Für jedes Placemark ein Tupel mit Name, lat, long, Ordnername, Pfad
# Hänge das Tupel an eine Liste von Tupeln an

entries = []
placemarks = dom.getElementsByTagName("Placemark")

for i in placemarks:
    longitude = i.getElementsByTagName("longitude")[0].firstChild.data
    latitude = i.getElementsByTagName("latitude")[0].firstChild.data
    try:
        name = i.getElementsByTagName("name")[0].firstChild.data
    except:
        name = ""
    parent = i.parentNode
    foldername = parent.getElementsByTagName("name")[0].firstChild.data
    path = subfolders(parent) 
    entries.append((name, latitude, longitude, foldername, path)) # Liste von Tupeln


df = pd.DataFrame(entries, columns=('name', 'latitude', 'longitude', 'folder', 'path'))
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude, crs="EPSG:4326"))

Nun können wir mit gdf.head() einen Blick auf unseren GeoDataFrame werfen, ihn als CSV speichern mit gdf.to_csv("travel.csv") und das CSV wieder öffnen mit gdf = gpd.read_file("travel.csv").

Jetzt könnten wir die Placemarks zum Beispiel auf einer einfachen Weltkarte darstellen und auch eine konvexe Hülle („Bounding Box“) um sie herum zeichnen.

# Verwende natural earth dataset als Grundkarte
natworld = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

# Für die convexe Hülle brauchen wir eine Geomtrie mit allen Placemarks
combined = gdf.dissolve()

# Plot
fig, ax = plt.subplots(figsize=(10,5))
natworld.plot(ax=ax, color="darkgrey", edgecolor="lightgrey")
gdf.plot(ax=ax, color="blue", marker=".")
combined.convex_hull.plot(ax=ax, color="none", edgecolor="red")

Ergebnis:

Grafik speichern:

fig.savefig("bounding-box.png")

GPS track (GPX-Datei) in Geopandas öffnen

Es ist einfach, einen kompletten Ordner mit GPS-Tracks (GPX-Dateien) in Geopandas zu öffnen, z.B. zum Plotten von Karten oder um sie als Shapefile zu exportieren.

Wer in Python programmiert, möchte vielleicht GPS-Tracks (aufgenommen mit einem GPS-Gerät oder einem Mobiltelefon und einer Tracking-App) in einem GeoDataFrame nutzen. Es ist einfach, die Tracks aller GPX-Dateien eines bestimmten Ordners zu importieren.

 

import geopandas as gpd
import os
folder = "gpx/"

Ich möchte einen GeoDataFrame mit 2 Spalten: 1) Dateiname und 2) Geometrie des Tracks.

# Leerer GeoDataFrame
track = gpd.GeoDataFrame(columns=['name', 'geometry'], geometry='geometry')

Versuch, alle gpx-Dateien im Ordner zu öffnen. Geopandas verwendet Fiona zum Öffnen/Parsen von gpx. Die Ebene „Tracks“ enthält nur den Track ohne Wegpunkte/Zeitstempel (gut genug, um eine Karte zu zeichnen).

for file in os.listdir(folder):
    if file.endswith(('.gpx')):
        try:
            gdf = gpd.read_file(folder + file, layer='tracks')
            track = track.append(gdf[['name', 'geometry']])
        except:
            print("Error", file)

track.sort_values(by="name", inplace=True)
track.reset_index(inplace=True, drop=True)

Jetzt können wir den GeoDataframe plotten, bearbeiten oder exportieren.

# Als Shapefile speichern
track.to_file(folder + 'track.shp')

# Einfacher Plot
track.plot()

Mit PyGMT können hübsche Karten der GPS-Tracks geplottet werden. Wie GPS-Tracks ausgewertet und auf einer interaktiven Karte gezeichnet werden können, erkläre ich in Runkeeper GPS-Tracks mit Python analysieren (Teil 1) und Runkeeper GPS-Tracks mit Python analysieren (Teil 2).

Vulkane pro Staat

In welchen Ländern gibt es die meisten Vulkane? Ich lerne derzeit die Programmiersprache Python, vor allem um damit Daten zu visualisieren und Karten zu zeichnen. Da mir die Beispiele in den Tutorials zu langweilig waren, habe ich mir interessantere Fragestellungen vorgenommen.

Vulkane der Welt als Heatmap
Vulkane der Welt als Heatmap

Ein interessanter Datensatz ist die Liste der aktiven Vulkane (d. h. mind. 1 Eruption im Holozän) von https://volcano.si.edu/, die man als Excel-Tabelle herunterladen kann.

Als Choropleth-Karte wollte ich die Anzahl der Vulkane pro Staat. Es war dann gar nicht so einfach, den Datensatz mit anderen Datensätzen (den Shapefiles der Länder etc.) zu kombinieren…

Anzahl der Vulkane pro Staat
Choropleth-Karte der Anzahl der Vulkane pro Staat

Dass die USA mit 147 Vulkanen an erster Stelle stehen, hat mich überrascht. Davon sind allerdings 26 % auf den Aleuten und 22 % auf der Alaska-Halbinsel … Auf Platz 2 folgt Indonesien mit 123 Vulkanen, Russland auf Platz 3 mit 99 Vulkanen, dann kommen Japan und Chile.

Da ich den Datensatz nun in Form gebracht habe, kann ich leicht noch Antworten auf andere Fragen finden, die ich mir noch nie gestellt habe. Fun Fact: Tonga hat die meisten Vulkane pro Quadratkilometer! Kein Wunder, dass hier vor allem kleine Inselstaaten auftauchen. Nur ein einziger Nichtinselstaat schafft es in die Top Ten, und zwar El Salvador.

Country Volcanoes per km2
1 Tonga 0.008021
2 Saint Kitts and Nevis 0.007663
3 Dominica 0.006631
4 Grenada 0.002907
5 Saint Vincent and the Grenadines 0.002571
6 Saint Lucia 0.001623
7 Vanuatu 0.001148
8 Comoros 0.000922
9 El Salvador 0.000903
10 Samoa 0.000679
11 Cape Verde 0.000496
12 Iceland 0.000291
13 Japan 0.000283
14 Guatemala 0.000197
15 Costa Rica 0.000196
16 Fiji 0.000164
17 Solomon Islands 0.000141
18 Armenia 0.000134
19 Philippines 0.000133
20 Nicaragua 0.000131
21 Ecuador 0.000122
22 Portugal 0.000119
23 Chile 0.000109
24 Equatorial Guinea 0.000107
25 Djibouti 0.000100