⚙️ How-To 3D & buildings Raster & rasterio 3D & Buildings · step 11 Notebook By · Updated

How to estimate rooftop solar potential in Python

On this page

Problem statement

Rooftop solar potential is four multiplications and a great deal of care about what each one means. Usable roof area times irradiance times panel efficiency times a performance ratio gives annual yield — and every term has a default that is wrong by 20% for somebody.

The one that catches GIS people is the first. A pitched roof's area is not its footprint: at 30° pitch the roof is 15% larger than the polygon it sits on, and at 45° it is 41% larger. Half of that extra area faces north and is worthless. Using the footprint as the roof area is therefore wrong twice, in opposite directions, by different amounts.

This guide computes usable area, orientation and irradiance from a real LoD2 model — a tile of The Hague's 3D city model — and a typical meteorological year for the same place from the EU's PVGIS service.

Quick answer

import numpy as np

def roof_area_3d(points):
    """Planar polygon area in 3D from the Newell normal."""
    n = np.zeros(3)
    for i in range(len(points)):
        a, b = points[i], points[(i + 1) % len(points)]
        n += np.cross(a, b)
    return np.linalg.norm(n) / 2, n / np.linalg.norm(n)

area, normal = roof_area_3d(surface_vertices)
normal = normal if normal[2] > 0 else -normal                            # point it up
tilt = np.degrees(np.arccos(abs(normal[2])))
azimuth = (np.degrees(np.arctan2(normal[0], normal[1])) + 360) % 360   # 180 = south
print(f"{area:.1f} m² roof plane, tilt {tilt:.1f}°, azimuth {azimuth:.0f}° | footprint {area * np.cos(np.radians(tilt)):.1f} m²")
8.3 m² roof plane, tilt 48.1°, azimuth 14° | footprint 5.5 m²

Compute area, tilt and azimuth per roof surface — not per building — because a gable roof has two faces with opposite azimuths and only one of them is worth panelling.

A pitched roof with its two planes, normals, tilt and azimuth. two planes, two answerssouth, usablenorth, notarea = footprint / cos(tilt)the multiplier15° → 1.04×30° → 1.15×45° → 1.41×60° → 2.00×1 / cos(tilt)Half the extra area faces the wrong way spatialworkflow.io
Every roof plane is its own site; a building-level average hides the north face.

Roof surfaces, 3D area, tilt and azimuth, filters

1. Get the roof surfaces from the semantics

An LoD2 model labels them, so the roof is a filter on the semantic surface type:

roof_faces, buildings = 0, set()
for oid, obj in d["CityObjects"].items():
    for g in obj.get("geometry", []):
        sem = g.get("semantics")
        if g["type"] == "Solid" and sem:
            for si in sem["values"][0]:
                if si is not None and sem["surfaces"][si]["type"] == "RoofSurface":
                    roof_faces += 1
                    buildings.add(oid)
print(f"{roof_faces:,} RoofSurface faces across {len(buildings):,} buildings — {roof_faces / len(buildings):.1f} per building")
3,004 RoofSurface faces across 1,990 buildings — 1.5 per building

One and a half faces per building is what you expect from a city with many flat roofs: 42% of these planes, and half the roof area, are flat (example 1).

2. Compute area in 3D, not from the footprint

The Newell method gives both the area and the normal of a planar polygon in one pass. The relationship to the footprint is area_3d = area_2d / cos(tilt):

for tilt in (15, 30, 45, 60):
    print(f"{tilt}° pitch: roof is {1 / np.cos(np.radians(tilt)) - 1:.1%} larger than its footprint")
15° pitch: roof is 3.5% larger than its footprint
30° pitch: roof is 15.5% larger than its footprint
45° pitch: roof is 41.4% larger than its footprint
60° pitch: roof is 100.0% larger than its footprint

3. Derive tilt and azimuth from the normal

Tilt is the angle between the normal and vertical; azimuth is the compass direction the normal points, with 180° meaning south in the northern hemisphere. A flat roof has tilt ≈ 0 and an undefined azimuth, which is correct — a flat roof is mounted at whatever tilt you choose.

4. Filter to usable surfaces

The usual criteria: area above a threshold (8–10 m² for a meaningful array), tilt below about 60°, and azimuth within roughly ±90° of south in the northern hemisphere. Then subtract an allowance for setbacks, chimneys, dormers and plant — commonly 20–30% of the gross plane.

5. Get irradiance for the plane, not for the horizontal

Horizontal irradiance underestimates a tilted south-facing plane and overestimates a north-facing one. pvlib transposes global horizontal irradiance onto an arbitrary plane; with a typical meteorological year, that is an hourly calculation over a year.

6. Account for shading

Neighbouring buildings shade roofs, especially in dense terraces and especially in winter. A sky view factor computed at roof level, or a full shadow accumulation, is the difference between a plausible estimate and a defensible one.

7. Convert to yield honestly

yield_kWh = usable_area_m2 × irradiance_kWh_m2 × module_efficiency × performance_ratio. Module efficiency is 0.18–0.22 for current silicon; performance ratio — inverter, wiring, temperature, soiling — is 0.75–0.85. Publish both.

The chain from roof planes to annual yield. five multiplications3D plane areanot the footprintorientationfiltertilt and azimuthsetbackallowance20–30%plane-of-arraykWh/m²from a TMY× efficiency ×PR0.20 × 0.80Each term has a default that is wrong by tens of percent for somebody spatialworkflow.io
Five multiplications, each with a default that is wrong by tens of percent for somebody.

Usable roof planes after setbacks

The data used on this page

The LoD2 tile of The Hague's 3D city model published as a CityJSON sample dataset (DenHaag_01.city.json, CityJSON 2.0, EPSG:7415), and a typical meteorological year for the tile's location from PVGIS 5.3 (European Commission Joint Research Centre).

import json
from pathlib import Path

import numpy as np
import pandas as pd
import requests

if not Path("DenHaag_01.city.json").exists():
    Path("DenHaag_01.city.json").write_bytes(requests.get(
        "https://3d.bk.tudelft.nl/opendata/cityjson/3dcities/v2.0/DenHaag_01.city.json", timeout=300).content)
d = json.loads(Path("DenHaag_01.city.json").read_text())
V = np.asarray(d["vertices"], dtype="float64") * d["transform"]["scale"] + d["transform"]["translate"]

if not Path("tmy.json").exists():
    Path("tmy.json").write_text(requests.get("https://re.jrc.ec.europa.eu/api/v5_3/tmy",
                                             params={"lat": 52.077, "lon": 4.315, "outputformat": "json"}, timeout=120).text)
hourly = pd.DataFrame(json.loads(Path("tmy.json").read_text())["outputs"]["tmy_hourly"])
tmy = pd.DataFrame({"ghi": hourly["G(h)"].values, "dni": hourly["Gb(n)"].values, "dhi": hourly["Gd(h)"].values},
                   index=pd.date_range("2026-01-01", periods=len(hourly), freq="1h", tz="UTC"))

def first_roof_face():
    for obj in d["CityObjects"].values():
        for g in obj.get("geometry", []):
            sem = g.get("semantics")
            if g["type"] == "Solid" and sem:
                for face, si in zip(g["boundaries"][0], sem["values"][0]):
                    if si is not None and sem["surfaces"][si]["type"] == "RoofSurface":
                        pts = V[face[0]]
                        n = np.cross(pts[1] - pts[0], pts[2] - pts[0])
                        if np.linalg.norm(n) > 0 and abs(n[2]) / np.linalg.norm(n) < 0.9:   # a pitched one
                            return pts
surface_vertices = first_roof_face()
print(d["version"], len(d["CityObjects"]), "city objects,", len(V), "vertices | TMY hours:", len(tmy),
      f"| annual GHI {tmy.ghi.sum() / 1000:,.0f} kWh/m²")
2.0 2498 city objects, 22997 vertices | TMY hours: 8760 | annual GHI 1,087 kWh/m²

Example 1 — roof planes with area, tilt and azimuth

import numpy as np, pandas as pd

def roof_planes(city_objects, V):
    rows = []
    for oid, obj in city_objects.items():
        for g in obj.get("geometry", []):
            sem = g.get("semantics")
            if g["type"] != "Solid" or not sem:
                continue
            surfaces, values = sem["surfaces"], sem["values"][0]
            for face, si in zip(g["boundaries"][0], values):
                if si is None or surfaces[si].get("type") != "RoofSurface":
                    continue
                pts = np.array([V[i] for i in face[0]])
                if len(pts) < 3:
                    continue
                n = np.zeros(3)
                for i in range(len(pts)):
                    n += np.cross(pts[i], pts[(i + 1) % len(pts)])
                area = np.linalg.norm(n) / 2
                if area <= 0:
                    continue
                u = n / np.linalg.norm(n)
                if u[2] < 0:
                    u = -u
                tilt = float(np.degrees(np.arccos(np.clip(abs(u[2]), 0, 1))))
                azim = float((np.degrees(np.arctan2(u[0], u[1])) + 360) % 360)
                rows.append({"id": oid, "area_3d_m2": float(area),
                             "footprint_m2": float(area * np.cos(np.radians(tilt))),
                             "tilt_deg": tilt, "azimuth_deg": azim if tilt > 5 else np.nan})
    return pd.DataFrame(rows)

planes = roof_planes(d["CityObjects"], V)
print(f"{len(planes):,} roof planes over {planes.id.nunique():,} buildings")
print(f"flat (tilt < 5°): {(planes.tilt_deg < 5).mean():.1%} of planes, "
      f"{planes.area_3d_m2[planes.tilt_deg < 5].sum() / planes.area_3d_m2.sum():.1%} of roof area")
print(f"total 3D roof area {planes.area_3d_m2.sum():,.0f} m², "
      f"footprint equivalent {planes.footprint_m2.sum():,.0f} m²")
3,004 roof planes over 1,990 buildings
flat (tilt < 5°): 42.4% of planes, 51.3% of roof area
total 3D roof area 94,437 m², footprint equivalent 82,609 m²

Example 2 — usable area after filtering and a setback allowance

import numpy as np

def usable(planes, min_area=10, max_tilt=60, south_window=90, setback=0.25,
           hemisphere="north"):
    p = planes.copy()
    flat = p.tilt_deg < 5
    ideal_azimuth = 180 if hemisphere == "north" else 0
    delta = (p.azimuth_deg - ideal_azimuth).abs()
    delta = np.minimum(delta, 360 - delta)

    keep = (p.area_3d_m2 >= min_area) & (p.tilt_deg <= max_tilt) & (flat | (delta <= south_window))
    p["usable_m2"] = np.where(keep, p.area_3d_m2 * (1 - setback), 0.0)
    return p

u = usable(planes)
print(f"planes kept: {(u.usable_m2 > 0).sum():,} of {len(u):,}")
print(f"usable area: {u.usable_m2.sum():,.0f} m² "
      f"({u.usable_m2.sum() / u.area_3d_m2.sum():.1%} of gross roof area)")
planes kept: 1,379 of 3,004
usable area: 50,510 m² (53.5% of gross roof area)

The setback is the single largest judgement in the whole calculation and the one most often left out. State it.

Example 3 — plane-of-array irradiance and annual yield

import pvlib

LAT, LON = 52.077, 4.315                                   # the tile, in WGS 84

def annual_poa(tilt, azimuth, weather):
    sp = pvlib.solarposition.get_solarposition(weather.index, LAT, LON)
    poa = pvlib.irradiance.get_total_irradiance(
        surface_tilt=tilt, surface_azimuth=azimuth,
        solar_zenith=sp.apparent_zenith, solar_azimuth=sp.azimuth,
        dni=weather.dni, ghi=weather.ghi, dhi=weather.dhi)
    return float(poa.poa_global.sum() / 1000)              # kWh/m² per year

def yield_kwh(area_m2, poa_kwh_m2, efficiency=0.20, performance_ratio=0.80):
    return area_m2 * poa_kwh_m2 * efficiency * performance_ratio

clear = pvlib.location.Location(LAT, LON).get_clearsky(tmy.index)
print(f"{'plane':12} {'clear-sky':>10} {'TMY':>7} {'ratio':>6}  kWh/m²/yr generated (TMY)")
for tilt, azim, label in [(35, 180, "south, 35°"), (35, 90, "east, 35°"), (35, 0, "north, 35°"), (5, 180, "flat")]:
    ideal, real = annual_poa(tilt, azim, clear), annual_poa(tilt, azim, tmy)
    print(f"{label:12} {ideal:10.0f} {real:7.0f} {real / ideal:6.0%}  {yield_kwh(1, real):6.0f}")
plane         clear-sky     TMY  ratio  kWh/m²/yr generated (TMY)
south, 35°         2029    1237    61%     198
east, 35°          1432    1005    70%     161
north, 35°          813     747    92%     120
flat               1655    1127    68%     180

The clear-sky model is an upper bound: it assumes no cloud all year. Against PVGIS's typical year for The Hague, a south-facing plane received 61% of its clear-sky figure — clear sky overstated it by nearly two thirds — which is why the TMY goes into any number quoted to anyone. The north-facing plane kept 92% of its (small) clear-sky figure, because on a cloudy day it loses little direct sun it was never getting, and diffuse light arrives from the whole sky.

The footprint is not the roof

Why the footprint is not the roof area

A plane tilted at angle θ projects onto a horizontal footprint of area × cos(θ). Inverting that, the roof is footprint / cos(θ) — 1.15× at 30°, 1.41× at 45°. On a gable roof both faces have the same tilt, so the total roof area is the footprint divided by the cosine, split equally between two opposite azimuths.

Why per-plane rather than per-building

A gable roof running east–west has a south face that is excellent and a north face that is nearly useless; averaging them gives a building that looks mediocre and is not. A building-level number also hides the case that matters most for planning — a large flat roof, which is unconstrained in azimuth and can be tilted optimally.

Why flat roofs are a different calculation

A flat roof has no inherent tilt or azimuth, so panels are mounted on frames at a chosen tilt. That means the usable fraction is set by row spacing rather than by the roof shape: rows must be far enough apart not to shade each other in winter, which typically limits coverage to 40–60% of the roof. Applying a pitched-roof setback to a flat roof overstates it considerably.

Why shading is the term people skip

Irradiance models assume an unobstructed horizon. In a terrace, in a courtyard, or next to anything taller, that is false for a large part of the year, and the error is largest in winter when the sun is low. A sky view factor computed at roof height is a cheap correction; a full shadow accumulation is the honest one.

The parameters in a rooftop solar estimate and their usual ranges. publish every termtermtypicalwhy it variesmodule efficiency0.18–0.22technology and ageperformance ratio0.75–0.85inverter, heat, soilingsetback, pitched20–30%chimneys, dormers, edgesusable, flat roof40–60%row spacing, not setbacksazimuth window±90° of southa policy choiceminimum plane area8–10 m²a meaningful arraythe setback is the largest judgement and the one most often omittedNone of these is a property of the roof spatialworkflow.io
Six judgement calls between the geometry and the headline figure.

Flat-roof azimuth, inward normals and dormers

  • Azimuth is undefined for a flat roof. Set it to NaN rather than 0.
  • Normals may point inwards. Flip so the z component is positive.
  • Holes reduce roof area. Rooflights and courtyards are interior rings.
  • Dormers are separate small planes and usually below the area threshold.
  • Structural capacity limits retrofits and is not in any GIS layer.
  • Heritage and planning constraints exclude roofs the geometry says are fine.
  • Clear-sky irradiance is an upper bound. Use a TMY for anything quoted.
  • Panels have a minimum sensible array size. A 6 m² plane is rarely worth it.

FAQ

Is a roof's area the same as its footprint?

No. A tilted plane has area footprint / cos(tilt) — 15% more at 30° and 41% more at 45°.

Should I compute solar potential per building or per roof plane?

Per plane. A gable roof has a good south face and a poor north face, and averaging them describes neither.

How do I get tilt and azimuth from a roof surface?

From the surface normal: tilt is the angle from vertical, azimuth is atan2(nx, ny) in degrees, with 180° meaning south.

What fraction of a roof is usable?

Typically 70–80% of a pitched plane after setbacks, chimneys and dormers, and 40–60% of a flat roof once row spacing is allowed for.

Which irradiance should I use?

Plane-of-array, transposed from a typical meteorological year. Clear-sky irradiance is an upper bound and overstates western European sites substantially.

Does shading matter much?

Yes, and most in winter when the sun is low. An unshaded irradiance model on a dense terrace is optimistic by a wide margin.