12 Geodaten mit Python

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.

Installation

Die Installation der in diesem Kapitel benötigten Module klappt hoffentlich so:

uv pip install pandas geopandas folium

oder, falls Ihr ohne virtuelle Umgebung arbeiten wollt:

pip install pandas geopandas folium

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

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.

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 - Shapely kennt keine Koordinatensysteme, distance()
# liefert nur eine "unitless" Zahl in der Einheit der Eingabekoordinaten.
# Da p und p2 in LV95 (Meter) vorliegen, ist das Ergebnis hier in Metern.
p2 = Point(2601500, 1201200)
print(f"Distanz: {p.distance(p2):.1f} m")

# Puffer um einen Punkt (Kreis mit Radius 5 km): auch buffer() nimmt die Distanz
# einfach in der Einheit der Koordinaten entgegen - hier 5000 Meter
buffer = p.buffer(5000)

# .area ist ebenfalls unitless, also (Koordinaten-Einheit)² - hier m². Deshalb /1e6 für km²
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?
# how="left" behält alle Punkte aus gdf, auch wenn sie in keiner Zone liegen.
# predicate="within" prüft, ob die Geometrie aus gdf vollständig in der aus zonen liegt.
# Dokumentation: https://geopandas.org/en/stable/docs/reference/api/geopandas.sjoin.html
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

# swisstopo-Pixelkarte als Basemap
karte = folium.Map(
    location=[mitte.y, mitte.x],
    zoom_start=10,
    tiles="https://wmts.geo.admin.ch/1.0.0/ch.swisstopo.pixelkarte-farbe/default/current/3857/{z}/{x}/{y}.jpeg",
    attr="swisstopo",
)

# 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/Katalog: STAC-Katalog von swissBOUNDARIES3D listet alle Versionen und deren Download-Links auf
• Download: nur als gezipptes GeoPackage, Shapefile, File-GDB oder XTF - kein einzelnes GeoJSON mehr

import geopandas as gpd

# Kantonsgrenzen direkt von swisstopo laden (benötigt Internetverbindung).
# /vsizip/vsicurl/ lässt GDAL die ZIP-Datei direkt über HTTP lesen,
# ohne sie vorher manuell herunterzuladen und zu entpacken.
url = (
    "/vsizip/vsicurl/https://data.geo.admin.ch/ch.swisstopo.swissboundaries3d/"
    "swissboundaries3d_2026-01/swissboundaries3d_2026-01_2056_5728.gpkg.zip"
)
kantone = gpd.read_file(url, layer="tlm_kantonsgebiet")
print(kantone[["kantonsnummer", "name"]].head())

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