3D Tiles and the web delivery of 3D city models explained

Problem statement

A city-scale 3D model is gigabytes. A browser has a few hundred megabytes of memory and a frame budget of sixteen milliseconds. Something has to decide which buildings to send, at what detail, in what order — and that something is a tiling scheme, not a file format.

3D Tiles is the OGC standard for this: a spatial tree of tiles, each with a bounding volume and a geometric error, and a client that walks the tree until the error is below what the current view can resolve. Understanding it is mostly understanding that one number, because everything about the experience — what loads, when, and how much — follows from it.

This guide covers the tileset structure, how geometric error drives loading, the payload formats, and how a Python-produced city model gets into one.

Quick answer

A tileset is a JSON tree; each node has a bounding volume, a geometric error and optionally a content URI:

{
  "asset": {"version": "1.1"},
  "geometricError": 500,
  "root": {
    "boundingVolume": {"region": [0.074496, 0.909346, 0.074669, 0.909442, 46.4, 80.9]},
    "geometricError": 100,
    "refine": "REPLACE",
    "content": {"uri": "l0/tile.glb"},
    "children": [
      {"boundingVolume": {"...": "..."}, "geometricError": 20,
       "content": {"uri": "l1/0.glb"}}
    ]
  }
}

The client computes the on-screen error of each tile from its geometricError and the camera distance, and descends only where that exceeds a threshold. A tileset whose geometric errors are wrong will either load everything at once or never load any detail.

Geometric error, refinement, bounding volumes, payload

1. Understand geometric error

It is the error, in metres, introduced by rendering this tile instead of its children. A root tile that replaces a whole city has a large value; a leaf holding real geometry has zero or near zero. The client projects it to screen space and compares against a Screen Space Error budget — CesiumJS's maximumScreenSpaceError defaults to 16 pixels.

import math

def refine_distance(geometric_error, sse=16, screen_px=1080, fov_deg=60):
    """Distance below which a tile's projected error exceeds the budget, so it refines."""
    return geometric_error * screen_px / (2 * sse * math.tan(math.radians(fov_deg) / 2))

for ge in (200, 50, 20, 1):
    print(f"geometricError {ge:>3} m refines closer than {refine_distance(ge):,.0f} m")

geometricError 200 m refines closer than 11,691 m
geometricError  50 m refines closer than 2,923 m
geometricError  20 m refines closer than 1,169 m
geometricError   1 m refines closer than 58 m

2. Choose REPLACE or ADD

REPLACE means a child's geometry replaces the parent's — used for level-of-detail pyramids. ADD means children add to the parent — used when the parent holds a coarse subset and children fill in. Most building tilesets use REPLACE.

3. Pick a bounding volume type

region is a geographic box in radians plus heights, box is an oriented box in the tileset's coordinate system, and sphere is a centre and radius. region is easiest to produce from geographic data; box gives tighter culling for a city block.

4. Know the payload

Tile content is glTF (.glb) in 3D Tiles 1.1, with the older .b3dm container still widely encountered. glTF is Y-up and uses a right-handed coordinate system, so a Z-up city model needs a rotation — (x, y, z) becomes (x, z, −y), as in Example 1 — or the city lies on its side.

5. Put the attributes somewhere the client can read

3D Tiles 1.1 uses EXT_mesh_features and EXT_structural_metadata to attach per-feature properties, so a click on a building can show its id, height and use. Without them the tileset is scenery.

6. Georeference the tileset

A tileset has a transform placing its local coordinates on the globe. The usual pattern is model coordinates in a local ENU frame plus a 4×4 transform to ECEF. Getting this wrong puts the city in the ocean or underground.

7. Generate rather than hand-write

py3dtiles produces tilesets from point clouds; conversion tools exist for CityJSON and CityGML. The examples below write the glTF and the JSON directly, which is short enough for one tile and shows every part a converter would otherwise hide. Hand-writing a tileset is reasonable for a single small tile and unreasonable beyond that.

8. Consider whether you need 3D Tiles at all

For a few thousand buildings in one neighbourhood, a single glTF or a GeoJSON extruded by the map library is simpler, smaller in total and far less work. 3D Tiles earns its complexity at city scale.

A minimal and a two-level tileset

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): meshes holds one trimesh per building solid, fan-triangulated from its CityJSON surfaces, with ids and ridge heights alongside. The coordinates are Dutch RD with NAP heights (EPSG:7415). Every tileset below is checked with Cesium's 3d-tiles-validator (npx 3d-tiles-validator --tilesetFile …).

import json
import os

import numpy as np
import requests
import trimesh

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

meshes, ids, heights = [], [], []
for oid, o in d["CityObjects"].items():
    a = o.get("attributes", {})
    for g in o.get("geometry", []):
        if g["type"] != "Solid" or "RelativeRidgeHeight" not in a:
            continue
        faces = [[ring[0], ring[k], ring[k + 1]]
                 for shell in g["boundaries"] for surface in shell
                 for ring in surface[:1] for k in range(1, len(ring) - 1)]
        m = trimesh.Trimesh(V, faces, process=False)
        m.remove_unreferenced_vertices()
        meshes.append(m)
        ids.append(oid)
        heights.append(a["RelativeRidgeHeight"])
print(f"{len(meshes):,} building solids, {sum(len(m.faces) for m in meshes):,} triangles; "
      f"centre {np.vstack([m.bounds for m in meshes]).mean(axis=0).round(0).tolist()}")

1,990 building solids, 35,602 triangles; centre [78693.0, 457892.0, 10.0]

Example 1 — a minimal tileset around a glTF you produced

import json, pathlib
import numpy as np, pyproj, trimesh

pyproj.network.set_network_enabled(True)          # NAP → ellipsoid needs the NLGEO2018 grid
to_geo = pyproj.Transformer.from_crs("EPSG:7415", "EPSG:4979", always_xy=True)   # RD+NAP → lon, lat, h
to_ecef = pyproj.Transformer.from_crs("EPSG:4979", "EPSG:4978", always_xy=True)

def enu_frame(lon, lat, h):
    """Rows are east, north, up at a point; plus the point in ECEF."""
    lo, la = np.radians(lon), np.radians(lat)
    R = np.array([[-np.sin(lo), np.cos(lo), 0],
                  [-np.sin(la) * np.cos(lo), -np.sin(la) * np.sin(lo), np.cos(la)],
                  [np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)]])
    return R, np.array(to_ecef.transform(lon, lat, h))

def write_tile(meshes, glb_path, frame=None):
    """Mesh in a local east-north-up frame, rotated to glTF's Y-up; returns region and frame."""
    scene = trimesh.util.concatenate(meshes)
    lon, lat, h = to_geo.transform(*scene.vertices.T)
    if frame is None:
        frame = enu_frame(lon.mean(), lat.mean(), h.min())
    R, origin = frame
    enu = (np.c_[to_ecef.transform(lon, lat, h)] - origin) @ R.T
    yup = np.c_[enu[:, 0], enu[:, 2], -enu[:, 1]]            # (x, y, z) → (x, z, −y)
    trimesh.Trimesh(yup, scene.faces, process=False).export(glb_path)
    region = [*np.radians([lon.min(), lat.min(), lon.max(), lat.max()]), float(h.min()), float(h.max())]
    return region, frame

def root_transform(frame):
    R, origin = frame
    T = np.eye(4)
    T[:3, :3], T[:3, 3] = R.T, origin                         # local ENU → ECEF
    return T.flatten(order="F").tolist()                      # 3D Tiles matrices are column-major

region, frame = write_tile(meshes, "block.glb")
tileset = {"asset": {"version": "1.1"}, "geometricError": 200,
           "root": {"transform": root_transform(frame),
                    "boundingVolume": {"region": region}, "geometricError": 50,
                    "refine": "REPLACE", "content": {"uri": "block.glb"}}}
pathlib.Path("tileset.json").write_text(json.dumps(tileset, indent=2))

nap = np.vstack([m.vertices for m in meshes])[:, 2]
print("region:", [round(float(v), 6) for v in region[:4]], [round(v, 1) for v in region[4:]])
print(f"heights: NAP {nap.min():.1f} to {nap.max():.1f} m → ellipsoidal {region[4]:.1f} to {region[5]:.1f} m")
print(f"block.glb {pathlib.Path('block.glb').stat().st_size / 1e6:.2f} MB")

region: [0.074496, 0.909346, 0.074669, 0.909442] [46.4, 80.9]
heights: NAP 3.0 to 37.5 m → ellipsoidal 46.4 to 80.9 m
block.glb 0.69 MB

import subprocess

def validate(path):
    out = subprocess.run(["npx", "--yes", "3d-tiles-validator", "--tilesetFile", path],
                         capture_output=True, text=True).stdout
    report = json.loads(out[out.index("{"):])
    return report["numErrors"], report["numWarnings"], [i["type"] for i in report.get("issues", [])]

print(validate("tileset.json"))

(0, 0, ['CONTENT_VALIDATION_INFO'])

The heights in a region bounding volume are metres above the WGS84 ellipsoid, not above a national datum: this model's NAP heights become ellipsoidal ones about 43 m higher. The transform goes through ETRS89, which differs from WGS84 by under a metre in the Netherlands today — well inside any tile's geometric error.

Example 2 — a two-level tileset from a height threshold

import numpy as np, json, pathlib

def two_level_tileset(meshes, heights, tall_threshold=15, chunk=500):
    """Tall buildings in the root, everything else in children; one local frame for all."""
    heights = np.asarray(heights)
    tall = [m for m, h in zip(meshes, heights) if h >= tall_threshold]
    rest = [m for m, h in zip(meshes, heights) if h < tall_threshold]

    root_region, frame = write_tile(meshes, "all.glb")            # the frame and region for the whole set
    tall_region, _ = write_tile(tall, "tall.glb", frame)
    children = []
    for i in range(0, len(rest), chunk):
        region, _ = write_tile(rest[i:i + chunk], f"c{i // chunk}.glb", frame)
        children.append({"boundingVolume": {"region": region},
                         "geometricError": 0, "content": {"uri": f"c{i // chunk}.glb"}})
    return {"asset": {"version": "1.1"}, "geometricError": 200,
            "root": {"transform": root_transform(frame),
                     "boundingVolume": {"region": root_region},
                     "geometricError": 20, "refine": "ADD",
                     "content": {"uri": "tall.glb"}, "children": children}}

ts = two_level_tileset(meshes, heights)
pathlib.Path("tileset2.json").write_text(json.dumps(ts, indent=2))
print(f"{sum(h >= 15 for h in heights)} buildings of 15 m or more in the root, "
      f"{len(ts['root']['children'])} child tiles; validator:", validate("tileset2.json"))

boxes = [c["boundingVolume"]["region"] for c in ts["root"]["children"]]
overlap = lambda a, b: a[0] < b[2] and b[0] < a[2] and a[1] < b[3] and b[1] < a[3]
pairs = [(i, j) for i in range(len(boxes)) for j in range(i + 1, len(boxes)) if overlap(boxes[i], boxes[j])]
print(f"{len(pairs)} of {len(boxes) * (len(boxes) - 1) // 2} child pairs have overlapping regions")

26 buildings of 15 m or more in the root, 4 child tiles; validator: (0, 0, ['CONTENT_VALIDATION_INFO', 'CONTENT_VALIDATION_INFO'])
6 of 6 child pairs have overlapping regions

ADD refinement here means the tall buildings stay visible as you zoom in and the smaller ones appear alongside them, which is usually what a city view wants. The children are cut in file order, 500 solids at a time, and on this tile every one of the six pairs of child regions overlaps — valid, but a client cannot cull one without the others. Cutting by space (a quadtree on the footprints) is what makes the tree useful for culling.

Example 3 — check the tileset before serving it

import json, math

def audit_tileset(path):
    ts = json.loads(open(path).read())
    problems = []

    def walk(node, parent_error=math.inf, depth=0):
        ge = node.get("geometricError")
        if ge is None:
            problems.append(f"depth {depth}: no geometricError")
        elif ge >= parent_error:
            problems.append(f"depth {depth}: geometricError {ge} not less than parent {parent_error}")
        bv = node.get("boundingVolume", {})
        if "region" in bv:
            w, s, e, n, lo, hi = bv["region"]
            if w >= e or s >= n:
                problems.append(f"depth {depth}: region is not west/south/east/north")
            if abs(w) > math.pi or abs(s) > math.pi / 2:
                problems.append(f"depth {depth}: region looks like degrees, not radians")
            if hi < lo:
                problems.append(f"depth {depth}: heights reversed")
        for child in node.get("children", []):
            walk(child, ge if ge is not None else parent_error, depth + 1)

    walk(ts["root"], ts.get("geometricError", math.inf))
    return problems

import copy

broken = copy.deepcopy(json.loads(open("tileset.json").read()))
w, so, e, n, lo, hi = broken["root"]["boundingVolume"]["region"]
broken["root"]["boundingVolume"]["region"] = [*np.degrees([w, so, e, n]), lo - 43.4, hi - 43.4]  # degrees, NAP
broken["root"]["geometricError"] = 250                                                      # above the tileset's
pathlib.Path("broken.json").write_text(json.dumps(broken))

for path in ("tileset.json", "tileset2.json", "broken.json"):
    print(path, "audit:", audit_tileset(path) or "clean", "| validator:", validate(path))

tileset.json audit: clean | validator: (0, 0, ['CONTENT_VALIDATION_INFO'])
tileset2.json audit: clean | validator: (0, 0, ['CONTENT_VALIDATION_INFO', 'CONTENT_VALIDATION_INFO'])
broken.json audit: ['depth 0: geometricError 250 not less than parent 200', 'depth 0: region looks like degrees, not radians'] | validator: (4, 0, ['BOUNDING_VOLUME_INVALID', 'BOUNDING_VOLUME_INVALID', 'BOUNDING_VOLUME_INVALID', 'BOUNDING_VOLUME_INVALID'])

Degrees in a region is the classic tileset error, because the file is still well-formed JSON. Both checks catch it — Cesium's validator with four BOUNDING_VOLUME_INVALID errors. Neither catches the other two faults in broken.json the same way: the audit flags a root error larger than the tileset's, which the validator lets through, and nothing notices that the heights were left in NAP, 43.4 m too low. A datum is a fact about the data, not about the file.

Geometric error is the whole design

Why geometric error is the whole design

Everything a 3D Tiles client does follows from comparing a tile's geometric error, scaled by distance and field of view, against a pixel budget. A tileset whose errors decrease properly towards the leaves loads progressively and looks right at every zoom. Errors that are too large refine too early — at 200 m, a tile's children are wanted from 11.7 km away (step 1), so the whole tree loads at once; errors that are too small never trigger refinement, and the view stays at the coarse level however close the camera comes.

Why glTF is Y-up and why that matters

glTF was designed for real-time graphics, where Y-up right-handed is the convention. Geospatial data is Z-up. The conversion is a −90° rotation about X, and forgetting it produces a model that is not subtly wrong but lying on its side — which at least makes it an easy bug to spot.

Why heights are ellipsoidal in a region volume

The region bounding volume is defined in the WGS84 geographic frame with heights above the ellipsoid, because that is what the globe renderer uses. National vertical datums are offset from the ellipsoid by tens of metres — 43.4 m under this tile, and 38–53 m across most of western and southern Europe (floating or sunken buildings measures them) — so a model placed with orthometric heights sinks by that amount.

Why per-feature metadata is what makes it useful

A tileset without feature metadata is a picture of a city. With EXT_mesh_features and EXT_structural_metadata, each building carries an id and properties, so a click selects a feature, a style can colour by height or use, and the viewer becomes a map rather than scenery.

Radians, legacy containers and Draco

• region is in radians. Degrees is the commonest mistake.
• .b3dm is the legacy container. Still common; .glb is the current one.
• Draco compression is worth it for large meshes and costs decode time.
• Refinement mode is per node, not per tileset.
• External tilesets nest. A content.uri may point at another tileset.json.
• Screen space error is a client setting. Do not tune your errors around one client's default.
• Terrain and buildings are separate tilesets and must share a datum.
• Cesium and its derivatives are the main clients. Test in the one your users have.

Internal links

• How to publish a 3D building layer to a web map — the practical route
• How to export a 3D city model to glTF in Python — producing the payload
• A 3D viewer shows nothing or a black screen — when the tileset loads and nothing appears
• CityJSON and CityGML explained — the source model
• Vertical datums explained — the 43 m offset
• Web map tiles explained — the 2D equivalent of this idea
• Vector tiles explained — where extruded GeoJSON fits
• How to serve PMTiles with MapLibre — the simpler delivery route

FAQ

What is 3D Tiles?

An OGC standard for streaming large 3D geospatial datasets: a spatial tree of tiles, each with a bounding volume and a geometric error, whose payloads are glTF meshes or point clouds.

What does geometric error mean?

The error in metres introduced by rendering a tile instead of its children. The client projects it to screen space and refines only where it exceeds a pixel budget.

Do I need 3D Tiles for a few thousand buildings?

Probably not. A single glTF or a GeoJSON extruded by the map library is simpler and smaller in total. 3D Tiles pays off at city scale.

Why is my tileset lying on its side?

glTF is Y-up and geospatial data is Z-up. The payload needs a −90° rotation about the X axis.

Why is my city underground?

The region bounding volume uses heights above the WGS84 ellipsoid. A model in a national datum needs the geoid separation added — about 43 m in the Netherlands.

How do I make buildings clickable?

Attach per-feature metadata with EXT_mesh_features and EXT_structural_metadata. Without them the tileset has no features to select.