โš™๏ธ How-To 3D & buildings Raster & rasterio 3D & Buildings ยท step 11 of 20 Updated

How to estimate rooftop solar potential in Python

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.

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)
tilt = np.degrees(np.arccos(abs(normal[2])))
azimuth = (np.degrees(np.arctan2(normal[0], normal[1])) + 360) % 360   # 180 = south

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.

Scene showing a pitched roof with its two planes, normals, tilt and azimuth.
Every roof plane is its own site; a building-level average hides the north face.

Step-by-step solution

1. Get the roof surfaces from the semantics

An LoD2 model labels them. A real Hague tile had 3,004 RoofSurface faces across 1,990 buildings โ€” an average of 1.5 planes per building, which is what you would expect from a stock that is 63.5% flat-roofed.

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): 15% larger at 30ยฐ, 41% at 45ยฐ.

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.

Flow from roof planes through area, orientation filtering, irradiance and losses to annual yield.
Five multiplications, each with a default that is wrong by tens of percent for somebody.

Code examples

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%}")
print(f"total 3D roof area {planes.area_3d_m2.sum():,.0f} mยฒ, "
      f"footprint equivalent {planes.footprint_m2.sum():,.0f} 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)")

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, pandas as pd, numpy as np

def annual_poa(lat, lon, tilt, azimuth, tz="Europe/Amsterdam"):
    times = pd.date_range("2026-01-01", "2026-12-31 23:00", freq="1h", tz=tz)
    sp = pvlib.solarposition.get_solarposition(times, lat, lon)
    cs = pvlib.location.Location(lat, lon, tz=tz).get_clearsky(times)
    poa = pvlib.irradiance.get_total_irradiance(
        surface_tilt=tilt, surface_azimuth=azimuth,
        solar_zenith=sp.apparent_zenith, solar_azimuth=sp.azimuth,
        dni=cs.dni, ghi=cs.ghi, dhi=cs.dhi)
    return float(poa.poa_global.sum() / 1000)      # kWh/mยฒ per year, clear-sky

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

for tilt, azim, label in [(35, 180, "south, 35ยฐ"), (35, 90, "east, 35ยฐ"),
                          (35, 0, "north, 35ยฐ"), (5, 180, "flat")]:
    poa = annual_poa(52.104, 4.2726, tilt, azim)
    print(f"{label:14} {poa:7.0f} kWh/mยฒ/yr  โ†’  "
          f"{yield_kwh(1, poa):6.0f} kWh/mยฒ/yr of generation")

The clear-sky model overestimates a real year substantially โ€” western European sites lose roughly a third to cloud โ€” so replace get_clearsky with a typical meteorological year from PVGIS or a reanalysis before quoting numbers to anyone.

Explanation

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.

Table of rooftop solar parameters โ€” module efficiency, performance ratio, setbacks, usable flat-roof fraction, azimuth window and minimum plane area โ€” with typical ranges.
Six judgement calls between the geometry and the headline figure.

Edge cases or notes

  • 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.