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

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

Recentring is not optional. A model on the Dutch national grid sits around x = 78,600, y = 457,900; in float32 that is a representable spacing of 3.1 cm, and every vertex snaps to it.

Triage of three glTF export pitfalls โ€” Y-up, float32 precision and lost attributes โ€” with the fix for each.
All three produce a file that opens and looks wrong in a different way.

Step-by-step solution

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.

Two scenes showing a Z-up model and the same model after the minus ninety degree rotation about X to Y-up.
One rotation matrix; without it the city lies on its side and nothing else is wrong.

Code examples

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.array([all_bounds[:, 0].min(), all_bounds[:, 1].min(), 0.0])
        origin = np.round(all_bounds.mean(axis=0), 0)

    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

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:.4f} 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

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.

Explanation

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.

Bars of float32 representable spacing at 1,000, 100,000, 457,900 and 1,000,000 metres, rising from 0.06 mm to 63 mm.
At national-grid magnitudes the quantisation is centimetres, which is visible.

Edge cases or notes

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