How to read a CityJSON city model in Python
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.
Step-by-step solution
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.
Code examples
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 41.0
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; for a real model with concave roof faces, triangulate with mapbox_earcut or shapely in the surface's own plane. The 99.7% watertight figure above is with the simple fan, which is why five solids fail.
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)}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.
Explanation
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 wrong for a concave one, where some triangles fall outside the polygon. Real roof surfaces are often concave, and the result is a mesh whose volume is close but whose faces overlap — which is why five of the 1,990 Hague solids fail the watertight test under a fan and would pass under an ear-clipping triangulation in the face's own plane.
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.
Edge cases or notes
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.