How to compute building shadows for a given time in Python
On this page
Problem statement
A shadow is a projection: take the building's outline, push it away from the sun by a distance that depends on the sun's elevation, and union the result with the building. The geometry is straightforward. The two things that make shadow analysis wrong are the sun position and the azimuth convention, and neither raises an error.
The scale of the effect is worth stating before any code. At The Hague, a 10 m building casts a 7.0 m shadow at midsummer noon and a 40.1 m shadow at midwinter noon — a factor of 5.7. Any analysis that picks "a sunny day" without specifying the date is not an analysis.
This guide computes shadows for a real footprint layer at real sun positions, and gets the conventions right.
Quick answer
import numpy as np, pandas as pd, pvlib, geopandas as gpd
from shapely.affinity import translate
from shapely.ops import unary_union
def shadow(geom, height, elevation_deg, azimuth_deg):
"""Azimuth is degrees clockwise from north; the shadow points away from the sun."""
if elevation_deg <= 0:
return None
length = height / np.tan(np.radians(elevation_deg))
dx = -length * np.sin(np.radians(azimuth_deg))
dy = -length * np.cos(np.radians(azimuth_deg))
moved = translate(geom, dx, dy)
return unary_union([geom, moved, geom.union(moved).convex_hull])
sp = pvlib.solarposition.get_solarposition(
pd.DatetimeIndex(["2026-12-21 12:00"]).tz_localize("Europe/Amsterdam"), 52.104, 4.2726)
elev, azim = sp.apparent_elevation.iloc[0], sp.azimuth.iloc[0]
shadows = [shadow(g, h, elev, azim) for g, h in zip(b.geometry, b.height_m)]pvlib gives apparent elevation — corrected for atmospheric refraction — which is what you want, because refraction lifts the apparent sun by about half a degree near the horizon and that changes a long shadow noticeably.
What you get — from The data used on this page
1,990 footprints, heights 0.5-25.1 m (median 6.7); extent 678 x 616 m
Library sun positions, aware timestamps, azimuth conventions
1. Get the sun position from a library, not a formula
Solar position depends on the date, the time, the time zone, the latitude, the longitude and the equation of time. pvlib.solarposition.get_solarposition handles all of it and returns both true and apparent elevation.
2. Use a timezone-aware timestamp
A naive timestamp is interpreted as UTC by most libraries, which on a summer afternoon in western Europe is two hours out. At 15:00 on 21 June in The Hague that moved the sun from azimuth 213.9° to 251.4° and from 57.9° to 43.3° elevation (Example 4).
3. Get the azimuth convention right
pvlib returns azimuth in degrees clockwise from north: 90° is east, 180° is south, 270° is west. The shadow points away from the sun, so the offset is (-L·sin(az), -L·cos(az)). Getting the sign wrong puts every shadow on the sunny side, which looks entirely plausible on a map.
4. Skip times when the sun is below the horizon
elevation <= 0 means no shadow, not an infinite one. At 08:00 on the winter solstice in The Hague the apparent elevation is −7.07°, and at 09:00 it is 0.86° — which already gives a 1,679 m shadow from the tallest building in the tile.
5. Build the shadow as a swept polygon
Translating the footprint gives the far end; the union of the original, the translation and the convex hull of the two gives the swept region. For a convex footprint that is exact; for a concave one the convex hull slightly overstates the shadow, and the alternative is to sweep each edge.
6. Union the shadows to get shaded ground
Overlapping shadows must be unioned before any area is computed, or the total double-counts.
7. Subtract the buildings if you want shaded open ground
A building's own footprint is not shaded ground; it is a building.
8. Accumulate over a day for a duration map
Rasterise each timestep's shaded area onto a grid and sum, which gives hours of shade per cell — far more useful than one instant.
Shaded area and hours of shade through a winter day
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), in the Dutch national grid with NAP heights. The setup takes each building's ground surface as its footprint and its recorded ridge height above that ground as height_m, as a GeoDataFrame b in EPSG:28992.
import json
import os
import geopandas as gpd
import numpy as np
import requests
from shapely.geometry import Polygon
if not os.path.exists("model.city.json"):
reply = requests.get("https://3d.bk.tudelft.nl/opendata/cityjson/3dcities/v2.0/DenHaag_01.city.json", timeout=300)
reply.raise_for_status()
open("model.city.json", "wb").write(reply.content)
d = json.load(open("model.city.json"))
V = np.array(d["vertices"]) * d["transform"]["scale"] + d["transform"]["translate"]
rows = []
for oid, o in d["CityObjects"].items():
attrs = o.get("attributes") or {}
for g in o.get("geometry", []):
if g["type"] != "Solid" or not g.get("semantics") or "AbsoluteRidgeHeight" not in attrs:
continue
surfaces = g["semantics"]["surfaces"]
for face, si in zip(g["boundaries"][0], g["semantics"]["values"][0]):
if si is not None and surfaces[si].get("type") == "GroundSurface":
ground = float(np.mean([V[i][2] for i in face[0]]))
rows.append({"id": oid, "height_m": attrs["AbsoluteRidgeHeight"] - ground,
"geometry": Polygon([V[i][:2] for i in face[0]])})
b = gpd.GeoDataFrame(rows, crs=28992)
b = b[b.is_valid & (b.height_m > 0)].reset_index(drop=True)
minx, miny, maxx, maxy = b.total_bounds
print(f"{len(b):,} footprints, heights {b.height_m.min():.1f}-{b.height_m.max():.1f} m (median {b.height_m.median():.1f}); "
f"extent {maxx - minx:.0f} x {maxy - miny:.0f} m")1,990 footprints, heights 0.5-25.1 m (median 6.7); extent 678 x 616 m
Example 1 — sun positions through a winter day
import pandas as pd, pvlib
lat, lon = 52.104, 4.2726
times = pd.date_range("2026-12-21 08:00", "2026-12-21 16:00", freq="1h",
tz="Europe/Amsterdam")
sp = pvlib.solarposition.get_solarposition(times, lat, lon)
print(sp[["apparent_elevation", "azimuth"]].round(2).to_string()) apparent_elevation azimuth
2026-12-21 08:00:00+01:00 -7.07 119.54
2026-12-21 09:00:00+01:00 0.86 131.09
2026-12-21 10:00:00+01:00 6.81 143.37
2026-12-21 11:00:00+01:00 11.38 156.49
2026-12-21 12:00:00+01:00 13.99 170.33
2026-12-21 13:00:00+01:00 14.41 184.50
2026-12-21 14:00:00+01:00 12.58 198.52
2026-12-21 15:00:00+01:00 8.68 211.93
2026-12-21 16:00:00+01:00 3.11 224.51
Note that the maximum elevation is at 13:00 local time, not 12:00 — solar noon is displaced by the longitude offset within the time zone and by the equation of time.
Example 2 — shaded area through the day
import numpy as np
from shapely.ops import unary_union
results = []
for ts, row in sp.iterrows():
if row.apparent_elevation <= 0:
continue
polys = [shadow(g, h, row.apparent_elevation, row.azimuth)
for g, h in zip(b.geometry, b.height_m)]
shade = unary_union([p for p in polys if p is not None])
results.append({"time": ts.strftime("%H:%M"),
"elevation": round(row.apparent_elevation, 2),
"azimuth": round(row.azimuth, 2),
"shadow_m2": round(shade.area),
"max_length_m": round(b.height_m.max() / np.tan(np.radians(row.apparent_elevation)), 1)})
print(pd.DataFrame(results).to_string(index=False)) time elevation azimuth shadow_m2 max_length_m
09:00 0.86 131.09 581672 1678.9
10:00 6.81 143.37 230305 209.9
11:00 11.38 156.49 201021 124.6
12:00 13.99 170.33 193688 100.6
13:00 14.41 184.50 199085 97.6
14:00 12.58 198.52 212335 112.4
15:00 8.68 211.93 239181 164.2
16:00 3.11 224.51 339680 460.8
The footprints span 678 m by 616 m, 0.42 km², and at noon 45.4% of that is either building or shadow (Example 4). At 09:00 the shadow area is 0.58 km², more than the whole tile, because shadows run far beyond it.
Example 3 — hours of shade per square metre
import numpy as np, rasterio
from rasterio.features import rasterize
from rasterio.transform import from_origin
def shade_hours(buildings, sun_positions, res=2.0):
minx, miny, maxx, maxy = buildings.total_bounds
w, h = int((maxx - minx) / res), int((maxy - miny) / res)
tr = from_origin(minx, maxy, res, res)
hours = np.zeros((h, w), dtype="float32")
steps = [t for t, r in sun_positions.iterrows() if r.apparent_elevation > 0]
dt = 1.0 if len(steps) < 2 else (steps[1] - steps[0]).total_seconds() / 3600
for ts in steps:
r = sun_positions.loc[ts]
polys = [shadow(g, ht, r.apparent_elevation, r.azimuth)
for g, ht in zip(buildings.geometry, buildings.height_m)]
polys = [p for p in polys if p is not None]
if not polys:
continue
mask = rasterize(polys, out_shape=(h, w), transform=tr, fill=0, default_value=1,
dtype="uint8")
hours += mask * dt
return hours, tr
hours, tr = shade_hours(b, sp)
print(f"median shade hours: {np.median(hours):.1f}, max {hours.max():.1f}")median shade hours: 4.0, max 8.0
Example 4 — the numbers quoted on this page
from shapely.geometry import box
site = dict(latitude=52.104, longitude=4.2726)
noon = pd.DatetimeIndex(["2026-06-21 12:00", "2026-03-21 12:00", "2026-12-21 12:00"]).tz_localize("Europe/Amsterdam")
for ts, e in pvlib.solarposition.get_solarposition(noon, **site).apparent_elevation.items():
print(f"{ts:%d %b} 12:00 elevation {e:5.2f}° shadow of a 10 m building {10 / np.tan(np.radians(e)):5.1f} m")
# A naive timestamp is read as UTC: 15:00 on 21 June, local against naive
aware = pvlib.solarposition.get_solarposition(pd.DatetimeIndex(["2026-06-21 15:00"]).tz_localize("Europe/Amsterdam"), **site)
naive = pvlib.solarposition.get_solarposition(pd.DatetimeIndex(["2026-06-21 15:00"]), **site)
print(f"15:00 local: azimuth {aware.azimuth.iloc[0]:.1f}°, elevation {aware.apparent_elevation.iloc[0]:.1f}°; "
f"naive (UTC): azimuth {naive.azimuth.iloc[0]:.1f}°, elevation {naive.apparent_elevation.iloc[0]:.1f}°")
# Refraction: pvlib's SPA correction at standard pressure and 12 °C
for true_el in (0.5, 5.0, 14.0):
apparent_el = true_el + pvlib.spa.atmospheric_refraction_correction(1013.25, 12, true_el, 0.5667)
print(f"true {true_el:4.1f}° -> apparent {apparent_el:5.2f}°: 10 m casts {10 / np.tan(np.radians(true_el)):7.1f} m "
f"at true, {10 / np.tan(np.radians(apparent_el)):7.1f} m at apparent")
tile = box(*b.total_bounds)
for ts in ("12:00", "13:00"):
row = sp.loc[f"2026-12-21 {ts}"]
shade = unary_union([p for p in (shadow(g, h, row.apparent_elevation, row.azimuth)
for g, h in zip(b.geometry, b.height_m)) if p is not None])
print(f"{ts}: {shade.intersection(tile).area / tile.area:.1%} of the {tile.area / 1e6:.2f} km² tile shaded or built, "
f"{shade.difference(tile).area:,.0f} m² of shadow outside it")21 Jun 12:00 elevation 55.06° shadow of a 10 m building 7.0 m
21 Mar 12:00 elevation 37.19° shadow of a 10 m building 13.2 m
21 Dec 12:00 elevation 13.99° shadow of a 10 m building 40.1 m
15:00 local: azimuth 213.9°, elevation 57.9°; naive (UTC): azimuth 251.4°, elevation 43.3°
true 0.5° -> apparent 0.92°: 10 m casts 1145.9 m at true, 626.0 m at apparent
true 5.0° -> apparent 5.16°: 10 m casts 114.3 m at true, 110.7 m at apparent
true 14.0° -> apparent 14.07°: 10 m casts 40.1 m at true, 39.9 m at apparent
12:00: 45.4% of the 0.42 km² tile shaded or built, 4,348 m² of shadow outside it
13:00: 46.8% of the 0.42 km² tile shaded or built, 3,894 m² of shadow outside it
Refraction and the azimuth sign
Why apparent elevation rather than true
Atmospheric refraction bends light, lifting the apparent position of the sun near the horizon far more than higher up: in pvlib's model, by 0.42° at a true elevation of 0.5° and by 0.07° at 14°. At low sun that difference is large in shadow terms: at a true 0.5° the apparent elevation is 0.92°, and the shadow of a 10 m building shortens from 1,146 m to 626 m; at 14° it is 40.1 m against 39.9 m (Example 4). Use apparent elevation; it is what an observer sees.
Why the azimuth sign is the classic bug
Azimuth in pvlib is clockwise from north, and the shadow points in the opposite direction. The offset (-L·sin(az), -L·cos(az)) encodes that. Using + gives shadows on the sunny side, which on a map looks like a rendering choice rather than an error, and survives review depressingly often.
Why shadow length is so sensitive near sunrise and sunset
Length is h / tan(elevation), which goes to infinity as elevation goes to zero. At 09:00 on the solstice the elevation is 0.86° and the tallest building in the tile casts 1,679 m; an hour later at 6.81° it casts 210 m. Any analysis that averages across the day has to decide what to do with the first and last hour, because they dominate every total.
Why hours-of-shade beats a single instant
A shadow at one moment tells you very little; a garden shaded at 09:00 and sunlit from 10:00 is a sunny garden. Accumulating a mask over the day, weighted by the timestep, gives hours of direct sun per cell, which is the quantity planning policy and daylight guidance actually use.
Flat ground, terrain and trees
- Flat ground is assumed. On a slope, the shadow runs along the terrain and needs a raster method.
- Terrain casts shadows too. A hill to the south shades everything behind it.
- Trees are not buildings and are often the dominant shade in residential streets.
- Convex hull overstates concave footprints slightly; sweep edges if it matters.
- Use a fine timestep near sunrise and sunset, where the length changes fastest.
- Time zones and DST. Always localise the timestamps.
- Check the CRS is projected. Metres in, metres out.
- The shadow of a pitched roof is not the shadow of a prism to the ridge.
Internal links
- Shadows point the wrong way or have the wrong length — the conventions, when they go wrong
- How to estimate rooftop solar potential in Python — the same sun positions, different question
- How to compute sky view factor from a surface model — diffuse rather than direct light
- How to extrude building footprints into 3D in Python — the input geometry
- Where building heights come from, and how wrong they are — the height the length depends on
- How to create a hillshade in Python — the raster equivalent for terrain
- How to calculate a viewshed in Python — the related line-of-sight problem
- How to rasterize a vector layer in Python — building the hours-of-shade grid
FAQ
How do I calculate a building's shadow in Python?
Translate the footprint by height / tan(elevation) in the direction away from the sun, and union it with the original and the convex hull of the two.
Where do I get the sun's position?
pvlib.solarposition.get_solarposition(times, lat, lon) with timezone-aware timestamps. Use apparent_elevation, which includes refraction.
Which way does the shadow point?
Away from the sun. With azimuth measured clockwise from north, the offset is (-L·sin(az), -L·cos(az)).
How much does the date matter?
Enormously. A 10 m building at The Hague casts 7.0 m at midsummer noon and 40.1 m at midwinter noon.
What do I do about sunrise and sunset?
Skip elevations at or below zero, and use a fine timestep near the horizon — at 0.86° elevation a building in the test cast a 1,679 m shadow.
Is a single-instant shadow map useful?
Rarely. Accumulate a shaded mask over the day to get hours of shade per cell, which is what daylight guidance is written in terms of.