Geodaten in Python

Arbeiten mit Geodaten in Python

In diesem Blog-Post erläutern wir, wie in Python effizient mit Geodaten gearbeitet werden kann.

Geodaten können verschiedene Objekte symbolisieren – die wichtigsten sind die folgenden drei.

  • Punkte (z.B. einzelne Adressen oder Messpunkte)
  • Linien (z.B. Strecken)
  • Polygone (z.B. Länder oder PLZ-Gebiete)

Neben Koordinatenformaten können Geodaten auch als Adressen vorliegen.

Wir erklären in diesem Post:

  • wie Adressen in Koordinaten (Breitengrad und Längengrad) umgewandelt werden können (Geocodierung),
  • wie Punkte und Polygone auf einer interaktiven Karte dargestellt werden können (Visualisierung),
  • und wie Distanzen berechnet werden können und eine damit eine K-Nearest-Neighbors-Regression auf Geodaten gemacht werden kann (Distanzberechnung und Regression auf Geodaten)

Geocodierung (Geocoding)

Voraussetzungen

Die benötigten Bibliotheken für diesen Abschnitt lassen sich mit PyPI wie folgt installieren (hier für Python3):

pip3 install geopy
pip3 install ratelimit
pip3 install tqdm # für Fortschrittsbalken

Geocodierung mit GeoPy und Nominatim

Geocodierung bezeichnet die Umwandlung von Adressen in Koordinaten und umgekehrt die Umwandlung von Koordinaten in die dazugehörige Adresse (Reverse Geocoding).

Es gibt eine Reihe an frei verfügbaren Geocoding-APIs, die sich für kleinere Use Cases, also z.B. die Geocodierung von weniger als 10.000 Punkten pro Tag, eignen.

Die APIs manuell aufzurufen ist in den meisten Fällen jedoch gar nicht nötig.geopy ist eine hervorragende Python-Bibliothek für (unter anderem) Geocoding und Reverse Geocoding, die viele APIs unterstützt. Wir nutzen in diesem Beispiel die Nominatim-API, die auf OpenStreetMap (OSM) Daten basiert. Die OSM Daten unterliegen der Open Database License (ODbL).

Eine einzelne Adresse ist in geopy mit minimalem Aufwand geocodiert:

from geopy.geocoders import Nominatim
gc = Nominatim(user_agent="fintu-blog-geocoding-python")
gc.geocode("Unter den Linden 1, Berlin")
# Location(Kommandantenhaus, 1, Unter den Linden, Spandauer Vorstadt, Mitte, Berlin, 10117, Deutschland, (52.51720765, 13.3978343993255, 0.0))

Die Nominatim-API benötigt nicht unbedingt eine Adresse bestehend aus Straße, Hausnummer und Ort, sondern kennt auch viele Geschäftsadressen und Sehenswürdigkeiten. Wir nutzen die API, um eine Liste der 114 größten Fußballstadien in Deutschland mit Koordinaten zu versehen – diese stellen wir im nächsten Abschnitt auf einer Karte dar.

Die Nutzungsbedingungen der Nominatim-API sehen eine maximale Frequenz von einem Aufruf pro Sekunde vor. Wir limitieren daher die API-Aufrufe mit der Bibliothek ratelimit.

import pandas as pd
from ratelimit import limits, sleep_and_retry

# Fortschrittsbalken für Pandas Apply
from tqdm import tqdm
tqdm.pandas()

# Stadien Import
df = pd.read_csv("./data/stadien.csv")

@sleep_and_retry
@limits(1,1)
def rate_limited_geocode(query):
    return gc.geocode(query)

def geocode(row):
    lookup_query = row["Name"] + ", " + row["Stadt"]
    lookup_result = rate_limited_geocode(lookup_query)
    return lookup_result

df["geo_location"] = df.progress_apply(geocode, axis=1)

Die Tabelle enthält nun die vollständigen Adressen und Geokoordinaten der Stadien.

df.loc[0,"geo_location"]
# Location(Signal Iduna Park, 50, Strobelallee, Brünninghausen, Dortmund, Regierungsbezirk Arnsberg, Nordrhein-Westfalen, 44139, Deutschland, (51.49266925, 7.45083008801226, 0.0))

Ein paar der Adressen konnten nicht geocodiert werden – wir ignorieren diese für’s erste.

df[df.geo_location.isnull()].shape[0]
# 9
df = df[~df.geo_location.isnull()]

Zur Geocodierung von mehr als ein paar tausend Adressen am Tag oder für große Batch-Geocodierungen empfiehlt es sich, einen eigenen Nominatim-Server aufzusetzen. Eine Anleitung dazu findet sich hier.

Visualisierung

Voraussetzungen

Wir nutzen in diesem Abschnitt Jupyter Notebooks / Jupyter Lab. Eine Anleitung zur Installation von Jupyter findet sich hier. Die benötigten Bibliotheken lassen sich mit PyPI wie folgt installieren (hier für Python3):

pip install ipyleaflet
jupyter nbextension enable --py --sys-prefix ipyleaflet # can be skipped for notebook 5.3 and above

# Falls nocht nicht geschehen, müssen iPython Widgets aktiviert werden:
jupyter nbextension enable --py --sys-prefix widgetsnbextension

# Für Jupyter Lab muss zusätzlich die Jupyter Lab Extension installiert werden:
jupyter labextension install @jupyter-widgets/jupyterlab-manager
jupyter labextension install jupyter-leaflet

Nun lassen sich in Jupyter / Jupyter Lab interaktive Karten erzeugen:

from ipyleaflet import Map
m = Map(center=(52.204793, 360.121558),zoom=4, scroll_wheel_zoom=True)

Interaktive Karten mit leaflet.js und ipyleaflet

Wir können nun eine Karte der geocodierten Fußballstadien erstellen.

from ipyleaflet import Map, Marker, MarkerCluster
m = Map(center=(50.721283, 9.056725),zoom=5, scroll_wheel_zoom=True)

def create_marker(row):
    loc = row["geo_location"]
    name = row["Name"]
    lat_lon = (loc.latitude, loc.longitude)
    return Marker(location=lat_lon,
                    draggable=False,
                    title=name)

markers = df.apply(create_marker, axis=1)
marker_cluster = MarkerCluster(markers=tuple(markers.values))
m.add_layer(marker_cluster)

m

Zur besseren Übersichtlichkeit ändern wir den Hintergrund auf eine monotonere Karte. Wir ersetzen zudem die Marker durch Kreise, deren Radius die Kapazität (Anzahl Plätze) des Stadiums darstellt.

from ipyleaflet import Map, basemaps, basemap_to_tiles, CircleMarker, LayerGroup

background = basemap_to_tiles(basemaps.OpenStreetMap.DE)
m = Map(layers=(background, ),
        center=(50.721283, 9.056725),
        zoom=5,
        scroll_wheel_zoom=True)

biggest_stadium_size = df["Plätze"].max()

def create_marker(row):
    loc = row["geo_location"]
    name = row["Name"]
    size = int((row["Plätze"]/biggest_stadium_size)*15)
    lat_lon = (loc.latitude, loc.longitude)
    return CircleMarker(location=lat_lon,
                    draggable=False,
                    title=name,
                    radius=size,
                    weight=1)

markers = df.apply(create_marker, axis=1)
layer_group = LayerGroup(layers=tuple(markers.values))
m.add_layer(layer_group)

m

Distanzberechnung und Regression auf Geodaten

Voraussetzungen

In diesem Abschnitt nutzen wir wieder geopy sowie für eine K-Nearest-Neighbors-Regression die scikit-learn Bibliothek.

pip3 install geopy
pip3 install scikit-learn

K-Nearest-Neighbors Regression auf Geodaten

In diesem Beispiel verwenden einen Datensatz mit 74.219 geocodierten Messungen von Wasserpumpen in Tansania (Quelle: DrivenData). Die (teilweise fehlenden oder falschen) Geocodierungen wurden dazu schon bereinigt und ergänzt.

df = pd.read_csv("./data/tanzania.csv", low_memory=False)
df.shape[0]
# 74219
df.latitude.isnull().sum() + df.longitude.isnull().sum()
# 0

Die Spalte gps_height, die die Höhe des Messpunkts angibt, enthält allerdings noch fehlende Werte für ca. 1/3 aller Messpunkte.

df.gps_height.isnull().sum()
# 25637

Wir nutzen eine gewichtete K-Nearest-Neighbors Regression auf den Spalten latitude und longitude um die fehlenden Werte zu schätzen. Zu beachten ist, dass die fehlenden Werte für gps_height eine stark regionale Verteilung aufweisen, was das Ergebnis unserer Regression tendenziell verschlechtert.

from ipyleaflet import Map, basemaps, basemap_to_tiles, CircleMarker, LayerGroup

# Wir nutzen aus Performancegründen 1000 zufällig ausgewählte Messungen
df_sel = df.sample(n=1000, random_state=2000)

background = basemap_to_tiles(basemaps.OpenStreetMap.BlackAndWhite)
m = Map(layers=(background, ),
        center=(-6.57729778901325, 35.11030912399293),
        zoom=5,
        scroll_wheel_zoom=True)

def create_marker(row):
    lat_lon = (row["latitude"], row["longitude"])
    is_missing = np.isnan(row["gps_height"])
    if is_missing:
        color = "#FF0000"
    else:
        color = "#0033FF"
    return CircleMarker(location=lat_lon,
                    draggable=False,
                    fill_color=color,
                    fill_opacity=0.5,
                    radius=5,
                    stroke=False)

markers = df_sel.apply(create_marker, axis=1)
layer_group = LayerGroup(layers=tuple(markers.values))
m.add_layer(layer_group)

m

Die roten Kreise markieren Messpunkte mit fehlendem Wert für gps_height, die blauen Kreise markieren Messpunkte, in denen gps_height vorhanden ist.

Die gewichtete K-Nearest-Neighbors Regression weist auf den vorhandenen Daten ein sehr gutes Ergebnis aus.

from geopy.distance import great_circle
from sklearn import neighbors
from sklearn.model_selection import train_test_split

# n_neighbors Regressionsparameter
n_neighbors = 15

# Aufsplittung der Daten nach Vorhandensein von `gps_height`
df_regr = df[~df.gps_height.isnull()][["latitude", "longitude", "gps_height"]]

# Aufsplitten in Training/Test Datensätze
df_train, df_test = train_test_split(df_regr, test_size=0.2, random_state=42)

# Aufsplitten in X und y
def split_x_y(df):
    X = df[["latitude", "longitude"]].as_matrix()
    y = df["gps_height"].as_matrix()
    return X,y
X_train, y_train = split_x_y(df_train)
X_test, y_test = split_x_y(df_test)

# Definition der Regression
regr = neighbors.KNeighborsRegressor(n_neighbors = n_neighbors, weights="distance", n_jobs=-1)

# Training
regr.fit(X_train, y_train)

# Berechnen des R² der Regression
regr.score(X_test, y_test)
# 0.9977175260422764

Der Algorithmus berechnet die Distanz zwischen zwei Punkten als euklidische (L2-) Distanz in den Dimensionen Längengrad und Breitengrad. Wir haben hierbei jedoch noch nicht in Betracht gezogen, dass

  • Längengrad und Breitengrad einen unterschiedlichen Umrechnungsfaktor in tatsächliche Distanz in dieser Dimension (in km) haben
  • Der Umrechnungsfaktor für den Längengrad abhängig vom Breitengrad ist

Dies kann den Algorithmus in der Auswahl des „nächsten Nachbarn“ stören und zu schlechteren Ergebnissen führen. Der Effekt ist natürlich umso größer, je weiter auseinander die Punkte liegen. Auch wenn in unserem Beispiel die Varianz im Breitengradeher gering ist und die Punkte recht nah beieinander liegen, sowie unser Algorithmus mit k=15 sehr robust gegenüber Unschärfen in der Distanz ist, versuchen wir die Regression durch die Berechnung der „echten“ Distanz sowie Gewichtung nach Distanz zu verbessern. Dies ist insbesondere aufgrund der oben erwähnten regionalen Verteilung der Punkte mit fehlendem gps_height-Wert sinnvoll.

from geopy.distance import great_circle
from sklearn import neighbors
from sklearn.model_selection import cross_val_score

# n_neighbors Regressionsparameter
n_neighbors = 15

# Definition der Distanzfunktion
def dist_km(a,b):
    return great_circle(a,b).kilometers

# Aufsplittung der Daten nach Vorhandensein von `gps_height`
df_regr = df[~df.gps_height.isnull()][["latitude", "longitude", "gps_height"]]

# Aufsplitten in X und y
X = df_regr[["latitude", "longitude"]].as_matrix()
y = df_regr["gps_height"].as_matrix()

# Definition der Regression
regr = neighbors.KNeighborsRegressor(n_neighbors = n_neighbors, weights="distance", n_jobs=-1)
regr_realdist = neighbors.KNeighborsRegressor(n_neighbors = n_neighbors, weights="distance", metric=dist_km, n_jobs=-1)

# Berechnung der Genauigkeit der Regression über Cross-Validation (dies kann einige Minuten dauern!)
score = cross_val_score(regr, X, y, cv=5, verbose=2)
score_realdist = cross_val_score(regr_realdist, X, y, cv=5, verbose=2)

df_scores = pd.DataFrame(np.matrix([score, score_realdist]).T, columns=["Länge/Breite Metrik", "KM-Metrik"])
%matplotlib inline
df_scores.boxplot()

Der Test zeigt, dass es für die „dichten“ Punkte mit vorhandenem gps_height-Wert keinen signifikanten Unterschied macht, ob die euklidische Distanz auf Längen- und Breitengrad oder die tatsächliche Distanz in Kilometern verwendet wird.

Da die Punkte, deren gps_height-Wert wir schätzen wollen, jedoch teilweise sehr weit von den im Training des Algorithmus verwendeten Punkten liegen, nutzen wir nichtsdestotrotz die tatsächliche Distanz um die nicht vorhandenen Werte zu schätzen.

# Training des Algorithmus auf dem vollen verfügbaren Datensatz
regr_realdist.fit(X, y)

df_missing = df[df.gps_height.isnull()][["id","latitude","longitude"]]
X_missing = df_missing[["latitude", "longitude"]].as_matrix()
y_missing = regr_realdist.predict(X_missing)
df_missing["gps_height"] = y_missing

df_not_missing = df[~df.gps_height.isnull()]
df_imputation = pd.concat([df_not_missing,df_missing])[["id","gps_height"]]

df_imputed = df.copy()
del df_imputed["gps_height"]
df_imputed = df_imputed.join(df_imputation.set_index("id"), on="id")