Digitalfotos geokodieren mit Google und Python (Software)

Martin Vogel , Dortmund / Bochum, Sun, 19.05.2019, 20:10 (vor 1796 Tagen)

Smartphones sind erstklassige Datensammler. Sie ermöglichen Google und Apple, bemerkenswert genaue Protokolle über die Bewegungen ihrer Kundinnen und Kunden zu führen. Das kann man bedenklich finden und in gewissem Maße auch beschränken, man kann die Daten aber auch für seine eigenen Zwecke nutzen. Google erlaubt es sogar, sämtliche seit Beginn der Nutzung gesammelten Bewegungsdaten herunterzuladen.

https://takeout.google.com/settings/takeout

Das Zusammenstellen dauert ein paar Minuten. Wer die Zeit nicht durch Warten verbringen möchte, kann sich den Link auf die fertige Google-Takeout-Datei auch per Mail schicken lassen. Packt man sie aus, findet man eine Datei mit dem Namen „Standortverlauf.json“, die durchaus mehrere hundert Megabyte groß sein kann. In ihr ist jeder einzelne Datenpunkt als Geokoordinate enthalten, einschließlich einem Wert, der besagt, wie genau die Messung war, manchmal auch um eine Höhenangabe und weitere Daten ergänzt.

Was kann man man nun sinnvolles mit den Daten anfangen? Als jemand, der nicht nur mit dem Smartphone, sondern gelegentlich auch mit einer „richtigen“ Digitalkamera fotografiert, finde ich es immer schade, dass ausgerechnet die besseren Fotos von meiner Bildverwaltung DigiKam nicht mit ihrem korrekten Standort auf einer Karte angezeigt werden. Sofern die Uhr in der Kamera richtig gestellt ist, lässt sich aber über den Standortverlauf für jedes Bild herausfinden, wo denn das Smartphone zu dem Zeitpunkt war. Idealerweise hat man das beim Fotografieren dabei und weiß daher, welche Koordinaten ins Bild einzutragen sind.

Von Hand ist das natürlich etwas mühselig. Allein das Öffnen der riesigen JSON-Datei dauert auf einem mittleren PC eine gefühlte Ewigkeit. Sinnvoller ist es, diese Aufgabe zu automatisieren. DigiKam verfügt dazu über die Funktion „Geo-Korrelator“, die einer Gruppe von Fotos die Geokoordinaten einer Wegaufzeichnung im GPX-Format zuweist. Hier besteht sogar noch die Möglichkeit, Abweichungen durch eine nicht ganz korrekte Uhrzeit der Kamera zu korrigieren.

Jetzt müssen wir nur die JSON-Datei aus dem Google-Takeout ins GPX-Format übertragen und schon können wir tausende von Digitalfotos auf dem Globus verorten (vorausgesetzt, das Handy war dabei, siehe oben). Das unten aufgelistete Python-Programm erledigt genau diese Übersetzung. Anstelle einer Riesen-GPX-Datei werden jedoch handlichere Dateien für einzelne Jahre angelegt. Bei mir waren die GPX-Jahresprotokolle jeweils etwa um die 30 MB groß und recht flott von DigiKams Geokorrelator geladen.

Aus einem Foto-Rundgang durch die mittelalterliche Stadt Lucca wurde so beispielsweise dieser Lageplan:
[image]
Karte: OpenStreetMap, Lizenz: CC-BY-SA 2.0

Die Zahlen geben an, wie viele Fotos an dem jeweiligen Kamerastandort vom Geokorrelator verortet wurden.

Das hier vorgestellte Pythonprogramm übernimmt jedoch nicht alle Koordinaten, sondern nur die mit ausreichender Qualität. GPS-Messungen, deren angegebene Genauigkeit mit einer möglichen Abweichung von mehr als 100 Metern angegeben wird, werden ausgefiltert. Ebenso werden plötzliche Sprünge der Koordinaten um mehrere Kilometer verworfen, die aus ungeklärten Gründen immer wieder in Googles Standortverlauf verzeichnet sind. Um die Sprünge herauszufiltern, betreibe ich allerdings keinen großen Aufwand, sondern verwerfe einfach alle Punkte, wenn seit der letzten gültigen Messung eine Geschwindigkeit von mehr als etwa 150 km/h erforderlich gewesen wäre, um den neuen Punkt zu erreichen. Wer regelmäßig bei höheren Geschwindigkeiten fotografiert, mag das Programm an der Stelle anpassen.

#!/usr/bin/env python3
import datetime
from math import radians, sin, cos, acos
 
print("Konvertiere Standortverlauf.json aus dem Google Takeout zu gpx …")
Startjahr = 2012
Endjahr = 2020
 
def Geschwindigkeit(tim0,lat0,lon0,tim,lat,lon):
    """ Schätzung der Geschwindigkeit zwischen zwei Streckenpunkten """
    # Zeitabstand zwischen den beiden Messungen?
    s = tim - tim0
    # Kein Zeitunterschied?
    if not s:
        return 0
    # Etwas sphärische Geometrie …
    # Herleitung der Formel unten:
    # https://www2.informatik.hu-berlin.de/~sanftleb/Kugelabstand.pdf
    # Winkel in Radiant statt Grad
    lat0, lon0, lat, lon = map(radians, (lat0, lon0, lat, lon))
    # Mittlerer Erdradius: 6371000 Meter
    r = 6371000
    # Das acos-Argument kann durch Rundungsfehler größer als 1 werden und
    # wird daher auf Werte <= 1 begrenzt.
    m = r*acos(min(1,cos(lat0)*cos(lat)*cos(lon0-lon)+sin(lat0)*sin(lat)))
    # Geschwindigkeit in m/s
    return m / s
 
# Eingabedatei aus Google Takeout lesen
with open("Standortverlauf.json") as json:
    tim0 = 0 # Zeitstempel des letzten Koordinatenpunktes
    K = []   # Liste der Koordinatenpunkte
    nk = 0   # Anzahl gefunder Koordinaten
    nzu = 0  # Anzahl verworfener, zu ungenauer Punkte
    nzs = 0  # Anzahl verworfener Punkte (Ausreisser)
    for zeile in json:
        zs = zeile.split(":")
        if "acc" in zs[0]:
            nk += 1
            # acc ist immer der letzte Eintrag zu einem Trackpunkt.
            # lat, lon tms sind da bekannt.
            acc = int(zs[1].split(",")[0])
            # Nur Punkte mit weniger als 100 Metern angegebener Ungenauigkeit
            # verwenden
            if acc < 100:
                if not tim0: # erster Durchlauf
                    tim0, lat0, lon0 = tim, lat, lon
                # nur Punkte verwenden, wenn sie nicht offensichtlich Glitches
                # sind.
                # Kriterium: v < 150 km/h bzw. 42 m/s
                if Geschwindigkeit(tim0,lat0,lon0,tim,lat,lon) < 42:
                    K.append([tim,lat,lon])
                    tim0, lat0, lon0 = tim, lat, lon    
                else:
                    nzs += 1
            else:
                nzu += 1
 
        elif "lat" in zs[0]:
            lat = float(zs[1].split(",")[0])/1e7
        elif "lon" in zs[0]:
            lon = float((zs[1].split(",")[0]))/1e7
        elif "time" in zs[0]:
            tim = float(zs[1].split('"')[1])/1e3
 
print(f"{nk:9} Koordinaten gefunden.")
print(f"{nzu:9} Koordinaten zu ungenau (>100 m).")
print(f"{nzs:9} Koordinaten zu schnell (>42 m/s).")
print(f"{len(K):9} Koordinaten werden umgewandelt.")
 
# Um die Ausgabedatei nicht übermäßig groß werden zu lassen, wird sie
# jahreweise unterteilt.
for jahr in range(Startjahr,Endjahr+1):
    print(jahr)
    js = f"{jahr}-"
    with open(f"Standortverlauf {jahr}.gpx","w") as gpx:
        tim0 = 0
        gpx.write('<?xml version="1.0" encoding="UTF-8" ?>\n')
        gpx.write('<gpx xmlns="http://www.topografix.com/GPX/1/1" version="1.1" '
                  'creator="Standortverlauf.py" '
                  '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">\n')
        gpx.write('<trk>\n')
        for tim, lat, lon in K:
            # ISO-formatierte Zeit
            tif = datetime.datetime.fromtimestamp(tim).isoformat()
            if js in tif:
                # 20 Minuten Pause? Neuer Track.
                if tim > tim0 + 1200: 
                    if tim0:  # jedoch nicht beim ersten Eintrag
                        gpx.write('</trkseg>')
                    gpx.write('<trkseg>\n')
                gpx.write(f'<trkpt lat="{lat}" lon="{lon}">')
                gpx.write(f'<time>{tif}</time></trkpt>\n')
                tim0 = tim
        gpx.write('</trkseg>\n</trk>\n</gpx>\n')