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.