How to Compute a Horizon Profile from a DEM in Python

On this page

Problem statement

PVGIS will return a horizon profile for any point, computed from its own 90 m elevation model. For a site in a narrow valley, a finer DEM can give a noticeably different skyline, and a local study may need horizons at hundreds of candidate points or from a surface model that includes forest. Computing the profile yourself is a short ray-marching loop: for each compass direction, walk outwards across the DEM and keep the steepest angle to the terrain.

Done properly it needs three details that are easy to skip — a metric grid, the Earth's curvature beyond a few kilometres, and enough radius to reach the ridges that matter.

Quick answer

For every azimuth, sample the DEM along a ray and take the maximum elevation angle:

import numpy as np, rasterio, requests, pandas as pd
from rasterio.warp import reproject, Resampling, transform as warp_xy
from rasterio.transform import from_origin

LAT, LON = 46.635, 8.594                     # Andermatt
CRS, RES, RADIUS = "EPSG:32632", 30.0, 20_000

url = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N46_00_E008_00_DEM/"
       "Copernicus_DSM_COG_10_N46_00_E008_00_DEM.tif")
(xs,), (ys,) = warp_xy("EPSG:4326", CRS, [LON], [LAT])
n = int(2 * RADIUS / RES) + 1
grid_transform = from_origin(xs - RADIUS, ys + RADIUS, RES, RES)
dem = np.full((n, n), np.nan, dtype="float32")
with rasterio.open(url) as src:
    reproject(rasterio.band(src, 1), dem, dst_transform=grid_transform, dst_crs=CRS,
              resampling=Resampling.bilinear)

def horizon(dem, row0, col0, azimuths, radius=RADIUS, step=RES, R=6_371_000, k=0.13):
    z0 = dem[row0, col0] + 1.5                          # eye height above ground
    d = np.arange(step, radius, step)
    drop = d ** 2 / (2 * R) * (1 - k)                   # curvature minus refraction
    out = []
    for az in azimuths:                                  # 0 = north, clockwise
        r = np.rint(row0 - d * np.cos(np.radians(az)) / step).astype(int)
        c = np.rint(col0 + d * np.sin(np.radians(az)) / step).astype(int)
        ok = (r >= 0) & (r < dem.shape[0]) & (c >= 0) & (c < dem.shape[1])
        z = dem[r[ok], c[ok]] - drop[ok]
        out.append(np.degrees(np.arctan((z - z0) / d[ok])).max())
    return np.clip(np.array(out), 0, None)

centre = n // 2
az = np.arange(0, 360, 7.5)
mine = pd.Series(horizon(dem, centre, centre, az), index=az)
print(f"site elevation {dem[centre, centre]:.0f} m; highest skyline {mine.max():.1f}° "
      f"at azimuth {mine.idxmax():.1f}° (0 = north)")
print(mine.loc[[0, 90, 157.5, 180, 202.5, 270]].round(1).to_string())

site elevation 1441 m; highest skyline 25.8° at azimuth 165.0° (0 = north)
0.0      10.4
90.0     17.5
157.5    25.2
180.0    23.3
202.5    14.0
270.0    15.7

The site sits at the bottom of a valley with ridges to the south; the skyline there is over 20°, higher than the noon sun in December.

Metric DEM, a reach from terrain, curvature and refraction

1. Reproject the DEM to a metric CRS around the site

Distances along the ray must be in metres. A local UTM zone is enough; bilinear resampling avoids stair-stepped ridges.

2. Choose the radius from the terrain

In the Alps, ridges 15–20 km away can matter. On a plain, 5 km is plenty. The cost of a larger radius is linear in the number of samples per ray.

3. Correct for curvature and refraction

At distance $d$ the Earth's surface drops $d^2/2R$ below a flat plane — 31 m at 20 km. Atmospheric refraction bends light back by about 13% of that. Ignoring both makes distant ridges look slightly too high.

4. Sample densely enough along each ray

A step equal to the cell size catches every ridge. Coarser steps can jump over a narrow crest and underestimate the horizon.

5. Compare with an independent profile

PVGIS's printhorizon is the obvious reference. Differences of a few degrees are expected from resolution; differences of tens of degrees mean a wrong azimuth convention or a wrong centre cell.

6. Feed the profile to the yield model

Either pass it to PVGIS as userhorizon, or remove the beam component in your own model whenever the sun is below the skyline.

The horizon angle is the steepest sight line along the ray, after the Earth has curved away.

A ray-marched skyline against PVGIS

Example 1 — against PVGIS's 90 m profile

j = requests.get("https://re.jrc.ec.europa.eu/api/v5_3/printhorizon", params={
    "lat": LAT, "lon": LON, "outputformat": "json"}, timeout=120).json()
pvgis = pd.DataFrame(j["outputs"]["horizon_profile"])
pvgis["az_north"] = (pvgis.A + 180) % 360            # PVGIS: 0 = south
pvgis = pvgis.drop_duplicates("az_north").set_index("az_north").H_hor.sort_index()

both = pd.DataFrame({"30 m DEM": mine, "PVGIS": pvgis.reindex(mine.index)}).dropna()
both["diff"] = both["30 m DEM"] - both["PVGIS"]
print(f"{len(both)} directions: mean difference {both['diff'].mean():+.1f}°, "
      f"largest {both['diff'].abs().max():.1f}° at azimuth {both['diff'].abs().idxmax():.1f}°")
print(both.loc[[90, 135, 180, 225, 270]].round(1).to_string())

48 directions: mean difference -0.0°, largest 2.8° at azimuth 150.0°
       30 m DEM  PVGIS  diff
90.0       17.5   16.8   0.7
135.0      18.4   19.9  -1.5
180.0      23.3   22.5   0.8
225.0       7.5    7.6  -0.1
270.0      15.7   15.7  -0.0

Forty lines of ray-marching reproduce PVGIS's skyline with no bias and a worst direction 2.8° off, towards the south-south-east where the nearest slopes are steepest and the finer DEM resolves them better. That is the check that makes the code trustworthy for the points PVGIS cannot serve — hundreds of candidate cells, or a surface model that includes forest.

Example 2 — what the profile does to a year of beam

import pvlib

times = pd.date_range("2023-01-01", "2024-01-01", freq="10min", tz="UTC", inclusive="left")
sun = pvlib.solarposition.get_solarposition(times, LAT, LON)
day = sun.apparent_elevation > 0
skyline = np.interp(sun.azimuth, np.append(az, 360), np.append(mine.values, mine.values[0]))
blocked = day & (sun.apparent_elevation < skyline)

print(f"daylight: {day.sum() / 6:,.0f} h; sun behind terrain: {blocked.sum() / 6:,.0f} h "
      f"({blocked.sum() / day.sum():.0%})")
monthly = blocked.groupby(times.month).sum() / day.groupby(times.month).sum()
print("share of daylight hours shaded, by month:")
print(" ".join(f"{m}:{v:.0%}" for m, v in monthly.items()))

daylight: 4,446 h; sun behind terrain: 1,320 h (30%)
share of daylight hours shaded, by month:
1:59% 2:27% 3:21% 4:24% 5:23% 6:23% 7:23% 8:24% 9:22% 10:20% 11:51% 12:68%

This counts hours, not energy — the blocked hours are the low-sun ones that carry least beam — but it shows the seasonal structure: winter mornings and afternoons disappear behind the ridges.

Example 3 — many sites at once

rows = []
for dy_km in (-4, -2, 0, 2, 4):
    for dx_km in (-4, 0, 4):
        r = centre - int(dy_km * 1000 / RES)
        c = centre + int(dx_km * 1000 / RES)
        h = horizon(dem, r, c, az, radius=15_000)
        rows.append({"north km": dy_km, "east km": dx_km, "elev m": round(float(dem[r, c])),
                     "mean horizon °": h.mean().round(1),
                     "south horizon °": h[az == 180][0].round(1)})
print(pd.DataFrame(rows).to_string(index=False))

 north km  east km  elev m  mean horizon °  south horizon °
       -4       -4    2150            14.2             23.7
       -4        0    2686            19.2             33.9
       -4        4    2390             9.5              3.2
       -2       -4    1602            20.1             16.1
       -2        0    2140            20.9             17.2
       -2        4    1835            21.6             15.7
        0       -4    2502            14.9              1.1
        0        0    1441            15.4             23.3
        0        4    2619             9.4              0.3
        2       -4    1773            17.0             31.0
        2        0    1772            17.8              9.7
        2        4    2062            18.8             26.8
        4       -4    1940            22.4             12.0
        4        0    1039            24.8             27.0
        4        4    2590            10.5              0.4

Within 4 km the mean horizon ranges from 9° on an exposed summit to 25° at the bottom of the valley, and the southern horizon — the one that decides winter — from 0.3° to 34°. For siting in mountainous terrain, the horizon belongs in the suitability score at the resolution of the candidate grid, not as one value per region.

Same valley, two DEMs, no bias between them and a worst direction 2.8° apart.

Why ray-marching is enough

The horizon in one direction is the maximum over a line of the angle to each terrain point. Nothing about it needs a viewshed algorithm: a viewshed answers "which cells are visible", the horizon only "how high is the skyline". One array lookup per step makes it fast enough for thousands of points.

Why the eye height matters a little

At the valley floor, adding 1.5 m changes the angle to a ridge 5 km away by 0.02°. It matters more for a site on a flat plain with low obstacles nearby, and for panels on a roof.

Why a surface model can overstate the horizon

Copernicus GLO-30 is a digital surface model: tree tops and roofs are part of it. For a site in a clearing, the forest edge is a real obstruction. For a site in open ground surrounded by forest at 2 km, the trees add a fraction of a degree. Knowing which model you used is part of reading the result.

Horizon is a property of the point, not of the region.

Azimuths, tile edges and nodata

Azimuth conventions: this code uses 0 = north; PVGIS 0 = south; pvlib 0 = north. Convert once.
DEM tile edges: a radius that crosses a tile boundary needs a mosaic; rasterio.merge handles it.
Nodata over lakes or at tile edges comes through as NaN; np.nanmax or a fill avoids a NaN horizon.
Resolution in azimuth: 7.5° steps can miss a narrow gap; use 1–2° where the gaps matter.
Performance: for a raster of horizons, vectorise over azimuths or use a GPU; a pure Python loop is fine for tens of points.

Internal links

Terrain shading and the horizon explained — why the skyline matters
How to model PV output with pvlib — where to apply the shading
How to reproject a raster with rasterio — the reprojection step
How to merge and mosaic rasters with rasterio — sites near tile edges
How to create a hillshade in Python — shading from a single sun position
How to build a solar suitability map — horizon as a scored factor

FAQ

How do I calculate a horizon profile from a DEM?

For each azimuth, step outwards along a ray, compute the elevation angle to each terrain sample after the curvature correction, and keep the maximum.

How far should the ray go?

Far enough to reach the ridges that can rise above a few degrees — 15–20 km in high mountains, a few kilometres elsewhere.

Do I need to correct for Earth curvature?

Beyond a few kilometres, yes: the drop is 31 m at 20 km, reduced by about 13% for refraction.

Why does my profile differ from PVGIS?

Resolution (30 m against 90 m), surface versus terrain model, and azimuth conventions.

Can I give my own horizon to PVGIS?

Yes, with the userhorizon parameter: a list of heights at equally spaced azimuths starting from north.

Is a viewshed the same as a horizon?

No. A viewshed says which cells are visible from a point; a horizon profile is only the skyline angle in each direction.