CityJSON coordinates come out tiny or at the origin

Problem statement

You load a CityJSON model, take the vertex array, and the numbers are wrong. Sometimes they are integers in the hundreds of thousands with no decimal part; sometimes they start at zero; sometimes the whole city is a few units across. Plotted, the model is at the origin of the coordinate system rather than in the city it describes.

The cause is always the same: CityJSON stores vertices as integers with a per-file transform, and the transform has not been applied. On a real Hague tile the raw vertices run from [0, 0, 0] to [787364, 671848, 35018]; after the transform they run from [78248.67, 457604.59, 2.46] to [79036.02, 458276.44, 37.48] โ€” the Dutch national grid, which is where the city actually is.

Quick answer

import json, numpy as np

d = json.loads(open("model.city.json").read())
V = np.asarray(d["vertices"], dtype="float64")
if "transform" in d:
    V = V * np.asarray(d["transform"]["scale"]) + np.asarray(d["transform"]["translate"])

That is the whole fix. The transform on the reference tile is:

{'scale': [0.001, 0.001, 0.001], 'translate': [78248.66, 457604.591, 2.463]}

A scale of 0.001 means the integers are millimetres, and the translate puts the local origin on the national grid.

Triage of four symptoms of a missed or mis-applied CityJSON transform and the fix for each.
Four symptoms, one root cause and two near misses.

Step-by-step solution

1. Confirm the symptom against the transform

raw = np.asarray(d["vertices"])
print("raw:", raw.min(axis=0), "โ†’", raw.max(axis=0))
print("transform:", d.get("transform"))

Integer vertices starting near zero plus a translate that looks like a national grid coordinate is the signature.

2. Apply scale then translate, in that order

vertex * scale + translate. Reversing them applies the translation at the wrong magnitude and produces coordinates that are wrong by a factor of the scale.

3. Do not round-trip through float32

The transform exists partly so coordinates are exact integers. Converting to float32 anywhere in the chain reintroduces the precision problem the format was designed to avoid โ€” at a northing of 457,900 m, float32 has a representable spacing of 3.1 cm.

4. Apply it once

Transforming an already-transformed array multiplies a national grid coordinate by 0.001 and adds the translate again, giving a model a few hundred metres from the origin. That is the "city is tiny" symptom rather than the "city is at zero" one.

5. Check the CRS separately

metadata.referenceSystem is an OGC URI, and its trailing number is the EPSG code. โ€ฆ/EPSG/0/7415 is EPSG:7415 โ€” RD New with NAP heights. Reading the transform correctly and then assuming EPSG:4326 puts the model in the Gulf of Guinea.

6. Watch for the same trap in other formats

The pattern โ€” integers plus a scale and offset โ€” is not unique to CityJSON. LAS and LAZ store scaled integer coordinates with exactly the same structure, and laspy applies them for you when you use las.x rather than las.X.

7. Write the transform back out

If you produce a CityJSON file, either write a transform and integer vertices, or write no transform and floats. Writing floats and a transform means every reader applies it to numbers that are already correct.

Diagram showing integer vertices multiplied by scale and added to translate to give national grid coordinates.
Scale first, then translate; the reference tile's scale of 0.001 means the integers are millimetres.

Code examples

Example 1 โ€” a loader that cannot get it wrong

import json, numpy as np

def load_cityjson(path):
    d = json.loads(open(path).read())
    V = np.asarray(d["vertices"], dtype="float64")
    tr = d.get("transform")
    if tr:
        V = V * np.asarray(tr["scale"], dtype="float64") \
            + np.asarray(tr["translate"], dtype="float64")
    crs = None
    ref = (d.get("metadata") or {}).get("referenceSystem")
    if ref:
        crs = int(str(ref).rstrip("/").rsplit("/", 1)[-1])
    return {"doc": d, "vertices": V, "epsg": crs}

m = load_cityjson("DenHaag_01.city.json")
print(m["epsg"], m["vertices"].min(axis=0).round(2), m["vertices"].max(axis=0).round(2))
7415 [ 78248.67 457604.59      2.46] [ 79036.02 458276.44     37.48]

Returning a dict rather than a bare array is deliberate: it makes it awkward to pass the untransformed vertices around by accident.

Example 2 โ€” a sanity check that catches it every time

import numpy as np

def check_vertices(V, epsg=None, expected_bounds=None):
    problems = []
    if np.issubdtype(V.dtype, np.integer):
        problems.append("vertices are integers โ€” the transform has not been applied")
    if np.allclose(V.min(axis=0)[:2], 0, atol=1.0):
        problems.append("minimum x/y is near zero โ€” untransformed local coordinates?")
    extent = V.max(axis=0) - V.min(axis=0)
    if extent[2] > extent[:2].max():
        problems.append(f"z extent {extent[2]:.1f} exceeds the horizontal extent โ€” scale mismatch?")
    if expected_bounds:
        w, s, e, n = expected_bounds
        if not (w <= V[:, 0].min() and V[:, 0].max() <= e):
            problems.append("x is outside the expected bounds for this CRS")
    return problems

print(check_vertices(m["vertices"], epsg=m["epsg"]))

Example 3 โ€” write CityJSON with a transform, correctly

import numpy as np, json

def to_cityjson_vertices(V, precision=0.001):
    """Return integer vertices plus the transform that restores them."""
    translate = V.min(axis=0)
    scaled = np.round((V - translate) / precision).astype("int64")
    return scaled.tolist(), {"scale": [precision] * 3, "translate": translate.tolist()}

verts, transform = to_cityjson_vertices(V)
back = np.asarray(verts) * np.asarray(transform["scale"]) + np.asarray(transform["translate"])
print("max round-trip error:", np.abs(back - V).max())

A precision of 0.001 gives a maximum round-trip error of half a millimetre, which is far below any survey accuracy and shrinks the file substantially compared with writing doubles as text.

Explanation

Why CityJSON stores integers

A coordinate written as 78248.663 is nine characters of JSON; the same value as an integer millimetre count is six. Across millions of vertices that is a large fraction of the file. Storing integers also makes the coordinates exact โ€” there is no decimal round-trip to lose the last digits, which is a real problem for text formats generally.

Why the symptoms differ

Not applying the transform at all leaves integers starting near zero โ€” the model is at the coordinate origin and is the right size in the wrong units. Applying it twice multiplies real coordinates by the scale, shrinking the model to about a thousandth of its size and moving it near the translate. Applying translate before scale scales the translated values, which puts the model near the origin and makes it tiny. The three look different on a plot, which is how you tell them apart.

Why float32 undoes the benefit

The format goes to some trouble to store exact coordinates. Converting the transformed array to float32 for a mesh library reintroduces quantisation at the level the format was avoiding โ€” 3.1 cm at the northings in the reference tile. Keep float64 until the last possible moment, and recentre on a local origin before any float32 conversion.

Why LAS has the same shape of bug

LAS stores X, Y and Z as 32-bit integers with per-file scale factors and offsets, for exactly the same reasons. laspy exposes both las.X (raw integers) and las.x (scaled doubles), and the capitalisation is the only thing distinguishing them โ€” which produces the identical symptom with the identical cause.

Two panels showing the scaled-integer coordinate pattern in CityJSON and in LAS, with the accessor that applies the transform in each.
In laspy the only difference between right and wrong is the capital letter.

Edge cases or notes

  • The transform is optional. A file without one has float vertices already.
  • Scale can differ per axis. Do not assume all three are equal.
  • geographicalExtent is already in real coordinates. Use it as a check.
  • cjio applies the transform when it reads and rewrites one when it writes.
  • CityJSONSeq repeats the transform in its first line; each feature refers to it.
  • Rounding before scaling loses precision. Scale, then round.
  • np.int64 overflow is possible for a very fine scale over a large extent.
  • Record the CRS including its vertical part. EPSG:7415 is not EPSG:28992.

FAQ

Why are my CityJSON coordinates integers?

Because the format stores them that way. Multiply by transform.scale and add transform.translate to get real coordinates.

Why is my model at the origin?

The transform has not been applied. Untransformed vertices start at or near [0, 0, 0].

Why is my model tiny?

The transform has been applied twice, or translate was applied before scale. Both shrink real coordinates by the scale factor.

What does a scale of 0.001 mean?

The integers are thousandths of the CRS unit โ€” millimetres for a metric CRS.

Does the transform affect the CRS?

No. The CRS is in metadata.referenceSystem as an OGC URI whose trailing number is the EPSG code.

Is this the same problem as laspy's X versus x?

Yes, exactly. LAS stores scaled integers with a per-file offset, and las.X gives the raw integers while las.x gives the scaled coordinates.