🧭 Concept 3D & buildings Files & formats 3D & Buildings · step 3 of 20 Updated

CityJSON and CityGML explained

Problem statement

CityGML is the OGC standard for semantic 3D city models: comprehensive, XML, and heavy enough that a city-scale file is measured in tens of gigabytes and parsed by almost nothing. CityJSON is the same data model encoded as JSON, about an order of magnitude smaller, and readable with json.loads.

If you are working in Python, CityJSON is the format you want, and the interesting part is not the encoding — it is the data model both share. Buildings split into parts, geometries as solids with indexed vertices, surfaces labelled by semantic type, and a transform that makes the coordinates look wrong until you apply it.

This guide covers the model, the file structure, and the three things that surprise people on first contact.

Quick answer 

import json, numpy as np

d = json.loads(open("model.city.json").read())
print(sorted(d.keys()))
print(d["version"], d["metadata"].get("referenceSystem"))

# Vertices are integers; apply the transform to get real coordinates
V = np.array(d["vertices"]) * np.array(d["transform"]["scale"]) \
    + np.array(d["transform"]["translate"])
print(V.min(axis=0).round(2), "→", V.max(axis=0).round(2))
['CityObjects', 'appearance', 'metadata', 'transform', 'type', 'version', 'vertices']
2.0 https://www.opengis.net/def/crs/EPSG/0/7415
[ 78248.67 457604.59      2.46] → [ 79036.02 458276.44     37.48]

Without the transform, the same vertices read [0, 0, 0] to [787364, 671848, 35018] — plausible-looking numbers in the wrong place, which is why CityJSON coordinates come out tiny or at the origin exists.

Stack of the top-level members of a CityJSON file: type, version, metadata, transform, vertices, CityObjects and appearance.
Seven top-level members; two of them decide whether your coordinates are right.

Step-by-step solution

1. Read the version and the CRS first

CityJSON 1.0, 1.1 and 2.0 differ in small ways that break code — lod is a number in 1.x and a string in 2.0, and the metadata members were reorganised. The CRS lives in metadata.referenceSystem as an OGC URI: …/EPSG/0/7415 means EPSG:7415, which is RD New plus NAP heights.

2. Apply the transform

Vertices are stored as integers with a scale and translate, which is how the format gets its size advantage. Every coordinate you use must be vertex * scale + translate.

3. Understand the object hierarchy

CityObjects is a flat dictionary keyed by id. Buildings reference their parts through children, and parts reference their parent. A real tile held 844 Building objects, 1,653 BuildingPart objects and one TINRelief — the terrain.

4. Read the geometry as indices, not coordinates

A Solid's boundaries is nested arrays of vertex indices: shell → surface → ring → index. The outer ring is the first array in each surface; any others are holes.

# boundaries[shell][surface][ring][index]
outer_ring = geometry["boundaries"][0][0][0]
coords = V[outer_ring]

5. Use the semantics

geometry["semantics"] has a surfaces list of types and a values array parallel to the boundaries, mapping each surface to its type. That is how you get roof area without inferring it from normals.

6. Know what CityGML adds

CityGML has everything CityJSON has plus ADEs (application domain extensions), richer appearance modelling, and the full XML apparatus. If you must consume CityGML, convert it: citygml-tools and cjio both do CityGML ↔ CityJSON, and working in CityJSON afterwards is faster in every sense.

7. Use cjio for the operations you would otherwise write

Subsetting, reprojecting, upgrading versions, extracting an LoD and validating are all one command in cjio, and all are fiddly to write correctly against the nested boundary arrays.

Diagram of the nested boundaries array from shell to surface to ring to vertex index.
Four levels of nesting, and the innermost values are indices into one shared vertex list.

Code examples

Example 1 — walk the model and summarise it 

import json, collections, numpy as np

def summarise(path):
    d = json.loads(open(path).read())
    V = np.array(d["vertices"], dtype="float64")
    if "transform" in d:
        V = V * np.array(d["transform"]["scale"]) + np.array(d["transform"]["translate"])

    types = collections.Counter(o["type"] for o in d["CityObjects"].values())
    lods, geoms, sem = collections.Counter(), collections.Counter(), collections.Counter()
    for o in d["CityObjects"].values():
        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(axis=0) - V.min(axis=0)).round(1).tolist(),
            "lods": lods.most_common(), "geometry_types": geoms.most_common(),
            "semantic_surfaces": sem.most_common()}

print(json.dumps(summarise("DenHaag_01.city.json"), indent=2, default=str))
"objects": 2498, "vertices": 22997, "extent_m": [787.4, 671.8, 35.0],
"lods": [["2", 1991]], "geometry_types": [["Solid", 1990], ["CompositeSurface", 1]],
"semantic_surfaces": [["WallSurface", 11215], ["RoofSurface", 3004], ["GroundSurface", 1990]]

Example 2 — extract roof surfaces using the semantics 

import numpy as np
from shapely.geometry import Polygon

def surfaces_of_type(obj, V, wanted="RoofSurface"):
    out = []
    for g in obj.get("geometry", []):
        sem = g.get("semantics")
        if not sem:
            continue
        surfaces, values = sem["surfaces"], sem["values"]
        shells = g["boundaries"] if g["type"] == "Solid" else [g["boundaries"]]
        vals = values if g["type"] != "Solid" else values
        for shell, shell_vals in zip(shells, vals if g["type"] == "Solid" else [values]):
            for face, si in zip(shell, shell_vals):
                if si is None or surfaces[si].get("type") != wanted:
                    continue
                out.append(np.array([V[i] for i in face[0]]))
    return out

roofs = [r for o in d["CityObjects"].values() for r in surfaces_of_type(o, V)]
print(f"{len(roofs):,} roof surfaces")

Example 3 — the 3D area of a surface, which is not its footprint area 

import numpy as np

def polygon_area_3d(points):
    """Area of a planar polygon in 3D, from the Newell normal."""
    n = np.zeros(3)
    for i in range(len(points)):
        a, b = points[i], points[(i + 1) % len(points)]
        n += np.cross(a, b)
    return np.linalg.norm(n) / 2

roof_area = sum(polygon_area_3d(r) for r in roofs)
print(f"total roof surface area: {roof_area:,.0f} m²")

A pitched roof's 3D area exceeds its footprint by 1/cos(pitch) — 15% at 30°, 41% at 45°. Using the footprint for a solar calculation underestimates by exactly that factor.

Explanation

Why the transform exists

Storing coordinates as doubles in JSON is expensive: 78248.663 is nine characters, and a city model has millions of them. CityJSON stores integers and a per-file scale and translate, which cuts the file size substantially and makes the numbers exact rather than subject to decimal round-tripping. The cost is that every reader must apply the transform, and a reader that does not produces a model at the origin.

Why buildings are split into parts

The data model says a Building may have BuildingPart children, each with its own geometry, when the building is not a single mass — a house with a lower extension, a block with a tower. Any code that iterates CityObjects and looks only at type == "Building" will miss most of the geometry; in the Hague tile that is 1,653 of 2,497 objects.

Why semantics beat normals

Classifying a surface as roof or wall from its normal fails on nearly vertical roofs, on leaning walls and on surfaces whose winding is inconsistent. The model already knows, because whoever produced it knew. The semantics block is a parallel array to boundaries, which is awkward to index and worth the trouble.

Why CityJSON rather than CityGML in Python

Size and parse cost. The same model is roughly six to ten times smaller as CityJSON, and json.loads is two orders of magnitude faster than an XML parse with namespace handling. Nothing in the data model is lost; CityJSON is a complete encoding of the CityGML model, and round-tripping through citygml-tools is routine.

Bars of semantic surface counts in a real LoD2 tile: 11,215 wall surfaces, 3,004 roof surfaces and 1,990 ground surfaces.
One ground surface per solid, 1.5 roof planes each, and the rest walls.

Edge cases or notes

  • lod is a string in 2.0. Compare as strings, or normalise.
  • CompositeSurface is not a solid. Terrain and some parts use it; it has no volume.
  • Holes are the second and later rings. Ignoring them overstates roof areas.
  • appearance holds textures and materials. Usually ignorable for analysis.
  • geographicalExtent is six numbers, in the model's own CRS.
  • CityJSONSeq exists for streaming: one object per line.
  • cjio validates. Run it before trusting a model you did not produce.
  • The CRS URI includes the vertical datum. EPSG:7415 is not EPSG:28992.

FAQ

What is CityJSON?

A JSON encoding of the CityGML data model: semantic 3D city objects with indexed vertices, solids, and labelled surfaces. It is far smaller than CityGML XML and readable with the standard library.

Why are my CityJSON coordinates tiny?

Because vertices are stored as integers and need the file's transform applied: vertex * scale + translate.

What is the difference between a Building and a BuildingPart?

A Building may be split into parts where its roof or mass changes. The geometry usually lives on the parts, so iterating only over Building objects misses most of the model.

How do I get roof surfaces out of a model?

Use the semantics block: its surfaces list gives the types and values maps each boundary surface to one. Do not classify by surface normal.

Should I use CityGML or CityJSON in Python?

CityJSON. It encodes the same model, is much smaller, and parses with json. Convert CityGML with citygml-tools if that is what you are given.

Is a roof's area the same as its footprint?

No. A pitched roof's 3D area is its footprint divided by the cosine of the pitch — 15% larger at 30° and 41% larger at 45°.