overpass.choverpass.ch
EUR
Automatic (USD) CHFEURUSD
Guide

Draw the summits seen from your hotel

Guests ask which mountain is which. A panorama with the names answers that, on your website or on a board on the terrace. A Python script draws it for any place: the skyline from our elevation data, the names and heights of the summits from OpenStreetMap.

Panorama from Mürren to the east and south: Eiger, Mönch and Jungfrau with their names and heights

Two pictures made with the script below: the village of Mürren to the east and south, and the mountain station of the Niederhorn to the south. On a phone, swipe them sideways.

Three steps

  1. Create a free account and, in the dashboard, create a key of the type “All APIs”. Keep it on your computer, where the script runs. Your website does not need it.
  2. Copy the script below into a file panorama.py, put your key where it says YOUR_API_KEY, and set where you stand (LAT, LON) and the directions of the view (FROM, TO).
  3. Run it with Python 3 and the package requests (pip install requests, then python3 panorama.py). In less than a minute, panorama.svg is ready.

Put the picture on your website like a photo, <img src="panorama.svg" alt="…">, or open it in a browser or a drawing program to print it. SVG stays sharp at any size, from a phone to a board of several metres. The website needs no key and makes no requests to us.

The script

# Draws the summits seen from a place as an SVG picture: the terrain from the
# overpass.ch Elevation API, the names of the summits from its Overpass API.
import html, math, time, requests

KEY = "YOUR_API_KEY"
API = f"https://api.overpass.ch/v1/{KEY}"

LAT, LON = 46.5592, 7.8925  # where you stand: the terrace
EYE = 3                     # metres above the ground
FROM, TO = 55, 205          # the view, as compass directions in degrees (across north: 300 to 420)
RADIUS = 80_000             # how far to look, in metres
STEP = 0.25                 # degrees between two lines of sight
STRETCH = 2                 # heights drawn twice as high as they look; 1 for true scale
OUT = "panorama.svg"
BATCH = 10_000              # locations per request: 1,000 on the free plan

EARTH = 6_371_000
REFRACTION = 0.13  # the air bends the light a little around the Earth

def move(bearing, d):
    """The point d metres from where you stand, in the direction 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):
    """Heights of many points, BATCH per request; the sea counts as 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:  # more requests a minute than the plan allows: wait
            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

# Each line of sight has a point every 90 m (the grid of the data) up to
# 10 km, then 1% further each time.
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):
    """Degrees above (or below) eye level, after the curve of the Earth."""
    drop = d * d * (1 - REFRACTION) / (2 * EARTH)
    return math.degrees(math.atan2(height - here - drop, d))

# For each line of sight, the highest angle up to each distance
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)

# Named summits from OpenStreetMap that rise above everything in front of
# them and stand on the skyline
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, for the coarse grid near summits
    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", "ele": ele, "az": az, "a": a})

# The picture: 10 pixels a degree, the far hills light and the near ones
# dark. The names stand in one row above the highest ridge, 17 pixels apart;
# the summits that look highest get theirs first, and a name moves aside by
# up to two places when its own is taken.
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  # how far right the slanted name reaches
    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: "NE", 90: "E", 135: "SE", 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">'
           'Terrain: Copernicus WorldDEM™-90 © DLR e.V., © Airbus DS · Summits: © OpenStreetMap contributors</text>')
svg.append("</svg>")
with open(OUT, "w", encoding="utf-8") as f:
    f.write("\n".join(svg))
print(f"{len(h):,} lookups, {len(seen)} summits on the skyline, {len(placed)} named: {OUT}")

How it works

From where you stand, the script sends out a line of sight every quarter of a degree and looks up the height of the ground along it: every 90 m up to 10 km, then further apart up to 80 km. On each line, the highest angle above your eyes is the skyline in that direction. The farther the ground, the more it sinks behind the curve of the Earth; the script takes that into account, and the bending of light in the air, which lifts far mountains a little.

The ground is drawn in five layers by distance, the far ones light and the near ones dark, as the eye sees them in haze. For the names, an Overpass query fetches every summit in OpenStreetMap with a name and a height within 80 km. A summit gets its name when nothing in front of it is higher and it stands on the skyline. The ones that look highest come first, and a name moves aside by up to two places when its own place is taken.

Settings

  • FROM and TO are compass directions: 90 is east, 180 south, 270 west. For a view across north, go past 360, for example 300 to 420. A full circle is 0 to 360.
  • EYE is the height of your eyes above the ground. On the third floor of a hotel, set about 10.
  • STRETCH draws the heights twice as high as they look, as drawn panoramas usually do. With 1, the picture has the proportions of a photo, and the mountains look low.
  • RADIUS is how far the script looks. On clear days the eye sees more than 100 km; beyond that, the picture gets crowded with names.
  • BATCH is the number of points per request. The free plan takes 1,000 and 30 requests a minute: set 1000, and the script waits when it has to and needs about 7 minutes.
  • For names in another language, change p["tags"]["name"] to p["tags"].get("name:en", p["tags"]["name"]).

What it costs

A picture of 150 degrees uses about 190,000 lookups and one Overpass request; a full circle about 460,000. The script runs once per picture. The free plan includes 1 million lookups a month, for personal and non-commercial use. For a hotel, restaurant or cable car, the Elevation plan costs €5 a month plus VAT and covers 50 million lookups, about 250 pictures.

Good to know

  • The heights of the summits come from OpenStreetMap; the skyline comes from a 90 m grid, in which sharp summits read a little lower. Near the viewpoint, buildings and trees are only partly in the data.
  • The picture shows every mountain in line of sight. On the day, haze or clouds can hide some of them.
  • A summit without a name or height in OpenStreetMap gets no label. Anyone can add them on openstreetmap.org; run the script again a few days later.
  • The script writes the credits into the corner of the picture: the elevation data and OpenStreetMap. Keep them when you publish it.
Get a free API keyAn elevation profile for your tour pages