How to test whether a window can see a landmark

Problem statement

A viewshed answers "what can be seen from this point over this terrain". Urban visibility is a different question: can a specific window, at a specific height on a specific facade, see a specific thing — and the obstructions are buildings, not hills.

The answer is a line-of-sight test in three dimensions: cast a ray from the observer to the target and check whether any building geometry lies across it. It is also the basis for protected-view analysis, for valuing an outlook, and for testing whether a proposed building blocks something it should not.

This guide does it two ways — exactly against meshes, and quickly against a surface model — and covers the cases where the two disagree.

Quick answer

import numpy as np, trimesh

scene = trimesh.util.concatenate(building_meshes)
intersector = trimesh.ray.ray_triangle.RayMeshIntersector(scene)

def visible(observer_xyz, target_xyz, epsilon=0.5):
    direction = np.asarray(target_xyz) - np.asarray(observer_xyz)
    distance = np.linalg.norm(direction)
    direction = direction / distance
    origin = np.asarray(observer_xyz) + direction * epsilon   # step off our own facade
    hits = intersector.intersects_location([origin], [direction])[0]
    if len(hits) == 0:
        return True
    nearest = np.linalg.norm(hits - origin, axis=1).min()
    return nearest >= distance - epsilon

The epsilon step-off is essential. An observer placed exactly on a facade hits its own building at distance zero and is reported as seeing nothing.

1,990 meshes, 1,990 single-polygon footprints

Observers on facades at storey height, then cast rays

1. Place the observer properly

A window is on a facade, at a floor height, facing outwards. Generate observers by walking the footprint boundary, offsetting outwards by half a metre along the outward normal, and setting z to ground plus floor height. Check the ring's direction first: a CityJSON ground surface faces downwards, so seen from above its ring runs clockwise — all 1,990 on the Hague tile do — and the "outward" normal of a clockwise ring points into the building. orient(poly, 1) makes every ring counter-clockwise before the walk.

from shapely.geometry.polygon import orient


def facade_points(poly, spacing=3.0, offset=0.5):
    line = orient(poly, 1).exterior                  # force CCW so the normal points out
    for d in np.arange(0, line.length, spacing):
        p = line.interpolate(d)
        q = line.interpolate(min(d + 0.1, line.length))
        nx, ny = (q.y - p.y), -(q.x - p.x)          # outward for a CCW ring
        n = np.hypot(nx, ny) or 1.0
        yield (p.x + offset * nx / n, p.y + offset * ny / n)

2. Set the observer height from the storey, not the building

Floor 3 of a building whose ground is at 5.96 m NAP with 3.20 m storeys — the figure calibrated for these streets in the heights guide — puts the eye at 5.96 + 3 × 3.20 + 1.5 ≈ 17.1 m. Using the building's top height puts every observer on the roof.

3. Decide what counts as the target

A landmark is not a point. For a tower, a small set of points on its silhouette is a better target than its centroid, and "visible" then becomes "at least one target point visible" or "more than n of them", which is a decision to make explicitly.

4. Choose mesh or raster

• Mesh ray casting is exact, handles overhangs, and costs a ray per pair.
• DSM profile sampling is far faster, works for millions of pairs, and cannot represent anything the surface model cannot.

For a handful of protected views, use meshes. For "which of 40,000 dwellings can see the sea", use the DSM.

5. Include the terrain

In anything but flat ground, the terrain blocks more views than the buildings do. A combined surface — DTM plus building heights — is the right input for the raster method.

6. Add the atmosphere for long views only

Earth curvature drops a target by about 0.0785 d² metres for d in kilometres, reduced to roughly 0.87 of that by standard refraction. Below about 5 km it is under 2 m and usually ignorable; at 20 km it is 27 m and dominant.

7. Report a fraction, not a boolean

"42% of the facade points on this building can see the target" is more useful than a yes or no, and is what a valuation or a planning argument actually needs.

Which windows can see a landmark

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, 1,990 solids). The setup parses it (d, and V with the transform applied), fan-triangulates every solid into a trimesh mesh (meshes, with their CityJSON ids in ids), builds footprints from each solid's ground surface with its recorded eaves height, with its ground level and ridge height (b, whose index is the CityJSON id), and passes every mesh to the scene as building_meshes. The target in Example 1 is a stated point, 40 m above NAP near the tile's eastern edge — the analysis does not depend on what stands there.

import json
import os

import geopandas as gpd
import numpy as np
import requests
import trimesh
from shapely.geometry import Polygon
from shapely.ops import unary_union

if not os.path.exists("DenHaag_01.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("DenHaag_01.city.json", "wb").write(reply.content)
d = json.load(open("DenHaag_01.city.json"))
V = np.array(d["vertices"]) * d["transform"]["scale"] + d["transform"]["translate"]
vertices = V

meshes, ids, rows, geometry = [], [], [], None
for oid, o in d["CityObjects"].items():
    for g in o.get("geometry", []):
        if g["type"] != "Solid":
            continue
        geometry = geometry or g
        tris = [[ring[0], ring[i], ring[i + 1]]
                for shell in g["boundaries"] for surface in shell
                for ring in surface[:1] for i in range(1, len(ring) - 1)]
        m = trimesh.Trimesh(vertices=V, faces=np.array(tris), process=True)
        m.remove_unreferenced_vertices()
        meshes.append(m)
        ids.append(oid)
        sem = g["semantics"]
        ground = [Polygon(V[s[0]][:, :2]) for shell, vals in zip(g["boundaries"], sem["values"])
                  for s, v in zip(shell, vals) if v is not None and sem["surfaces"][v]["type"] == "GroundSurface"]
        rows.append({"id": oid, "height_eaves_m": o["attributes"]["RelativeEavesHeight"],
                     "height_m": o["attributes"]["RelativeRidgeHeight"],
                     "ground_z": o["attributes"]["AbsoluteRidgeHeight"] - o["attributes"]["RelativeRidgeHeight"],
                     "geometry": unary_union(ground)})
footprints = gpd.GeoDataFrame(rows, crs=28992).set_index("id")
mesh = next(m for m in meshes if not m.is_watertight).copy()
b = footprints[footprints.geom_type == "Polygon"]                # the table the examples walk
building_meshes = meshes
print(f"{len(meshes):,} meshes, {len(b):,} single-polygon footprints")

1,990 meshes, 1,990 single-polygon footprints

Example 1 — facade observers, and the fraction that can see a target

import numpy as np, geopandas as gpd, trimesh

def facade_observers(b, storey, storey_height=3.20, eye=1.5, spacing=3.0, offset=0.5):
    rows = []
    for row in b.itertuples():
        z = row.ground_z + storey * storey_height + eye
        if z > row.ground_z + row.height_m:
            continue                                  # storey does not exist
        for x, y in facade_points(row.geometry, spacing, offset):
            rows.append({"id": row.Index, "x": x, "y": y, "z": z})
    return gpd.GeoDataFrame(rows, geometry=gpd.points_from_xy(
        [r["x"] for r in rows], [r["y"] for r in rows]), crs=b.crs)

obs = facade_observers(b, storey=3)
target = np.array([78_900.0, 458_050.0, 40.0])

vis = np.array([visible((o.x, o.y, o.z), target) for o in obs.itertuples()])
obs["visible"] = vis
by_building = obs.groupby("id")["visible"].mean()
print(f"{len(obs):,} facade observers on {obs.id.nunique():,} buildings")
print(f"buildings with any view: {(by_building > 0).sum():,}")
print(f"buildings with more than half their facade in view: {(by_building > 0.5).sum():,}")
inside = obs.within(gpd.GeoSeries(b.loc[obs.id].geometry.values, index=obs.index))
print(f"observers inside their own footprint: {int(inside.sum())}")

2,436 facade observers on 170 buildings
buildings with any view: 148
buildings with more than half their facade in view: 62
observers inside their own footprint: 3

On the Hague tile, 2,436 third-floor observers on 170 buildings cast rays at the target; 148 buildings have at least one view and 62 have more than half their facade in view. The three observers left inside their own footprint sit in concave corners, where a half-metre step along the normal lands inside a neighbouring wall.

Example 2 — the fast raster method

import numpy as np, rasterio

def visible_on_surface(surface, transform, obs_xyz, target_xyz, step=1.0):
    """Sample the surface along the line and compare against the sight line."""
    o, t = np.asarray(obs_xyz, float), np.asarray(target_xyz, float)
    d = np.linalg.norm(t[:2] - o[:2])
    n = max(int(d / step), 2)
    xs = np.linspace(o[0], t[0], n)
    ys = np.linspace(o[1], t[1], n)
    zs_line = np.linspace(o[2], t[2], n)

    inv = ~transform
    cols, rows = inv * (xs, ys)
    rows, cols = rows.astype(int), cols.astype(int)
    ok = (rows >= 0) & (rows < surface.shape[0]) & (cols >= 0) & (cols < surface.shape[1])
    z_surface = np.full(n, -np.inf)
    z_surface[ok] = surface[rows[ok], cols[ok]]

    interior = slice(1, -1)                            # ignore the end points
    return bool(np.all(z_surface[interior] <= zs_line[interior]))


# a 1 m surface: ground plus ridge height inside each footprint, ground level elsewhere
from rasterio.features import rasterize
x0, y0, x1, y1 = b.total_bounds
transform = rasterio.transform.from_origin(np.floor(x0), np.ceil(y1), 1.0, 1.0)
shape = (int(np.ceil(y1) - np.floor(y0)), int(np.ceil(x1) - np.floor(x0)))
surface = rasterize(((g, z) for g, z in zip(b.geometry, b.ground_z + b.height_m)),
                    out_shape=shape, transform=transform, fill=float(b.ground_z.median()), dtype="float64")
raster_vis = np.array([visible_on_surface(surface, transform, (o.x, o.y, o.z), target) for o in obs.itertuples()])
print(f"mesh and raster agree on {(raster_vis == vis).mean():.1%} of {len(obs):,} observers; "
      f"raster says visible {raster_vis.sum():,}, mesh {vis.sum():,}")

mesh and raster agree on 88.4% of 2,436 observers; raster says visible 839, mesh 1,066

On the same 2,436 observers, a 1 m surface built from the footprints agreed with the mesh on 88.4%, and was the stricter of the two: 839 visible against 1,066. Vectorising one profile is already fast; for millions of pairs, batch the profiles into a 2D array and compare in one operation.

Example 3 — earth curvature and refraction for long views

import numpy as np

def apparent_drop(distance_m, k=0.13, R=6_371_000.0):
    """Drop of a distant target below the tangent plane, with standard refraction."""
    return (1 - k) * distance_m ** 2 / (2 * R)

for km in (1, 2, 5, 10, 20, 50):
    print(f"{km:3d} km: {apparent_drop(km * 1000):7.2f} m")

  1 km:    0.07 m
  2 km:    0.27 m
  5 km:    1.71 m
 10 km:    6.83 m
 20 km:   27.31 m
 50 km:  170.70 m

Below about 5 km the correction is smaller than the uncertainty in most building heights. Beyond 10 km it is the dominant term and omitting it makes distant views look possible when they are not.

Stepping off the facade, and target sets

Why the epsilon step-off matters so much

A ray starting exactly on a triangle intersects it at t = 0, and intersects_location duly reports a hit at zero distance. Every observer on a facade is then reported as blind. Stepping half a metre along the ray before casting removes the self-intersection without materially changing the geometry — and half a metre is also roughly the offset a window has from the facade plane.

Why targets should be sets of points

A landmark's visibility is not binary. A cathedral tower may have its spire visible and its body hidden; a hill may be visible over one roof and not another. Sampling a handful of points on the target's silhouette and reporting the fraction visible turns an argument about a definition into a number.

Why the raster method can disagree with the mesh

A DSM cannot represent an overhang, an arcade or a bridge, and it quantises building edges to the cell size. Near a boundary the two methods disagree by roughly the cell size projected along the ray, which for a 1 m DSM at a 200 m view is a fraction of a degree — irrelevant for a city-wide statistic and decisive for a single protected view.

Why trees are usually the answer nobody wants

Street trees block more urban views than buildings do, they change with the season, and no city model contains them reliably. A visibility analysis on buildings alone is an upper bound, and saying so is better than presenting it as a measurement.

Uneven windows, glass and curvature

• Step off the facade. Otherwise every observer is blind.
• Check the storey exists. Floor 8 of a four-storey building is not an observer.
• Windows are not evenly spaced. Facade sampling is a proxy.
• Glass and courtyards produce views that a solid model denies.
• Terrain first. In hilly cities it dominates the buildings.
• Batch the rays. intersects_location accepts arrays; per-ray Python loops are slow.
• Curvature matters beyond 5 km, and refraction varies with the weather.
• Vegetation is missing. Report the result as an upper bound.

Internal links

• How to calculate a viewshed in Python — the terrain-only version of this question
• How to compute sky view factor from a surface model — the same horizon geometry, integrated
• How to extrude building footprints into 3D in Python — producing the meshes
• How to read a CityJSON city model in Python — a model with real roofs
• 3D spatial data models explained: 2.5D, meshes and solids — what each model can block
• How to compute building shadows for a given time in Python — the same obstruction problem with a moving source
• Where building heights come from, and how wrong they are — the uncertainty in every ray
• How to extract raster values at points with rasterio — sampling the surface profile

FAQ

How do I test line of sight between two points in 3D?

Cast a ray from the observer to the target and check whether any geometry intersects it closer than the target. trimesh.ray.ray_triangle.RayMeshIntersector does this against a mesh.

Why does every observer report as blind?

Because the ray starts on its own facade and hits it at distance zero. Step half a metre along the ray before casting.

Should I use meshes or a surface model?

Meshes for a small number of exact tests, a DSM profile for millions of pairs. The two disagree near building edges by about the cell size.

How do I place window observers?

Walk the footprint boundary at a fixed spacing, offset outwards along the outward normal, and set z from the ground height plus the storey height times the floor number.

When does earth curvature matter?

Beyond about 5 km, where the drop exceeds 1.7 m. At 20 km it is 27 m and dominates the result.

Do I need to include trees?

Ideally yes, and you usually cannot. A buildings-only result is an upper bound on visibility and should be labelled as one.