A building mesh is not watertight and its volume is wrong
On this page
Problem statement
mesh.volume returns a number. It always returns a number β the calculation is a sum over faces, and it does not require the faces to enclose anything. On an open mesh the result is plausible, wrong, and accompanied by no warning at all.
On a real LoD2 tile, 1,985 of 1,990 solids triangulated to watertight meshes and five did not. Those five reported volumes like everything else. The 1,985 watertight solids hold 683,065 mΒ³; with the five filled in as prisms to the eaves the tile totals 683,115 mΒ³ (Example 3). Here the difference is small, and the point is that nothing in the numbers says which five to distrust.
This guide diagnoses why a mesh leaks and fixes the causes that are fixable.
Quick answer
Assert, then diagnose:
import trimesh
assert mesh.is_watertight, (
f"not watertight: {len(mesh.faces)} faces, "
f"{len(trimesh.repair.broken_faces(mesh))} broken, "
f"Euler number {mesh.euler_number}")
print(f"volume {abs(mesh.volume):,.1f} mΒ³, area {mesh.area:,.1f} mΒ²")
euler_number is the quickest diagnostic. A closed surface of genus 0 β a building without a hole through it β has Euler number 2. Anything else means a hole, a duplicated shell or a non-manifold edge.
Diagnose the break, repair, check the triangulation
1. Find out what kind of broken it is
import trimesh
def diagnose(mesh):
return {
"watertight": bool(mesh.is_watertight),
"winding_consistent": bool(mesh.is_winding_consistent),
"volume_sign": float(mesh.volume),
"euler_number": int(mesh.euler_number),
"faces": len(mesh.faces),
"broken_faces": len(trimesh.repair.broken_faces(mesh)),
"duplicate_faces": len(mesh.faces) - len(trimesh.grouping.unique_rows(
np.sort(mesh.faces, axis=1))[0]),
"unreferenced_vertices": len(mesh.vertices) - len(np.unique(mesh.faces)),
}
print(diagnose(mesh))
{'watertight': False, 'winding_consistent': True, 'volume_sign': 11.855092693275461, 'euler_number': 3, 'faces': 8, 'broken_faces': 6, 'duplicate_faces': 0, 'unreferenced_vertices': 0}
2. Try the automatic repairs first
Two problems fix themselves: inconsistent winding and duplicated or degenerate faces.
import trimesh
mesh.update_faces(mesh.unique_faces()) # trimesh 4: was remove_duplicate_faces()
mesh.update_faces(mesh.nondegenerate_faces()) # trimesh 4: was remove_degenerate_faces()
mesh.remove_unreferenced_vertices()
trimesh.repair.fix_normals(mesh)
trimesh.repair.fix_winding(mesh)
trimesh.repair.fill_holes(mesh)
print("watertight after the automatic repairs:", mesh.is_watertight)
watertight after the automatic repairs: True
fill_holes only closes simple boundary loops. It will not invent a roof.
3. Check the triangulation if the mesh came from a city model
Fan triangulation β splitting a ring from its first vertex β is correct for convex rings and can go wrong for concave ones. On the reference tile it was not the cause: 478 solids with a concave ring came out watertight (Example 4). The five failures each had two rings that repeat a vertex index, one closed explicitly by repeating its first vertex and one with the same vertex twice in a row β 10 of the tile's 16,209 rings. Fanning those produces zero-area triangles, and the degenerate faces break the manifold test. Dropping degenerate faces (step 2) fixes all five, and so does ear-clipping each surface in its own plane:
import numpy as np, mapbox_earcut, shapely
def triangulate_planar(points):
"""Ear-clip a planar 3D ring after projecting to its own plane."""
n = np.zeros(3)
for i in range(len(points)):
n += np.cross(points[i], points[(i + 1) % len(points)])
n /= np.linalg.norm(n)
u = np.cross(n, [0, 0, 1] if abs(n[2]) < 0.9 else [1, 0, 0])
u /= np.linalg.norm(u)
v = np.cross(n, u)
flat = np.c_[points @ u, points @ v]
idx = mapbox_earcut.triangulate_float64(flat, np.array([len(flat)]))
return np.asarray(idx).reshape(-1, 3)
4. Check for a missing cap
An extrusion built by hand from wall quads has no top or bottom. mesh.euler_number will be 0 rather than 2, and fill_holes usually closes it if the boundary is a simple loop.
5. Check for coincident but unmerged vertices
Two faces that share an edge geometrically but use different vertex indices leave a crack. mesh.merge_vertices() closes it; process=True on construction does the same.
6. Check for a non-manifold edge
An edge used by more than two faces β where two building parts share a wall, or a mesh has an internal partition β is not repairable by welding. It needs the duplicate geometry removed.
7. If it cannot be repaired, report it rather than measuring it
A volume from an open mesh is not a conservative estimate or an approximation. Exclude the mesh, count the exclusions, and fall back to a prism volume where you have a footprint and a height.
Auditing and repairing a whole model
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, and sets geometry to the first solid's geometry and mesh to the first mesh that is not watertight.
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"],
"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()
print(f"{len(meshes):,} meshes, {sum(m.is_watertight for m in meshes):,} watertight")
1,990 meshes, 1,985 watertight
Example 1 β audit a whole model
import numpy as np, trimesh, pandas as pd
def audit_meshes(meshes, ids):
rows = []
for m, i in zip(meshes, ids):
rows.append({
"id": i,
"watertight": bool(m.is_watertight),
"winding_ok": bool(m.is_winding_consistent),
"euler": int(m.euler_number),
"faces": len(m.faces),
"volume_m3": float(abs(m.volume)),
})
df = pd.DataFrame(rows)
print(f"{len(df):,} meshes, {df.watertight.sum():,} watertight "
f"({df.watertight.mean():.1%})")
print(df[~df.watertight][["id", "euler", "faces"]].head(10).to_string(index=False))
return df
audit = audit_meshes(meshes, ids)
1,990 meshes, 1,985 watertight (99.7%)
id euler faces
GUID_233F8BDE-9657-4981-AE9F-85EEC9F0385C_19 3 8
GUID_233F8BDE-9657-4981-AE9F-85EEC9F0385C_18 3 8
GUID_C78A3205-27CC-4AD7-B240-8D2CFA89DB49_7 3 8
GUID_351BE36C-81E2-4C22-A7DD-642A0DE9D72D_16 3 8
GUID_351BE36C-81E2-4C22-A7DD-642A0DE9D72D_13 3 8
Example 2 β repair, then re-check
import trimesh
def repair(mesh):
m = mesh.copy()
m.merge_vertices()
m.update_faces(m.nondegenerate_faces())
m.update_faces(m.unique_faces())
m.remove_unreferenced_vertices()
trimesh.repair.fix_winding(m)
trimesh.repair.fix_inversion(m)
if not m.is_watertight:
trimesh.repair.fill_holes(m)
return m
fixed = [repair(m) for m in meshes]
before = sum(m.is_watertight for m in meshes)
after = sum(m.is_watertight for m in fixed)
print(f"watertight before {before:,}, after {after:,}")
watertight before 1,985, after 1,990
Report both numbers. Here the repair closed all five; one that fixes three of five is still a result; a repair that silently changes the geometry of the 1,985 that were already fine is a problem.
Example 3 β fall back to a prism, with the fallback recorded
import numpy as np, pandas as pd
def volumes_with_fallback(meshes, ids, footprints):
rows = []
for m, i in zip(meshes, ids):
if m.is_watertight:
rows.append({"id": i, "volume_m3": abs(m.volume), "method": "solid"})
else:
f = footprints.loc[i]
rows.append({"id": i, "volume_m3": f.geometry.area * f.height_eaves_m,
"method": "prism to eaves (mesh not watertight)"})
df = pd.DataFrame(rows)
print(df.method.value_counts().to_string())
return df
vols = volumes_with_fallback(meshes, ids, footprints)
print(f"total {vols.volume_m3.sum():,.0f} mΒ³")
method
solid 1985
prism to eaves (mesh not watertight) 5
total 683,115 mΒ³
Example 4 β is fan triangulation the cause?
from shapely.geometry import Polygon as P2
def ring_is_concave(points):
"""Concave if the ring, flattened into its own plane, is not its own convex hull."""
n = np.cross(points - points.mean(axis=0), np.roll(points, -1, axis=0) - points.mean(axis=0)).sum(axis=0)
n /= np.linalg.norm(n)
u = np.cross(n, [0, 0, 1] if abs(n[2]) < 0.9 else [1, 0, 0]); u /= np.linalg.norm(u)
flat = P2(np.c_[points @ u, points @ np.cross(n, u)])
return flat.convex_hull.area - flat.area > 1e-6
solids = {oid: g for oid, o in d["CityObjects"].items() for g in o.get("geometry", []) if g["type"] == "Solid"}
broken = [i for m, i in zip(meshes, ids) if not m.is_watertight]
for oid in broken:
g = solids[oid]
rings = [np.asarray(V[s[0]]) for shell in g["boundaries"] for s in shell]
tris = []
for shell in g["boundaries"]:
for s in shell:
ring = s[0]
local = triangulate_planar(V[ring])
tris += [[ring[a], ring[b], ring[c]] for a, b, c in local]
ear = trimesh.Trimesh(vertices=V, faces=np.array(tris), process=True)
print(f"{oid}: concave rings {sum(ring_is_concave(r) for r in rings)}, "
f"ear-clipped watertight {ear.is_watertight}, euler {ear.euler_number}")
concave_ok = sum(any(ring_is_concave(np.asarray(V[s[0]])) for shell in solids[i]["boundaries"] for s in shell)
for m, i in zip(meshes, ids) if m.is_watertight)
print(f"watertight solids that also have a concave ring: {concave_ok:,}")
GUID_233F8BDE-9657-4981-AE9F-85EEC9F0385C_19: concave rings 0, ear-clipped watertight True, euler 2
GUID_233F8BDE-9657-4981-AE9F-85EEC9F0385C_18: concave rings 0, ear-clipped watertight True, euler 2
GUID_C78A3205-27CC-4AD7-B240-8D2CFA89DB49_7: concave rings 0, ear-clipped watertight True, euler 2
GUID_351BE36C-81E2-4C22-A7DD-642A0DE9D72D_16: concave rings 0, ear-clipped watertight True, euler 2
GUID_351BE36C-81E2-4C22-A7DD-642A0DE9D72D_13: concave rings 0, ear-clipped watertight True, euler 2
watertight solids that also have a concave ring: 478
The method column is what stops the fallback from being invisible in the totals.
Why an open mesh still reports a volume
Why an open mesh still returns a volume
Trimesh computes volume as a sum of signed tetrahedron volumes formed by each triangle and the origin. For a closed, consistently oriented surface, the contributions outside the solid cancel exactly and the sum is the enclosed volume. With a hole they do not cancel, and the sum is some number related to the faces that happen to be present. It is not an underestimate, an approximation or a bound β it is arithmetic on an incomplete surface.
Why a repeated vertex produces the failures
A CityJSON ring must not repeat a vertex: the ring is closed implicitly. A ring that repeats one β [a, b, c, a], or [a, b, b, c] β is really a triangle stored as a quadrilateral, and fanning it from its first vertex produces one real triangle and one of zero area. The mesh looks right in a viewer, and the degenerate face breaks the edge bookkeeping that the manifold test relies on. That is the five-in-1,990 pattern observed here. Concave rings are the other classic culprit of fan triangulation β triangles that stick out of an L-shaped face β and on this tile they did no harm.
Why merging vertices is usually the first fix
Meshes assembled from independently triangulated faces often have coincident vertices with different indices, so adjacent faces do not share an edge. Every such pair is a crack. merge_vertices welds them within a tolerance, and it is by far the commonest single-line fix.
Why a non-manifold edge is a modelling problem, not a mesh problem
An edge shared by three or more faces means the geometry describes something that is not a surface β typically two building parts that share a wall, each modelled as a complete solid. The fix is in the model: either merge the parts into one solid, or keep them separate and do not concatenate them before measuring.
Euler numbers, winding and hole filling
euler_number2 is closed and genus 0. An arcade through a building gives 0.volumesign follows winding. Take the absolute value after checking consistency.process=Truemerges vertices on construction and is usually what you want.fill_holesonly closes simple loops. It cannot rebuild a missing roof.- Concatenating meshes creates non-manifold edges where they touch.
- Tolerances matter. Vertices a micrometre apart are coincident for a building.
- Report the failures. Five of 1,990 is a quality figure worth publishing.
- Prism fallback needs a height definition. Record which one.
Internal links
- How to calculate building volumes and floor area in Python β the calculation this protects
- How to read a CityJSON city model in Python β where the triangulation happens
- 3D spatial data models explained: 2.5D, meshes and solids β why watertightness is the defining property
- How to extrude building footprints into 3D in Python β producing watertight meshes by construction
- How to export a 3D city model to glTF in Python β where a broken mesh shows up as a hole
- How to fix invalid geometries in GeoPandas β the 2D equivalent
- Geometry validity explained β the same idea one dimension down
- A 3D viewer shows nothing or a black screen β when the normals are the problem
FAQ
Why does my mesh volume look wrong?
Because the mesh is probably not watertight. Trimesh returns a number regardless; check is_watertight before using it.
What does the Euler number tell me?
A closed surface with no holes through it has Euler number 2. Anything else indicates a boundary, a duplicate shell or a non-manifold edge.
How do I repair a mesh?
Merge vertices, remove degenerate and duplicate faces, fix winding, then fill simple holes. Re-check afterwards and report how many were fixed.
Why do some CityJSON solids fail?
On the reference tile, every failure came from a ring that repeats a vertex index, which fan triangulation turns into a zero-area triangle. Dropping degenerate faces fixes it; so does ear-clipping each face in its own plane, which also handles concave faces.
Can I just take the absolute value and carry on?
No. The sign is a winding problem; the magnitude is wrong for an open mesh whatever the sign.
What should I do with a mesh I cannot repair?
Exclude it from solid volumes, fall back to a prism from its footprint and height, and record which method each building used.