How to compute building shadows for a given time in Python
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.
Step-by-step solution
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 early and moves the sun by 30°.
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.
Code examples
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 tile is 787 m by 672 m, about 0.53 km², so at noon roughly 37% of it is in shadow. At 09:00 the shadow area exceeds a third of a square kilometre because shadows extend far beyond the tile.
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}")Explanation
Why apparent elevation rather than true
Atmospheric refraction bends light, lifting the apparent position of the sun by about 34 arcminutes at the horizon and by a few arcminutes higher up. At low sun that difference is large in shadow terms: at a true elevation of 0.5°, apparent elevation is closer to 1°, and the shadow of a 10 m building shortens from 1,146 m to 573 m. 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.
Edge cases or notes
- 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.