Anleitung

Die Gipfel vor Ihrem Hotel zeichnen

Gäste fragen, welcher Berg welcher ist. Ein Panorama mit den Namen beantwortet das, auf Ihrer Website oder auf einer Tafel auf der Terrasse. Ein Python-Skript zeichnet es für jeden Ort: die Silhouette aus unseren Höhendaten, die Namen und Höhen der Gipfel aus OpenStreetMap.

Panorama von Mürren nach Osten und Süden: Eiger, Mönch und Jungfrau mit ihren Namen und Höhen

Zwei Bilder, gezeichnet mit dem Skript unten: das Dorf Mürren nach Osten und Süden, und die Bergstation des Niederhorns nach Süden. Auf dem Handy wischen Sie sie seitwärts.

Drei Schritte

  1. Erstellen Sie ein kostenloses Konto, und erstellen Sie im Dashboard einen Schlüssel vom Typ «Alle APIs». Behalten Sie ihn auf Ihrem Computer, wo das Skript läuft. Ihre Website braucht ihn nicht.
  2. Kopieren Sie das Skript unten in eine Datei panorama.py, setzen Sie Ihren Schlüssel dort ein, wo IHR_API_SCHLÜSSEL steht, und tragen Sie ein, wo Sie stehen (LAT, LON) und in welche Richtungen Sie schauen (FROM, TO).
  3. Starten Sie es mit Python 3 und dem Paket requests (pip install requests, dann python3 panorama.py). In weniger als einer Minute ist panorama.svg fertig.

Setzen Sie das Bild wie ein Foto auf Ihre Website, <img src="panorama.svg" alt="…">, oder öffnen Sie es in einem Browser oder Zeichenprogramm, um es zu drucken. SVG bleibt in jeder Grösse scharf, vom Handy bis zur Tafel von mehreren Metern. Die Website braucht keinen Schlüssel und schickt uns keine Anfragen.

Das Skript

# Zeichnet die Gipfel, die man von einem Ort aus sieht, als SVG-Bild: das Gelände
# aus der Höhen-API von overpass.ch, die Namen der Gipfel aus ihrer Overpass API.
import html, math, time, requests

KEY = "IHR_API_SCHLÜSSEL"
API = f"https://api.overpass.ch/v1/{KEY}"

LAT, LON = 46.5592, 7.8925  # wo Sie stehen: die Terrasse
EYE = 3                     # Meter über dem Boden
FROM, TO = 55, 205          # die Aussicht, als Himmelsrichtungen in Grad (über Norden: 300 bis 420)
RADIUS = 80_000             # wie weit der Blick reicht, in Metern
STEP = 0.25                 # Grad zwischen zwei Sichtlinien
STRETCH = 2                 # Höhen doppelt so hoch gezeichnet, wie sie wirken; 1 für den wahren Massstab
OUT = "panorama.svg"
BATCH = 10_000              # Punkte pro Anfrage: 1000 im kostenlosen Tarif

EARTH = 6_371_000
REFRACTION = 0.13  # die Luft biegt das Licht ein wenig um die Erde

def move(bearing, d):
    """Der Punkt d Meter von Ihrem Standort entfernt, in Richtung bearing."""
    la, lo, b, a = math.radians(LAT), math.radians(LON), math.radians(bearing), d / EARTH
    la2 = math.asin(math.sin(la) * math.cos(a) + math.cos(la) * math.sin(a) * math.cos(b))
    lo2 = lo + math.atan2(math.sin(b) * math.sin(a) * math.cos(la), math.cos(a) - math.sin(la) * math.sin(la2))
    return [round(math.degrees(la2), 6), round(math.degrees(lo2), 6)]

def bearing_and_distance(lat, lon):
    la1, lo1, la2, lo2 = map(math.radians, (LAT, LON, lat, lon))
    b = math.atan2(math.sin(lo2 - lo1) * math.cos(la2),
                   math.cos(la1) * math.sin(la2) - math.sin(la1) * math.cos(la2) * math.cos(lo2 - lo1))
    h = math.sin((la2 - la1) / 2) ** 2 + math.cos(la1) * math.cos(la2) * math.sin((lo2 - lo1) / 2) ** 2
    return math.degrees(b) % 360, 2 * EARTH * math.asin(math.sqrt(h))

def heights(points):
    """Höhen vieler Punkte, BATCH pro Anfrage; das Meer zählt als 0 m."""
    out = []
    for i in range(0, len(points), BATCH):
        r = requests.post(f"{API}/elevation", json={"locations": points[i:i + BATCH]}, timeout=120)
        while r.status_code == 429:  # mehr Anfragen pro Minute, als der Tarif erlaubt: warten
            time.sleep(int(r.headers.get("Retry-After", "5")))
            r = requests.post(f"{API}/elevation", json={"locations": points[i:i + BATCH]}, timeout=120)
        r.raise_for_status()
        out += [0 if res["surface"] == "sea" else res["elevation"] or 0 for res in r.json()["results"]]
    return out

# Jede Sichtlinie hat alle 90 m einen Punkt (das Raster der Daten) bis 10 km,
# dann jedes Mal 1 % weiter.
dists = list(range(90, 10_000, 90))
while dists[-1] < RADIUS:
    dists.append(round(dists[-1] * 1.01))
rays = [FROM + i * STEP for i in range(round((TO - FROM) / STEP) + 1)]

h = heights([[LAT, LON]] + [move(az, d) for az in rays for d in dists])
here = h[0] + EYE

def angle(height, d):
    """Grad über (oder unter) Augenhöhe, nach der Krümmung der Erde."""
    drop = d * d * (1 - REFRACTION) / (2 * EARTH)
    return math.degrees(math.atan2(height - here - drop, d))

# Für jede Sichtlinie der höchste Winkel bis zu jeder Distanz
sky = []
for i in range(len(rays)):
    row, top = [], -90
    for j, d in enumerate(dists):
        top = max(top, angle(h[1 + i * len(dists) + j], d))
        row.append(top)
    sky.append(row)

# Gipfel mit Namen aus OpenStreetMap, die über alles vor ihnen ragen und auf
# der Silhouette stehen
query = f'[out:json];node["natural"="peak"]["name"]["ele"](around:{RADIUS},{LAT},{LON});out;'
r = requests.post(f"{API}/interpreter", data={"data": query}, timeout=180)
r.raise_for_status()
seen = []
for p in r.json()["elements"]:
    try:
        ele = float(p["tags"]["ele"])
    except ValueError:
        continue
    az, d = bearing_and_distance(p["lat"], p["lon"])
    az += 360 if az < FROM else 0
    if az > TO or d < 500:
        continue
    i = round((az - FROM) / STEP)
    front = [a for a, x in zip(sky[i], dists) if x < d - 300]
    a = angle(ele, d)
    slack = math.degrees(30 / d)  # 30 m, für das grobe Raster bei Gipfeln
    if (not front or a > front[-1] - slack) and a > sky[i][-1] - 0.5:
        seen.append({"name": p["tags"]["name"], "height": f"{ele:,.0f} m".replace(",", "'"), "ele": ele, "az": az, "a": a})

# Das Bild: 10 Pixel pro Grad, die fernen Hügel hell, die nahen dunkel. Die
# Namen stehen in einer Reihe über dem höchsten Grat, 17 Pixel auseinander; die
# Gipfel, die am höchsten wirken, bekommen ihren zuerst, und ein Name weicht um
# bis zu zwei Plätze aus, wenn seiner besetzt ist.
PX = 10
W = round((TO - FROM) * PX)
x = lambda az: (az - FROM) * PX
placed = []
for p in sorted(seen, key=lambda p: -p["a"]):
    reach = (len(p["name"]) + len(p["height"])) * 3.5  # wie weit der schräge Name nach rechts reicht
    for shift in (0, -17, 17, -34, 34):
        p["x"] = x(p["az"]) + shift
        if all(abs(p["x"] - o["x"]) >= 17 for o in placed) and 0 <= p["x"] and p["x"] + reach < W:
            placed.append(p)
            break
ridge = max([r[-1] for r in sky] + [p["a"] for p in placed])
low = min(r[-1] for r in sky) - 6 / STRETCH
top = ridge * PX * STRETCH + 20 + max([(len(p["name"]) + len(p["height"])) * 6 for p in placed] + [0])
H = round(top - low * PX * STRETCH) + 24
y = lambda a: top - a * PX * STRETCH
row = y(ridge) - 12

svg = [f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {W} {H}" font-family="system-ui, sans-serif" font-size="12">',
       '<defs><linearGradient id="sky" x1="0" y1="0" x2="0" y2="1">'
       '<stop offset="0" stop-color="#d5e4f3"/><stop offset="1" stop-color="#f6f9fc"/></linearGradient></defs>',
       f'<rect width="{W}" height="{H}" fill="url(#sky)"/>']
for far, colour in [(RADIUS, "#c4d3e2"), (30_000, "#9db1c6"), (12_000, "#7a90a9"), (4_000, "#566b84"), (1_500, "#3a4d63")]:
    j = max(k for k, d in enumerate(dists) if d <= far)
    ridge_line = " ".join(f"{x(az):.1f},{y(sky[i][j]):.1f}" for i, az in enumerate(rays))
    svg.append(f'<polygon points="0,{H} {ridge_line} {W},{H}" fill="{colour}"/>')
for p in placed:
    X = x(p["az"])
    svg.append(f'<polyline points="{X:.1f},{y(p["a"]) - 3:.1f} {X:.1f},{row + 10:.1f} {p["x"]:.1f},{row:.1f}" fill="none" stroke="#1d2b3a" stroke-opacity=".4"/>')
    svg.append(f'<text transform="translate({p["x"] + 3:.1f},{row - 4:.1f}) rotate(-60)" fill="#1d2b3a">'
               f'{html.escape(p["name"])} <tspan fill="#5b6b7c">{p["height"]}</tspan></text>')
sides = {0: "N", 45: "NO", 90: "O", 135: "SO", 180: "S", 225: "SW", 270: "W", 315: "NW"}
for az in range(math.ceil(FROM / 10) * 10, TO + 1, 10):
    svg.append(f'<text x="{x(az):.1f}" y="{H - 8}" text-anchor="middle" fill="#ffffff">{sides.get(az % 360, f"{az % 360}°")}</text>')
svg.append(f'<text x="{W - 6}" y="{H - 24}" text-anchor="end" font-size="9" fill="#ffffff" fill-opacity=".8">'
           'Gelände: Copernicus WorldDEM™-90 © DLR e.V., © Airbus DS · Gipfel: © OpenStreetMap-Mitwirkende</text>')
svg.append("</svg>")
with open(OUT, "w", encoding="utf-8") as f:
    f.write("\n".join(svg))
print(f"{len(h)} Abfragen, {len(seen)} Gipfel auf der Silhouette, {len(placed)} mit Namen: {OUT}")

So funktioniert es

Von Ihrem Standort aus schickt das Skript jedes Viertelgrad eine Sichtlinie los und fragt die Höhe des Geländes entlang dieser Linie ab: alle 90 m bis 10 km, danach in grösseren Abständen bis 80 km. Auf jeder Linie ist der höchste Winkel über Ihren Augen die Silhouette in dieser Richtung. Je weiter weg das Gelände ist, desto mehr sinkt es hinter die Krümmung der Erde; das Skript rechnet damit, und mit der Brechung des Lichts in der Luft, die ferne Berge etwas anhebt.

Das Gelände wird in fünf Schichten nach Distanz gezeichnet, die fernen hell und die nahen dunkel, wie das Auge sie im Dunst sieht. Für die Namen holt eine Overpass-Abfrage jeden Gipfel in OpenStreetMap mit Namen und Höhe im Umkreis von 80 km. Ein Gipfel bekommt seinen Namen, wenn nichts vor ihm höher ist und er auf der Silhouette steht. Die Gipfel, die am höchsten wirken, kommen zuerst, und ein Name weicht um bis zu zwei Plätze aus, wenn sein Platz besetzt ist.

Einstellungen

  • FROM und TO sind Himmelsrichtungen: 90 ist Osten, 180 Süden, 270 Westen. Für eine Aussicht über Norden gehen Sie über 360 hinaus, zum Beispiel 300 bis 420. Ein ganzer Kreis ist 0 bis 360.
  • EYE ist die Höhe Ihrer Augen über dem Boden. Im dritten Stock eines Hotels setzen Sie etwa 10.
  • STRETCH zeichnet die Höhen doppelt so hoch, wie sie wirken, wie es gezeichnete Panoramen meist tun. Mit 1 hat das Bild die Proportionen eines Fotos, und die Berge wirken niedrig.
  • RADIUS ist, wie weit das Skript schaut. An klaren Tagen sieht das Auge mehr als 100 km; darüber hinaus wird das Bild voll mit Namen.
  • BATCH ist die Zahl der Punkte pro Anfrage. Der kostenlose Tarif nimmt 1000 und 30 Anfragen pro Minute: Setzen Sie 1000, dann wartet das Skript, wenn es muss, und braucht etwa 7 Minuten.
  • Für Namen in einer anderen Sprache ändern Sie p["tags"]["name"] in p["tags"].get("name:fr", p["tags"]["name"]).

Was es kostet

Ein Bild von 150 Grad braucht etwa 190'000 Höhenabfragen und eine Overpass-Anfrage, ein ganzer Kreis etwa 460'000. Das Skript läuft einmal pro Bild. Der kostenlose Tarif enthält 1 Million Höhenabfragen im Monat, für private und nicht kommerzielle Nutzung. Für ein Hotel, ein Restaurant oder eine Bergbahn kostet der Tarif Elevation €5 im Monat zuzüglich MWST und umfasst 50 Millionen Höhenabfragen, etwa 250 Bilder.

Gut zu wissen

  • Die Höhen der Gipfel kommen aus OpenStreetMap; die Silhouette kommt aus einem 90-m-Raster, in dem spitze Gipfel etwas tiefer liegen. In der Nähe des Standorts sind Gebäude und Bäume nur zum Teil in den Daten.
  • Das Bild zeigt jeden Berg, der in Sichtlinie liegt. Am Tag selbst können Dunst oder Wolken einige davon verdecken.
  • Ein Gipfel ohne Namen oder Höhe in OpenStreetMap bekommt keine Beschriftung. Jede Person kann sie auf openstreetmap.org ergänzen; lassen Sie das Skript ein paar Tage später noch einmal laufen.
  • Das Skript schreibt die Quellen in die Ecke des Bildes: die Höhendaten und OpenStreetMap. Lassen Sie sie stehen, wenn Sie es veröffentlichen.