Terrain Shading and the Horizon Explained

Problem statement

A solar resource map gives each grid cell one irradiance number, computed for a flat surface under an open sky. In a valley the sun rises late over one ridge and sets early behind another, and in winter it may not clear the southern ridge at all. On a slope the ground itself is tilted, so a field on a north-facing hillside is a panel pointing the wrong way before any panel is installed.

Both effects are properties of the terrain, both can be computed from a digital elevation model, and both are routinely left out of first-pass siting because the resource map already looks like an answer.

Quick answer

The horizon profile is the elevation angle of the skyline in every compass direction, seen from one point. Whenever the sun is below that angle, the beam component is blocked and only diffuse light arrives. PVGIS computes it from its own DEM:

import requests, pandas as pd

def horizon(lat, lon):
    j = requests.get("https://re.jrc.ec.europa.eu/api/v5_3/printhorizon", params={
        "lat": lat, "lon": lon, "outputformat": "json"}, timeout=120).json()
    h = pd.DataFrame(j["outputs"]["horizon_profile"])
    return h.set_index("A")["H_hor"], j["inputs"]["location"]["elevation"]

prof, elev = horizon(46.635, 8.594)          # Andermatt, Swiss Alps
print(f"Andermatt, {elev:.0f} m: {len(prof)} directions, azimuth 0 = south")
print(prof.loc[[-90.0, -45.0, 0.0, 45.0, 90.0]].rename("horizon °").to_string())
print(f"highest skyline {prof.max():.1f}° at azimuth {prof.idxmax():.1f}")

Andermatt, 1443 m: 49 directions, azimuth 0 = south
A
-90.0    16.8
-45.0    19.9
 0.0     22.5
 45.0     7.6
 90.0    15.7
highest skyline 26.4° at azimuth -22.5

PVGIS uses azimuth 0 for south, negative towards east. The sun at noon in late December in Andermatt stands about 20° above the horizon — which the southern skyline beats.

A horizon with enough reach, against the sun path

1. Compute the horizon from a DEM with enough reach

Mountains 20 km away can shade a site. The horizon calculation needs a DEM extending well beyond the site — PVGIS uses a 90 m DEM; for a local study, a 30 m DEM with a 20–30 km radius is typical.

2. Compare the sun path with the skyline

For every hour of the year, check whether the sun's elevation is above the horizon angle at the sun's azimuth. If it is not, remove the beam component and keep the diffuse.

3. Report the loss by month, not only by year

Horizon losses concentrate in winter, when the sun is low. An annual loss of a few percent can mean a December with almost no direct sun.

4. Account for the slope of the ground separately

A slope acts like a free tilt. A south-facing 10° slope helps a flat-mounted system; a north-facing one hurts it. For ground-mounted rows, the slope also changes row spacing and self-shading.

5. Keep near-field shading out of the terrain horizon

Trees, buildings and neighbouring rows shade too, but they are not in a DEM. A digital surface model (DSM) captures some of them; row-to-row shading needs a layout model.

Irradiation with and without the horizon

Example 1 — irradiation with and without the horizon

def monthly(lat, lon, usehorizon):
    j = requests.get("https://re.jrc.ec.europa.eu/api/v5_3/seriescalc", params={
        "lat": lat, "lon": lon, "startyear": 2020, "endyear": 2020,
        "angle": 35, "aspect": 0, "components": 1, "usehorizon": usehorizon,
        "outputformat": "json"}, timeout=120).json()
    h = pd.DataFrame(j["outputs"]["hourly"])
    h["month"] = h["time"].str[4:6].astype(int)
    h["beam"] = h["Gb(i)"]
    h["total"] = h["Gb(i)"] + h["Gd(i)"] + h["Gr(i)"]
    return h.groupby("month")[["beam", "total"]].sum() / 1000

for name, (lat, lon) in {"Andermatt": (46.635, 8.594), "Toledo": (39.86, -4.02)}.items():
    open_sky = monthly(lat, lon, 0)
    shaded = monthly(lat, lon, 1)
    loss = 1 - shaded.total / open_sky.total
    print(f"{name}: annual {open_sky.total.sum():,.0f} → {shaded.total.sum():,.0f} kWh/m² "
          f"({1 - shaded.total.sum() / open_sky.total.sum():.1%} lost)")
    print("  loss by month:", " ".join(f"{v:.0%}" for v in loss))

Andermatt: annual 1,617 → 1,473 kWh/m² (8.9% lost)
  loss by month: 59% 17% 2% 2% 0% 0% 0% 1% 2% 6% 40% 41%
Toledo: annual 2,021 → 2,015 kWh/m² (0.3% lost)
  loss by month: 3% 0% 0% 0% 0% 0% 0% 0% 0% 0% 0% 3%

Toledo's horizon is essentially open: 0.3% a year, and 3% in the two darkest months. In Andermatt an annual loss of 8.9% hides a winter in which the tilted panel loses 40–59% of its irradiation. A resource map that ignores the horizon would rate the two sites on their open-sky numbers.

Example 2 — the slope of the ground as a free tilt

def pv(lat, lon, angle, aspect):
    j = requests.get("https://re.jrc.ec.europa.eu/api/v5_3/PVcalc", params={
        "lat": lat, "lon": lon, "peakpower": 1, "loss": 14,
        "angle": angle, "aspect": aspect, "outputformat": "json"}, timeout=120).json()
    return j["outputs"]["totals"]["fixed"]["E_y"]

flat = pv(39.86, -4.02, 0, 0)
print(f"flat ground: {flat:,.0f} kWh/kWp")
for aspect, name in ((0, "south"), (-90, "east"), (90, "west"), (180, "north")):
    row = [f"{pv(39.86, -4.02, slope, aspect) / flat - 1:+.0%}" for slope in (5, 10, 20)]
    print(f"{name:>6}-facing  5°: {row[0]:>5}   10°: {row[1]:>5}   20°: {row[2]:>5}")

flat ground: 1,383 kWh/kWp
 south-facing  5°:   +5%   10°:   +9%   20°:  +14%
  east-facing  5°:   -0%   10°:   -1%   20°:   -2%
  west-facing  5°:   -0%   10°:   -1%   20°:   -3%
 north-facing  5°:   -5%   10°:  -11%   20°:  -24%

A module laid parallel to a 20° north-facing slope produces 24% less than one on flat ground in the same place. For a ground-mounted plant, racking can correct some of that, at the cost of wider row spacing — which is why slope and aspect belong in the suitability map even inside the workable range.

Example 3 — when the sun is behind the ridge

import numpy as np, pvlib

times = pd.date_range("2020-12-21 06:00", "2020-12-21 17:00", freq="30min", tz="UTC")
sun = pvlib.solarposition.get_solarposition(times, 46.635, 8.594)
az_pvgis = sun["azimuth"] - 180                     # convert to PVGIS convention
skyline = np.interp(az_pvgis, prof.index, prof.to_numpy())
up = sun["apparent_elevation"] > 0
visible = sun["apparent_elevation"] > skyline
print(f"21 December, Andermatt: sun above the flat horizon for "
      f"{up.sum() / 2:.1f} h, above the skyline for {visible.sum() / 2:.1f} h")
print(f"maximum solar elevation {sun['apparent_elevation'].max():.1f}°, "
      f"skyline to the south {prof.loc[0.0]:.1f}°")

21 December, Andermatt: sun above the flat horizon for 8.5 h, above the skyline for 2.5 h
maximum solar elevation 20.0°, skyline to the south 22.5°

On the shortest day the sun is up for 8.5 hours and clears the skyline for 2.5 of them, through the lower gap to the south-west; at noon it is 2.5° below the southern ridge. That is invisible in any irradiance product that assumes an open sky.

Only the beam is removed

Why only the beam is removed

When the sun is behind a ridge, the site still sees most of the sky, so diffuse light keeps arriving (slightly reduced by the portion of sky the ridge blocks). That is why horizon losses are much smaller in cloudy climates, where diffuse light is a large share of the total, and much larger at clear high-altitude sites in winter.

Why the slope effect is asymmetric

A south-facing slope tilts a flat module towards the sun, which helps until the tilt passes the optimum. A north-facing slope tilts it away, and the loss grows faster than the gain because low winter sun hits the tilted surface at a grazing angle. East and west slopes sit in between and shift production towards morning or afternoon.

Why 90 m DEMs understate local shading

A coarse DEM smooths ridges and fills narrow valleys. For a site at the foot of a steep slope, the true skyline can be several degrees higher than a 90 m model suggests, and a few degrees near the winter sun path decide whether December has any direct sun at all.

Azimuth conventions, curvature and snow

• Azimuth conventions differ: PVGIS uses 0 = south, pvlib and most GIS tools 0 = north. Convert once, at the boundary.
• Earth curvature and refraction matter for distant horizons (over 20–30 km) and are included in good horizon tools.
• Snow changes the albedo and so the reflected component on tilted panels in exactly the months where horizon losses are worst.
• Horizon files can be supplied to PVGIS (userhorizon) when a better local profile exists from a DSM or a site survey.
• Bifacial modules see more of the ground and the sky behind them, and are less sensitive to beam loss than monofacial ones.

Internal links

• How to compute a horizon profile from a DEM — doing it yourself with 30 m data
• Solar irradiance explained — beam and diffuse components
• How to find the optimal tilt and orientation — the same geometry, chosen rather than imposed
• How to create a hillshade in Python — shading from a single sun position
• How to calculate slope and aspect from a DEM — the terrain inputs
• How to estimate solar potential on roofs — near-field shading on buildings

FAQ

What is a horizon profile?

The elevation angle of the skyline in each compass direction from one point. The sun's beam is blocked whenever it is below that angle.

How much energy does terrain shading cost?

On open plains, nearly nothing. In mountain valleys, a modest share of the annual total and most of the winter beam.

Does shading remove all the light?

No — only the direct beam. Diffuse skylight still arrives, slightly reduced by the blocked part of the sky.

Which azimuth convention does PVGIS use?

0 is south, negative is east, positive is west. pvlib uses 0 for north and 180 for south.

Does the ground slope matter for a ground-mounted plant?

Yes. It acts as a tilt, it changes row spacing, and on north-facing ground it reduces yield substantially.

Why does my DEM horizon differ from a site survey?

The DEM resolution smooths ridges, and it does not contain trees or buildings. Use a DSM or measured profile for local obstructions.