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.

Five reasons a building mesh is not watertight and the fix for each. why it leakscausefixcoincident but unmerged verticesmesh.merge_vertices()inconsistent windingtrimesh.repair.fix_windingmissing top or bottom capfill_holes, or extrude properlya ring that repeats a vertexdrop degenerate faces, or ear-clipnon-manifold shared wallremove the duplicate geometryThree repair automatically; two need the geometry fixed spatialworkflow.io
Five causes; three are repairable automatically and two need the geometry fixed.

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.

A closed prism, an open-topped one and a mesh with an unmerged crack. closed β€” Euler 2volume is meaningfulno cap β€” Euler 0volume is arithmeticcrackmerge_vertices closes it1,985 of 1,990 solids in the reference tile were watertight spatialworkflow.io
The Euler number distinguishes the cases before you look at the geometry.

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.

What to do with a mesh that cannot be repaired. when repair fails1assert is_watertightnever warn2try the automatic repairsreport before and after3exclude what still leakscount it4fall back to a prismrecord the method columnAn open mesh’s volume is not a bound, an estimate or an approximation spatialworkflow.io
The method column is what stops the fallback from vanishing into the totals.

Euler numbers, winding and hole filling

  • euler_number 2 is closed and genus 0. An arcade through a building gives 0.
  • volume sign follows winding. Take the absolute value after checking consistency.
  • process=True merges vertices on construction and is usually what you want.
  • fill_holes only 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.

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.