Skip to content

12 Geodaten mit Python

Jump to bottom
Jürgen Hansmann edited this page Jun 15, 2026 · 3 revisions

Geodaten mit Python

Bis jetzt haben wir Koordinaten als gewöhnliche Zahlen behandelt - Ostwerte in einer Variable, Nordwerte in einer anderen, eine for-Schleife drüber, fertig. Das funktioniert für einfache Aufgaben, stösst aber schnell an Grenzen: Was ist mit Polygonen? Mit Koordinatensystemen? Mit der Frage, welcher Punkt in welchem Gebiet liegt?

Mit GeoPandas und Shapely heben wir das auf ein anderes Level - und die Messpunkte, die uns seit Kapitel 4 begleiten, bekommen endlich ihre verdiente Karte.

Kurze Einführung: pandas

GeoPandas baut auf pandas auf - einem Python-Werkzeug für tabellarische Daten. Wir schauen uns pandas in Kapitel 14 genauer an. Hier nur so viel, wie wir für GeoPandas brauchen:

Ein DataFrame ist im Wesentlichen eine Tabelle - Zeilen und Spalten, wie eine CSV-Datei oder ein Excel-Sheet, aber direkt im Python-Code bearbeitbar. pd.read_csv() liest eine CSV-Datei in einen DataFrame ein, auf einzelne Spalten greift man mit df["spaltenname"] zu:

import pandas as pd

df = pd.read_csv("messpunkte.csv")  # lädt die CSV als Tabelle
print(df["name"])                   # gibt die ganze Spalte "name" aus
print(df["hoehe"].max())            # höchster Wert in der Spalte "hoehe"

Das war pandas in 30 Sekunden.

Installation

uv pip install geopandas folium

GeoPandas bringt automatisch Shapely (Geometrieoperationen) und pyproj (Koordinatentransformation) mit.

GeoDataFrame - pandas mit Geometrie

Erinnert ihr euch an messpunkte.csv aus Kapitel 8? Damals haben wir sie Zeile für Zeile mit open() eingelesen, den String mit .split(",") zerlegt und jeden Wert manuell umgewandelt. Das war lehrreich - aber ab jetzt gibt es einen eleganteren Weg.

Ein GeoDataFrame ist ein pandas DataFrame mit einer zusätzlichen Spalte für Geometrien (geometry). Jede Zeile ist ein Messpunkt, eine Linie oder ein Polygon - inklusive Koordinatensystem.

import geopandas as gpd
import pandas as pd

# messpunkte.csv einlesen (kennen wir schon aus Kapitel 8)
df = pd.read_csv("messpunkte.csv")

# DataFrame in GeoDataFrame umwandeln: Ost/Nord → Geometrie
gdf = gpd.GeoDataFrame(
    df,
    geometry=gpd.points_from_xy(df["ost"], df["nord"]),
    crs="EPSG:2056"  # LV95
)

print(gdf)
print(gdf.crs)

crs="EPSG:2056" teilt GeoPandas mit, dass unsere Koordinaten im Schweizer Landeskoordinatensystem LV95 vorliegen. EPSG ist ein weltweites Register von Koordinatensystemen - jedes hat eine eindeutige Nummer. LV95 ist 2056, WGS84 ist 4326. Diese beiden werden uns noch oft begegnen.

Koordinaten umprojizieren

Für die Visualisierung auf Web-Karten brauchen wir WGS84 (Längen-/Breitengrade). Das geht mit einer einzigen Zeile - was in Kapitel 11 noch einen HTTP-Request an einen swisstopo-Dienst brauchte, erledigt pyproj jetzt lokal, ohne Internetverbindung:

gdf_wgs84 = gdf.to_crs("EPSG:4326")

# Jetzt sind die Koordinaten in Länge/Breite
print(gdf_wgs84.geometry[0])  # POINT (7.xxx 46.xxx)

pyproj kennt mehr Koordinatensysteme als die meisten von uns - darunter auch ein paar, die man wirklich nicht braucht. EPSG:2056 und EPSG:4326 reichen für die Schweiz in 99% der Fälle.

Geometrien mit Shapely

Shapely liefert die Geometrieobjekte. Man kann damit Punkte, Linien und Polygone erstellen und geometrische Operationen ausführen.

from shapely.geometry import Point

# Einzelnen Punkt erstellen
p = Point(2600000, 1200000)
print(p.x, p.y)

# Abstand zwischen zwei Punkten (in Metern, da wir in LV95 arbeiten)
p2 = Point(2601500, 1201200)
print(f"Distanz: {p.distance(p2):.1f} m")

# Puffer um einen Punkt (Kreis mit Radius 5 km)
buffer = p.buffer(5000)
print(f"Pufferfläche: {buffer.area / 1e6:.2f} km²")

In Kapitel 7 haben wir eine distanz()-Funktion mit Pythagoras geschrieben - p.distance(p2) macht exakt dasselbe, nur kürzer.

Räumliche Abfragen

Erinnert ihr euch an die Bounding-Box-Übung aus Kapitel 5? Dort haben wir von Hand geprüft, ob ein Punkt innerhalb eines Rechtecks liegt - mit if und and. Ein räumlicher Join macht dasselbe, aber für beliebige Polygone und ganze Datensätze auf einmal:

# Messgebiete (Polygone) laden
zonen = gpd.read_file("messgebiete.geojson")

# Räumlichen Join: welcher Punkt liegt in welcher Zone?
result = gpd.sjoin(gdf, zonen, how="left", predicate="within")
print(result[["name_left", "name_right"]])

sjoin ist das Pendant zu einer SQL-Abfrage mit ST_Within - aber in zwei Zeilen Python, ohne Datenbank.

Auf einer Karte darstellen

import folium

# Zentrum der Karte auf den Mittelpunkt der Punkte setzen
mitte = gdf_wgs84.geometry.union_all().centroid
karte = folium.Map(location=[mitte.y, mitte.x], zoom_start=10)

# Alle Messpunkte als Marker hinzufügen
for _, row in gdf_wgs84.iterrows():
    folium.Marker(
        location=[row.geometry.y, row.geometry.x],
        tooltip=f"{row['name']}: {row['hoehe']} m.ü.M."
    ).add_to(karte)

karte.save("messpunkte.html")

Geodaten speichern

# Als GeoJSON exportieren (Standard für Web-Anwendungen)
gdf.to_file("messpunkte.geojson", driver="GeoJSON")

# Als GeoPackage (moderner, unterstützt mehrere Layer)
gdf.to_file("messpunkte.gpkg", driver="GPKG")

GeoPackage (.gpkg) ist heute der empfohlene Standard für den Austausch von Vektordaten und ersetzt in vielen Workflows das ältere Shapefile-Format. Das Shapefile wurde 1998 eingeführt. Es lebt noch.

Übung: Messpunkte auf der Karte geodaten.py

  1. Lade messpunkte.csv (aus Kapitel 8) und erstelle daraus einen GeoDataFrame in LV95.
  2. Projiziere die Punkte nach WGS84.
  3. Erstelle eine folium-Karte mit den Messpunkten. Zeige im Tooltip Name und Höhe an.
  4. Speichere die Karte als messpunkte.html und öffne sie im Browser.

Übung: Räumlicher Join mit Messgebieten geodaten_zonen.py

  1. Lade zusätzlich messgebiete.geojson als GeoDataFrame.
  2. Führe einen räumlichen Join (gpd.sjoin) durch: welcher Messpunkt liegt in welcher Zone?
  3. Gib aus, welche Punkte keiner Zone zugeordnet werden konnten.
  4. Zeige die Zonen und Punkte zusammen auf einer folium-Karte. Färbe die Zonen unterschiedlich ein (folium.GeoJson).

Tipp: Mit how="left" im sjoin bleiben auch Punkte ohne Zone im Ergebnis (mit NaN in der Zonenspalte).

Weiterführend: echte Kantonsgrenzen

Für Projekte mit realen Schweizer Geodaten bietet swisstopo freie Geobasisdaten an:

  • swissBOUNDARIES3D: Gemeinde-, Bezirks- und Kantonsgrenzen
  • API: data.geo.admin.ch
  • Download: GPKG oder GeoJSON, direkt mit gpd.read_file() ladbar
# Kantonsgrenzen direkt von swisstopo laden (benötigt Internetverbindung)
kantone = gpd.read_file(
    "https://data.geo.admin.ch/ch.swisstopo.swissboundaries3d-kanton-flaeche.fill/"
    "swissboundaries3d-kanton-flaeche.fill/"
    "swissboundaries3d-kanton-flaeche.fill_2056.json"
)
print(kantone[["KANTONSNUM", "NAME"]].head())

In Kapitel 14 lernt ihr, wie ihr solche Datensätze mit pandas noch weiter analysieren und filtern könnt.

Finde mehr Informationen und Ressourcen in Wie weiter?

Clone this wiki locally