How to read a CityJSON city model in Python
On this page
Problem statement
CityJSON is JSON, so reading it is json.loads. Getting anything useful out of it is not, because the geometry is four levels of nested indices into a shared vertex list, the coordinates need a transform applied, and the surface labels live in a parallel array that has to be zipped against the boundaries.
None of it is difficult, and all of it is easy to get subtly wrong — a reader that ignores the transform produces a model at the origin, one that ignores interior rings overstates roof areas, and one that looks only at Building objects misses most of the geometry.
This guide builds a reader that handles all three, tested against a real LoD2 tile of The Hague: 2,498 city objects, 1,990 solids and 22,997 vertices.
Quick answer
import json, numpy as np
def load_cityjson(path):
d = json.loads(open(path).read())
V = np.asarray(d["vertices"], dtype="float64")
if "transform" in d:
V = V * np.asarray(d["transform"]["scale"]) + np.asarray(d["transform"]["translate"])
return d, V
d, V = load_cityjson("DenHaag_01.city.json")
print(d["version"], d["metadata"].get("referenceSystem"), len(d["CityObjects"]), len(V))2.0 https://www.opengis.net/def/crs/EPSG/0/7415 2498 22997
Applying the transform is the first thing and the thing most readers forget: without it the vertices run from [0, 0, 0] to [787364, 671848, 35018] instead of sitting on the Dutch national grid.
Transform, CRS, every object, nested boundaries
1. Load and transform
Vertices are integers. Multiply by transform.scale and add transform.translate.
2. Read the CRS from the metadata
metadata.referenceSystem is an OGC URI. https://www.opengis.net/def/crs/EPSG/0/7415 means EPSG:7415 — RD New with NAP heights. Parse the trailing number rather than assuming.
crs = int(d["metadata"]["referenceSystem"].rsplit("/", 1)[-1])3. Walk every object, not only Buildings
CityObjects is flat. The Hague tile holds 844 Building objects, 1,653 BuildingPart objects and one TINRelief. The geometry is on the parts; a loop filtering on type == "Building" processes almost nothing.
4. Understand the boundaries nesting
For a Solid: boundaries[shell][surface][ring][vertex_index]. Shell 0 is the exterior; further shells are interior voids. Ring 0 of a surface is its outer ring; further rings are holes.
5. Zip the semantics against the boundaries
semantics.values mirrors the boundaries structure and holds an index into semantics.surfaces. For a Solid, values[shell][surface].
sem = geometry["semantics"]
for face, si in zip(geometry["boundaries"][0], sem["values"][0]):
kind = sem["surfaces"][si]["type"] if si is not None else None6. Decide what you want out
- Footprints + heights → take the
GroundSurfaceof each solid and the height attributes. - Roof surfaces → take the
RoofSurfacefaces and compute 3D areas. - Meshes → fan-triangulate every surface.
7. Use cjio for anything structural
Subsetting by bounding box, reprojecting, upgrading a version, extracting a single LoD and validating are all one command and all tedious to write against the nested arrays.
Footprints, heights and watertight meshes
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, CityJSON 2.0, 3.2 MB), kept under its own name. It covers about 790 × 670 m of the city in EPSG:7415 (the Dutch national grid with NAP heights).
import os
import requests
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)
print(os.path.getsize("DenHaag_01.city.json"))3150451
Example 1 — footprints and heights as a GeoDataFrame
import json, numpy as np, geopandas as gpd
from shapely.geometry import Polygon
def footprints(path):
d, V = load_cityjson(path)
crs = int(d["metadata"]["referenceSystem"].rsplit("/", 1)[-1])
rows = []
for oid, obj in d["CityObjects"].items():
attrs = obj.get("attributes") or {}
for g in obj.get("geometry", []):
if g["type"] != "Solid" or not g.get("semantics"):
continue
surfaces, values = g["semantics"]["surfaces"], g["semantics"]["values"][0]
for face, si in zip(g["boundaries"][0], values):
if si is None or surfaces[si].get("type") != "GroundSurface":
continue
outer = face[0]
if len(outer) < 3:
continue
holes = [[V[i][:2] for i in ring] for ring in face[1:]]
rows.append({
"id": oid, "parent": (obj.get("parents") or [None])[0],
"lod": str(g.get("lod")),
"roof_type": attrs.get("roofType"),
"ground_z": float(np.mean([V[i][2] for i in outer])),
"eaves": attrs.get("AbsoluteEavesHeight"),
"ridge": attrs.get("AbsoluteRidgeHeight"),
"geometry": Polygon([V[i][:2] for i in outer], holes),
})
gdf = gpd.GeoDataFrame(rows, crs=crs)
gdf["height_eaves_m"] = gdf["eaves"] - gdf["ground_z"]
gdf["height_ridge_m"] = gdf["ridge"] - gdf["ground_z"]
return gdf
b = footprints("DenHaag_01.city.json")
print(len(b), b.height_ridge_m.median().round(2), b.area.median().round(1))1990 6.68 24.4
Example 2 — solids as meshes, with the watertight check
import numpy as np, trimesh
def solids(path):
d, V = load_cityjson(path)
out = {}
for oid, obj in d["CityObjects"].items():
for g in obj.get("geometry", []):
if g["type"] != "Solid":
continue
tris = []
for shell in g["boundaries"]:
for surface in shell:
ring = surface[0]
for i in range(1, len(ring) - 1):
tris.append([ring[0], ring[i], ring[i + 1]])
if not tris:
continue
m = trimesh.Trimesh(vertices=V, faces=np.array(tris), process=True)
m.remove_unreferenced_vertices()
out[oid] = m
return out
meshes = solids("DenHaag_01.city.json")
wt = sum(m.is_watertight for m in meshes.values())
vols = np.array([abs(m.volume) for m in meshes.values()])
print(f"{len(meshes):,} solids, {wt:,} watertight ({wt/len(meshes):.1%})")
print(f"volume: median {np.median(vols):.1f} m³, total {vols.sum():,.0f} m³")1,990 solids, 1,985 watertight (99.7%)
volume: median 118.4 m³, total 683,115 m³
Fan triangulation is only correct for convex rings, and it turns a ring that repeats a vertex index into a zero-area triangle. The five failures here are the second kind: each has two rings that repeat a vertex (the watertight guide traces them). Dropping degenerate faces, or ear-clipping with mapbox_earcut in each surface's own plane, closes all five.
Example 3 — summarise an unfamiliar model before writing any analysis
import collections
def describe(path):
d, V = load_cityjson(path)
types = collections.Counter(o["type"] for o in d["CityObjects"].values())
lods, geoms, sem, attrs = (collections.Counter() for _ in range(4))
for o in d["CityObjects"].values():
for k in (o.get("attributes") or {}):
attrs[k] += 1
for g in o.get("geometry", []):
lods[str(g.get("lod"))] += 1
geoms[g["type"]] += 1
for s in (g.get("semantics") or {}).get("surfaces", []):
sem[s.get("type")] += 1
return {"version": d.get("version"),
"crs": d.get("metadata", {}).get("referenceSystem"),
"objects": len(d["CityObjects"]), "types": types.most_common(),
"vertices": len(V),
"extent_m": (V.max(0) - V.min(0)).round(1).tolist(),
"lods": lods.most_common(), "geometry": geoms.most_common(),
"semantics": sem.most_common(), "attributes": attrs.most_common(8)}
summary = describe("DenHaag_01.city.json")
print(summary["semantics"])
print(summary["attributes"])[('WallSurface', 11215), ('RoofSurface', 3004), ('GroundSurface', 1990)]
[('roofType', 1990), ('RelativeEavesHeight', 1990), ('RelativeRidgeHeight', 1990), ('AbsoluteEavesHeight', 1990), ('AbsoluteRidgeHeight', 1990)]
On the Hague tile this reports [('WallSurface', 11215), ('RoofSurface', 3004), ('GroundSurface', 1990)] and attributes roofType, RelativeEavesHeight, RelativeRidgeHeight, AbsoluteEavesHeight, AbsoluteRidgeHeight on all 1,990 solids. Five minutes of describing saves an afternoon of wrong assumptions.
Shared vertices and fan triangulation
Why the vertex list is shared
Adjacent buildings and adjacent faces of the same building share corners. A shared, indexed vertex list means each coordinate appears once, which is why a model with 1,990 solids and 16,000 surfaces has only 22,997 vertices. It also means a mesh built from a subset must either carry the whole vertex array or be re-indexed — remove_unreferenced_vertices does the latter.
Why fan triangulation is not quite enough
Splitting a ring into triangles from its first vertex is correct for a convex ring and can be wrong for a concave one, where some triangles fall outside the polygon. On the Hague tile that did no harm — 478 watertight solids have a concave ring — but ten rings repeat a vertex index, and fanning those produces zero-area triangles. Those are the five solids that fail the watertight test; dropping degenerate faces, or ear-clipping each face in its own plane, makes them pass.
Why the height attributes are worth using
AbsoluteEavesHeight and AbsoluteRidgeHeight are producer-supplied measurements in the model's vertical datum, and they are more reliable than anything you would derive from the geometry — they were computed from the point cloud the model was built from. Subtract the ground surface's mean z to get a height above ground.
Why to check parents and children
A BuildingPart records its Building in parents. Reporting per building rather than per part means grouping by that value; reporting per part is what you want for geometry. Mixing the two produces totals that are correct and counts that are not.
String levels, voids and semantics
lodis a string in CityJSON 2.0 and a number in 1.x.CompositeSurfacehas no shells. Its boundaries are one level shallower than a Solid's.MultiSurfaceis common for terrain.- Interior shells are voids. Rare in buildings, real in some models.
semantics.valuescan contain nulls for unlabelled surfaces.- Attributes vary by producer. Never assume a key exists.
cjiovalidates and subsets. Use it before writing a custom filter.- CityJSONSeq streams one object per line for very large models.
Internal links
- CityJSON and CityGML explained — the data model behind the reader
- CityJSON coordinates come out tiny or at the origin — the transform, missed
- Level of detail explained for 3D city models — reading the
lodfield - How to calculate building volumes and floor area in Python — the next step
- A building mesh is not watertight and its volume is wrong — the triangulation problem
- 3D spatial data models explained: 2.5D, meshes and solids — what a Solid is
- How to explore an unfamiliar spatial dataset in Python — the same discipline for vector data
- How to estimate rooftop solar potential in Python — using the roof surfaces
FAQ
How do I read CityJSON in Python?
json.loads the file, then multiply the vertex array by transform.scale and add transform.translate. Everything else is indexing into that array.
Why is the geometry an array of numbers?
They are indices into the shared vertex list. boundaries[shell][surface][ring][index] for a Solid.
How do I get building footprints out?
Take the GroundSurface faces identified by the semantics block, and build a polygon from the outer ring with the later rings as holes.
Why does my code find almost no geometry?
Because it is filtering on type == "Building". In a real tile, 1,653 of 2,497 objects were BuildingPart, and that is where the geometry lives.
Is fan triangulation good enough?
Almost. On a real tile it produced 1,985 watertight meshes out of 1,990; the failures are concave roof faces. Triangulate in the face's plane for a correct result.
Should I write my own subsetter?
Not for bounding-box subsets, reprojection, LoD extraction or validation — cjio does all of those and gets the nested arrays right.