⚙️ How-To 3D & buildings Terrain & elevation 3D & Buildings · step 14 of 20 Updated

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.

Scene showing an observer at a window, a ray to a landmark and an intervening building blocking it.
The test is one ray; the work is placing the observer correctly.

Step-by-step solution

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.

def facade_points(poly, spacing=3.0, offset=0.5):
    line = poly.exterior
    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 2.74 m storeys puts the eye at roughly 5.96 + 3 × 2.74 + 1.5 ≈ 15.7 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.

Comparison grid of mesh ray casting and DSM profile sampling across exactness, speed and what each can represent.
Two methods, and the choice is decided by how many rays you need.

Code examples

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=2.74, 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():,}")

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]))

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.30 m
 50 km:  170.60 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.

Explanation

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.

Bars of the apparent drop from earth curvature with refraction at 1, 5, 10, 20 and 50 kilometres, rising from 0.07 m to 170.6 m.
Beyond 10 km curvature is the dominant term in a long-view test.

Edge cases or notes

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

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.