How to compute sky view factor from a surface model

Problem statement

Sky view factor is the fraction of the hemisphere above a point that is sky rather than building or terrain. It runs from 1 in the open to near 0 in a narrow courtyard, and it is the standard proxy for urban heat island intensity, diffuse daylight, night-time cooling and street-level thermal comfort.

Computed from a surface model it is a horizon problem: in each of a set of directions, find the largest angle to the skyline, and integrate what is left. The two decisions that change the answer are how many directions you sample and how far you look, and neither has an obviously right value.

This guide computes SVF on a real derived DSM — 678 × 616 cells at 1 m over central The Hague, 19.5% built — and shows what the parameters do.

Quick answer

import numpy as np

def sky_view_factor(dsm, res, n_dirs=16, max_dist=80):
    """Fraction of sky visible, from horizon angles in n_dirs directions."""
    steps = int(max_dist / res)
    horizon = np.zeros((n_dirs,) + dsm.shape, dtype="float32")
    for k in range(n_dirs):
        a = 2 * np.pi * k / n_dirs
        dx, dy = np.cos(a), -np.sin(a)
        best = np.zeros(dsm.shape, dtype="float32")
        for s in range(1, steps + 1):
            sx, sy = int(round(dx * s)), int(round(dy * s))
            if sx == 0 and sy == 0:
                continue
            shifted = np.roll(np.roll(dsm, -sy, axis=0), -sx, axis=1)
            best = np.maximum(best, (shifted - dsm) / (s * res))
        horizon[k] = np.arctan(best)
    return (np.cos(horizon) ** 2).mean(axis=0)

On the test DSM, judged on the open ground at least 80 m from the edge (where np.roll's wrap-around cannot reach), this gave a median SVF of 0.839, a mean of 0.781 and a minimum of 0.009, with 12.1% of open ground below 0.5 and 35.8% below 0.7.

1,990 footprints, heights 0.5-25.1 m (median 6.7); extent 678 x 616 m

Surface model, resolution, directions, search distance

1. Build or obtain a surface model

A DSM from LiDAR is the right input. Without one, rasterise building heights onto a grid — that gives a built-form SVF that ignores trees, which is often what an urban morphology study wants anyway.

from rasterio.features import rasterize
from rasterio.transform import from_origin

dsm = rasterize(zip(b.geometry, b.height_m), out_shape=(h, w),
                transform=from_origin(minx, maxy, res, res), fill=0.0, dtype="float32")

2. Choose the resolution

1 m is the usual compromise for street-level work. Coarser grids close narrow streets: the share of interior open ground below 0.5 was 12.1% at 1 m, 9.8% at 2 m and 4.8% at 4 m, with the built fraction hardly changing (Example 4) — the deep, narrow spaces are exactly what a coarse grid loses. Finer grids cost four times as much per halving of the cell size.

3. Choose the number of directions

Eight directions leave gaps between the sampled azimuths that show as a star pattern around obstructions; 16 is the common default; 32 is smoother and cost 0.5 s against 0.2 s here. On this tile the distribution barely moved: at 80 m the interior median was 0.8410 with 8 directions, 0.8393 with 16 and 0.8380 with 32 (Example 2). The difference is in individual cells near corners, not in the summary.

4. Choose the search distance

Beyond a few times the tallest obstruction, further cells cannot raise the horizon. For a stock whose maximum height is 25 m, 80 m is generous. Too short a distance systematically overestimates SVF in wide streets with tall buildings at the end.

5. Use the correct integration formula

For a horizon angle γ in each of n equal sectors, the standard estimator is SVF = (1/n) Σ cos²(γ). Using 1 − sin(γ) or the mean of cos(γ) are both different quantities, and all three appear in the literature — say which you used.

6. Decide where you are measuring

SVF at ground level is what matters for thermal comfort and street daylight. SVF at roof level is what matters for solar. The same DSM gives both; the difference is whether you evaluate on built cells or only on open ground.

7. Report the distribution, not the mean

The mean of 0.781 on the test tile is dominated by open ground. The interesting numbers are the fractions below thresholds — 12.1% below 0.5 and 35.8% below 0.7 — which locate the courtyards and narrow streets.

Sky view factor end to end, and its parameters

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 — build the DSM and compute SVF end to end

import numpy as np, geopandas as gpd, rasterio
from rasterio.features import rasterize
from rasterio.transform import from_origin

def dsm_from_buildings(b, res=1.0, height_col="height_m"):
    minx, miny, maxx, maxy = b.total_bounds
    w = int(np.ceil((maxx - minx) / res))
    h = int(np.ceil((maxy - miny) / res))
    tr = from_origin(minx, maxy, res, res)
    dsm = rasterize(((g, v) for g, v in zip(b.geometry, b[height_col])),
                    out_shape=(h, w), transform=tr, fill=0.0, dtype="float32")
    return dsm, tr

dsm, tr = dsm_from_buildings(b, res=1.0)
print(f"DSM {dsm.shape[1]} x {dsm.shape[0]} at 1.0 m, "
      f"built fraction {(dsm > 0).mean():.1%}, max height {dsm.max():.1f} m")

svf = sky_view_factor(dsm, res=1.0, n_dirs=16, max_dist=80)
# np.roll wraps, so the outer 80 m see the far side of the tile: judge only the interior
edge = 80
interior = np.zeros(dsm.shape, bool)
interior[edge:-edge, edge:-edge] = True
ground = svf[(dsm == 0) & interior]
print(f"SVF on open ground, {interior.sum():,} interior cells: min {ground.min():.3f}, median {np.median(ground):.3f}, "
      f"mean {ground.mean():.3f}, max {ground.max():.3f}")
for thr in (0.5, 0.7, 0.9):
    print(f"  below {thr}: {(ground < thr).mean():6.1%}")

DSM 678 x 616 at 1.0 m, built fraction 19.5%, max height 25.1 m
SVF on open ground, 236,208 interior cells: min 0.009, median 0.839, mean 0.781, max 1.000
  below 0.5:  12.1%
  below 0.7:  35.8%
  below 0.9:  57.6%

Example 2 — how the parameters change the answer

import itertools, numpy as np

for n_dirs, max_dist in itertools.product((8, 16, 32), (40, 80, 160)):
    s = sky_view_factor(dsm, 1.0, n_dirs=n_dirs, max_dist=max_dist)
    g = s[(dsm == 0) & interior]
    print(f"{n_dirs:3d} dirs, {max_dist:3d} m: median {np.median(g):.4f}, "
          f"mean {g.mean():.4f}, below 0.7 {(g < 0.7).mean():.2%}")

  8 dirs,  40 m: median 0.8487, mean 0.7895, below 0.7 34.82%
  8 dirs,  80 m: median 0.8410, mean 0.7843, below 0.7 35.35%
  8 dirs, 160 m: median 0.8400, mean 0.7830, below 0.7 35.42%
 16 dirs,  40 m: median 0.8474, mean 0.7864, below 0.7 35.31%
 16 dirs,  80 m: median 0.8393, mean 0.7812, below 0.7 35.82%
 16 dirs, 160 m: median 0.8384, mean 0.7798, below 0.7 35.88%
 32 dirs,  40 m: median 0.8463, mean 0.7854, below 0.7 35.44%
 32 dirs,  80 m: median 0.8380, mean 0.7803, below 0.7 35.95%
 32 dirs, 160 m: median 0.8371, mean 0.7789, below 0.7 36.03%

Run this once on your own study area and pick the point where the numbers stop moving. Reporting the parameters is as important as reporting the values.

Example 3 — write it out, and sample it at points

import rasterio, geopandas as gpd

profile = {"driver": "GTiff", "height": svf.shape[0], "width": svf.shape[1],
           "count": 1, "dtype": "float32", "crs": b.crs, "transform": tr,
           "compress": "deflate", "predictor": 3, "tiled": True}
with rasterio.open("svf.tif", "w", **profile) as dst:
    dst.write(svf.astype("float32"), 1)
    dst.descriptions = ("sky view factor",)
    dst.update_tags(METHOD="16 directions, 80 m, cos^2 integration",
                    DSM="buildings rasterised at 1 m, no vegetation")

open_cells = np.argwhere((dsm == 0) & interior)[::5000]                   # every 5,000th open cell
xs, ys = rasterio.transform.xy(tr, open_cells[:, 0], open_cells[:, 1])
pts = gpd.GeoDataFrame(geometry=gpd.points_from_xy(xs, ys), crs=b.crs)
with rasterio.open("svf.tif") as src:
    pts["svf"] = [v[0] for v in src.sample(zip(pts.geometry.x, pts.geometry.y))]
print(len(pts), "points:", [round(float(v), 3) for v in pts.svf])

35 points: [1.0, 1.0, 0.956, 0.584, 0.504, 0.526, 0.97, 0.998, 1.0, 0.984, 0.725, 0.994, 0.983, 0.782, 0.985, 0.825, 0.625, 0.91, 0.903, 0.91, 0.866, 0.99, 0.632, 0.734, 0.889, 0.645, 0.562, 0.685, 0.942, 0.67, 0.685, 0.993, 0.53, 0.637, 0.67]

Putting the method in the tags matters more here than for most rasters, because three different SVF definitions are in common use and the values are not comparable between them.

Example 4 — edges, cost and resolution

import time

# What the wrap does: all open ground against the interior, and a zero-padded run
pad = 80
padded = sky_view_factor(np.pad(dsm, pad), 1.0, n_dirs=16, max_dist=80)[pad:-pad, pad:-pad]
for label, values in (("wrapped, every open cell", svf[dsm == 0]),
                      ("wrapped, interior only", svf[(dsm == 0) & interior]),
                      ("zero-padded, every open cell", padded[dsm == 0])):
    print(f"{label:30} median {np.median(values):.3f}  mean {values.mean():.3f}  "
          f"below 0.5 {(values < 0.5).mean():5.1%}  below 0.7 {(values < 0.7).mean():5.1%}")
outer = (dsm == 0) & ~interior
print(f"outer band: largest change from wrapping {np.abs(svf - padded)[outer].max():.3f}")

# Cost of the direction count
for n in (8, 16, 32):
    start = time.perf_counter()
    sky_view_factor(dsm, 1.0, n_dirs=n, max_dist=80)
    print(f"{n:2d} directions x 80 steps: {time.perf_counter() - start:.1f} s")

# Resolution: coarser grids close narrow streets
for res in (1.0, 2.0, 4.0):
    d, _ = dsm_from_buildings(b, res=res)
    e = int(80 / res)
    s_res = sky_view_factor(d, res, n_dirs=16, max_dist=80)[e:-e, e:-e]
    g = s_res[d[e:-e, e:-e] == 0]
    print(f"{res:.0f} m: built fraction {(d > 0).mean():5.1%}, open-ground median {np.median(g):.3f}, "
          f"below 0.5 {(g < 0.5).mean():5.1%}")

wrapped, every open cell       median 0.965  mean 0.851  below 0.5  7.8%  below 0.7 23.1%
wrapped, interior only         median 0.839  mean 0.781  below 0.5 12.1%  below 0.7 35.8%
zero-padded, every open cell   median 0.975  mean 0.854  below 0.5  7.8%  below 0.7 23.1%
outer band: largest change from wrapping 0.420
 8 directions x 80 steps: 0.1 s
16 directions x 80 steps: 0.2 s
32 directions x 80 steps: 0.5 s
1 m: built fraction 19.5%, open-ground median 0.839, below 0.5 12.1%
2 m: built fraction 19.5%, open-ground median 0.844, below 0.5  9.8%
4 m: built fraction 19.3%, open-ground median 0.856, below 0.5  4.8%

The cos-squared formula and the cross artefact

Why the cos² formula

For an unobstructed hemisphere, the fraction of sky in a sector with horizon angle γ is cos²(γ) under the assumption of uniform sky radiance weighted by the cosine of the zenith angle — which is the assumption that makes SVF proportional to diffuse irradiance received. Averaging over equal-azimuth sectors gives the whole-hemisphere estimate. Other definitions exist, including the simpler 1 − sin(γ), which corresponds to a uniform radiance sky with no cosine weighting.

Why the direction count shows up as a cross

With eight directions the sampled azimuths are the cardinal and diagonal ones, so a building corner that blocks the sky between two sampled directions is missed entirely. The artefact looks like a faint eight-armed star around obstructions. Sixteen directions reduces it below the noise for most purposes.

Why the search distance matters more than it looks

An obstruction at distance d and height h subtends atan(h/d). At 80 m a 25 m building subtends 17°, which removes about 9% of the sky in that sector — not negligible. Cutting the search at 40 m in a stock with 25 m buildings systematically raises SVF in open squares surrounded by tall blocks, which are exactly the locations a heat study cares about.

Why the shift-and-compare implementation is fast

Rolling the whole array by a fixed offset and taking a maximum is vectorised over every cell at once, so the cost is n_dirs × steps array operations rather than a per-cell ray cast. For a 678 × 616 grid with 16 directions and 80 steps that is 1,280 array passes, which ran in 0.2 s. Note that np.roll wraps at the edges, so the outer max_dist metres of the result see the far side of the tile — crop them, or pad the DSM first. Wrapping changed single edge cells by up to 0.420. Reporting every open cell instead of the interior moved the median from 0.839 to 0.965, but mostly because the tile's fringe is more open than its core: zero-padding gives 0.975 over the same cells. Whatever surrounds the tile is simply not in the data, so the honest statistic is the interior one.

Wrapping arrays, vegetation and terrain

• np.roll wraps. Pad the DSM by max_dist and crop the result.
• Vegetation is usually the dominant obstruction at street level and is absent from a building-only DSM.
• Terrain matters in hilly cities and belongs in the same surface.
• Evaluate on open ground. SVF on a rooftop cell is a different quantity.
• Nodata in a DSM propagates. Fill before computing.
• Report the definition. cos², 1 − sin and continuous integration give different numbers.
• Anisotropic SVF exists for directional sky radiance; this is the isotropic version.
• A coarse DSM closes narrow streets. Check the built fraction at your resolution.

Internal links

• How to compute building shadows for a given time in Python — direct rather than diffuse light
• How to estimate rooftop solar potential in Python — where roof-level SVF is the shading correction
• How to extrude building footprints into 3D in Python — producing the input
• LiDAR surfaces explained: DSM, DTM and CHM — the measured DSM
• How to calculate a viewshed in Python — the same horizon problem, one direction at a time
• How to rasterize a vector layer in Python — building the DSM
• How to create a hillshade in Python — another angle-from-a-surface calculation
• How to test whether a window can see a landmark — line of sight in 3D

FAQ

What is sky view factor?

The fraction of the sky hemisphere visible from a point: 1 in the open, near 0 in a narrow courtyard. It is the standard proxy for diffuse daylight and night-time cooling.

How do I compute it from a DSM?

Find the maximum horizon angle in each of n directions, then average $\cos^2$ of those angles across the directions.

How many directions should I use?

Sixteen is the common default. Eight leaves gaps between sampled azimuths around obstructions; thirty-two is smoother and cost more than twice as much here, for a median that moved in the third decimal place.

How far should I search?

Far enough that further cells cannot raise the horizon — roughly three times the tallest obstruction. For 25 m buildings, 80 m is generous.

What resolution do I need?

About 1 m for street-level work. At 4 m, the share of open ground below 0.5 fell from 12.1% to 4.8%: narrow streets close up and the SVF is wrong where it matters most.

Why do my edges look wrong?

np.roll wraps around the array. Pad the DSM by the search distance and crop the result.