How to export a 3D city model to glTF in Python

On this page

Problem statement

glTF is the format every 3D viewer reads: browsers, game engines, Blender, 3D Tiles payloads. Getting a city model into it is straightforward with trimesh, and three things go wrong every time.

The first is orientation: glTF is Y-up, geospatial data is Z-up, so a model exported without a rotation lies on its side. The second is precision: national grid coordinates are in the hundreds of thousands of metres, and glTF stores vertex positions as 32-bit floats, which at 500,000 have a spacing of about 3 cm β€” enough to make a building visibly ragged. The third is that a mesh has no attributes, so a model exported without an id mapping cannot be joined back to anything.

This guide exports a real city model with all three handled.

Quick answer

import numpy as np, trimesh

def to_gltf(meshes, out="model.glb", origin=None):
    scene = trimesh.util.concatenate(meshes)
    if origin is None:
        origin = scene.bounds.mean(axis=0)
    scene.apply_translation(-origin)                     # kill the float32 precision problem
    scene.apply_transform(trimesh.transformations.rotation_matrix(
        -np.pi / 2, [1, 0, 0]))                          # Z-up β†’ Y-up
    scene.export(out)
    return out, origin

path, origin = to_gltf(meshes)
print(path, origin.round(2).tolist())
model.glb [78681.73, 457939.06, 20.24]

Recentring is not optional. The Hague tile on the Dutch national grid is centred near x = 78,693, y = 457,892; at that northing float32 has a representable spacing of 3.1 cm (Example 2), and every vertex snaps to it. Recentred, the exported model's largest coordinate is 355 m.

Three glTF export pitfalls and the fix for each. three ways to export a wrong filepitfallfixZ-up data in a Y-up formatrotate βˆ’90Β° about Xnational grid coordinates infloat32subtract a local originno attributes on the meshnode names plus a sidecar id mapAll three produce a file that opens and is wrong in a different way spatialworkflow.io
All three produce a file that opens and looks wrong in a different way.

Meshes, grouping, a local origin, Y-up, materials

1. Build the meshes

Either from extruded footprints or from a CityJSON model's solids. Keep them as a list, one per building, until you have decided how to group them.

2. Decide the grouping

One mesh per building lets a viewer select buildings and costs a draw call each. One merged mesh draws fast and is a single object. The usual compromise is a mesh per block or per tile of a few hundred buildings.

3. Recentre on a local origin

Subtract a fixed origin β€” the model's centre, or a round number near it β€” and record it. Everything downstream needs to know it, including the 3D Tiles transform that puts the model back on the globe.

4. Rotate Z-up to Y-up

A βˆ’90Β° rotation about X. trimesh.transformations.rotation_matrix(-np.pi/2, [1,0,0]) is the whole fix.

5. Attach colours or materials

glTF carries vertex colours and PBR materials. Colouring by height, use or age is the cheapest way to make a model readable, and vertex colours survive every exporter.

import matplotlib
norm = matplotlib.colors.Normalize(vmin=0, vmax=30)
cmap = matplotlib.colormaps["viridis"]
for m, h in zip(meshes, heights):
    m.visual.vertex_colors = (np.array(cmap(norm(h))) * 255).astype("uint8")

6. Keep the id mapping

Write a JSON beside the .glb mapping mesh index or node name to the building id, so a selection in the viewer can be resolved back to the data.

7. Check the file before shipping it

Reload it, compare the bounds against the source, count the meshes, and look at it. A model that is 90Β° out is obvious in a viewer and invisible in a byte count.

A Z-up model and the same model after the rotation to Y-up. Z-up (as exported)ZXlies on its side in a viewerY-up (after rotation)YXrotation_matrix(βˆ’Ο€/2, [1,0,0])One rotation matrix; nothing else about the model is wrong spatialworkflow.io
One rotation matrix; without it the city lies on its side and nothing else is wrong.

Exporting with names, colours and an id map

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), in the Dutch national grid with NAP heights (EPSG:7415). The setup turns every solid that carries the city's recorded heights into a trimesh mesh β€” each surface's outer ring fanned into triangles β€” and keeps its CityJSON id and its ridge height above ground as meshes, ids and heights.

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 β€” export with names, colours and an id map

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

def export_city(meshes, ids, heights, out="city.glb", origin=None, vmax=30):
    scene = trimesh.Scene()
    cmap = matplotlib.colormaps["viridis"]
    norm = matplotlib.colors.Normalize(vmin=0, vmax=vmax)

    all_bounds = np.vstack([m.bounds for m in meshes])
    if origin is None:
        origin = np.round(all_bounds.mean(axis=0), 0)      # a round number near the centre

    rot = trimesh.transformations.rotation_matrix(-np.pi / 2, [1, 0, 0])
    for m, bid, h in zip(meshes, ids, heights):
        g = m.copy()
        g.apply_translation(-origin)
        g.apply_transform(rot)
        g.visual.vertex_colors = (np.array(cmap(norm(h))) * 255).astype("uint8")
        scene.add_geometry(g, node_name=str(bid), geom_name=str(bid))

    scene.export(out)
    pathlib.Path(out).with_suffix(".ids.json").write_text(json.dumps({
        "origin": origin.tolist(),
        "crs": "EPSG:7415",
        "axis_convention": "Y-up (rotated -90Β° about X from Z-up)",
        "nodes": {str(b): str(b) for b in ids},
    }, indent=2))
    return out, origin


out, city_origin = export_city(meshes, ids, heights)
print(out, city_origin.tolist(), f"{os.path.getsize(out) / 1e6:.2f} MB")
city.glb [78693.0, 457892.0, 10.0] 1.47 MB

The sidecar is what makes the export reversible. Without the origin and the CRS, a .glb is a shape with no location.

Example 2 β€” prove the precision problem to yourself

import numpy as np

for x in (1e3, 1e5, 4.579e5, 1e6):
    spacing = np.spacing(np.float32(x))
    print(f"x = {x:>9,.0f} m: float32 spacing {spacing:.6f} m")
x =     1,000 m: float32 spacing 0.000061 m
x =   100,000 m: float32 spacing 0.007812 m
x =   457,900 m: float32 spacing 0.031250 m
x = 1,000,000 m: float32 spacing 0.062500 m

At national-grid values the quantisation is centimetres. Recentred on a local origin the same model has sub-millimetre precision.

Example 3 β€” check the export

import numpy as np, trimesh

def check_gltf(path, source_meshes, origin):
    loaded = trimesh.load(path)
    meshes = list(loaded.geometry.values()) if hasattr(loaded, "geometry") else [loaded]

    src_bounds = np.vstack([m.bounds for m in source_meshes])
    src_size = src_bounds.max(axis=0) - src_bounds.min(axis=0)      # x, y, z extent
    out_size = loaded.bounds[1] - loaded.bounds[0]                  # x, z, y after rotation

    report = {
        "meshes": len(meshes),
        "expected_meshes": len(source_meshes),
        "source_extent_xyz_m": src_size.round(1).tolist(),
        "exported_extent_m": out_size.round(1).tolist(),
        "vertical_axis_looks_like": "Y" if out_size[1] < out_size[0] else "check this",
        "max_abs_coordinate": float(np.abs(loaded.bounds).max()),
    }
    if report["max_abs_coordinate"] > 1e4:
        report["warning"] = "coordinates are large β€” was the model recentred?"
    return report


for key, value in check_gltf("city.glb", meshes, city_origin).items():
    print(f"{key:26} {value}")
meshes                     1990
expected_meshes            1990
source_extent_xyz_m        [677.8, 615.6, 34.5]
exported_extent_m          [677.8, 34.5, 615.6]
vertical_axis_looks_like   Y
max_abs_coordinate         354.85699462890625

Comparing extents rather than coordinates is the trick: after a correct rotation the x extent is unchanged, the old z extent becomes the y extent, and the old y extent becomes z.

Y-up, float32 and attributes outside the mesh

Why glTF is Y-up

glTF came from real-time graphics, where the convention is a right-handed Y-up coordinate system because it matches how a camera looks at a scene. Geospatial data is Z-up because z is height. Neither is wrong, and the conversion is a fixed rotation β€” but nothing in either format records which convention a file uses, so the only signal is that the model looks wrong.

Why float32 is enough for graphics and not for geodesy

Graphics coordinates are usually within a few thousand units of the origin, where float32 has millimetre precision. A national grid coordinate is five or six digits before the decimal point, which consumes most of the mantissa. The fix is the same one used everywhere in 3D geospatial work: store geometry in a local frame and carry the transform separately, which is exactly what a 3D Tiles tileset transform does.

Why meshes cannot carry your attributes

glTF has no concept of feature attributes in its core. The 3D Tiles extensions EXT_mesh_features and EXT_structural_metadata add them, and most writers do not emit them. Until you use those, a node name plus a sidecar map is how a click in the viewer becomes a building id.

Why one mesh per building is usually right at city block scale

Draw calls cost, so a thousand separate meshes is slower than one merged mesh. But a merged mesh cannot be selected, styled or hidden per building, which is most of the point of a city model. A few hundred buildings per mesh, grouped spatially, keeps both properties acceptable.

float32 representable spacing at four coordinate magnitudes. float32 spacing in metres1,000 m0.06 mm100,000 m7.8 mm457,900 m31 mm1,000,000 m63 mm00.07Recentred on a local origin, the same model has sub-millimetre precision spatialworkflow.io
At national-grid magnitudes the quantisation is centimetres, which is visible.

Binary containers, Draco and vertex colours

  • .glb is the binary container. Prefer it to .gltf plus a .bin.
  • Draco compression shrinks large meshes considerably and needs decoder support.
  • Vertex colours are per vertex. For flat per-face colour, split the vertices.
  • Normals are computed if absent and are usually fine for buildings.
  • Units are metres by convention in glTF; there is no unit field.
  • Textures multiply the file size. A city model rarely needs them.
  • Node names must be unique or viewers will silently merge them.
  • Record the origin and CRS. A .glb on its own has no location.

FAQ

How do I export a 3D model to glTF in Python?

Build trimesh meshes, concatenate or add them to a Scene, recentre on a local origin, rotate βˆ’90Β° about X, and call .export("model.glb").

Why is my exported model lying on its side?

glTF is Y-up and geospatial data is Z-up. Apply rotation_matrix(-np.pi/2, [1,0,0]) before exporting.

Why does my model look ragged or faceted?

Float32 precision. At a northing of 457,900 m the representable spacing is 3.1 cm. Subtract a local origin before exporting and record it.

Can I put attributes in a glTF?

Not in the core format. Use node names plus a sidecar id map, or the 3D Tiles metadata extensions if your viewer supports them.

Should I merge all the buildings into one mesh?

Merge for speed, keep separate for selectability. A few hundred buildings per mesh, grouped spatially, is the usual compromise.

How do I check the export is right?

Reload it and compare extents rather than coordinates: after a correct rotation the x extent is unchanged and the old z extent becomes the y extent.