CityGML and GML Interchange

CityGML is the open standard for representing objects in a city — buildings, terrain, transport, vegetation — as thematically classified features with geometry at a declared level of generalisation, and reading or producing it in Python is where a building model becomes a city object. It belongs to the Core Format Fundamentals & Schema Mapping section because everything difficult about it is a schema question rather than a geometry one.

The difficulty is not parsing XML. It is that CityGML asks a question the source data usually cannot answer: at what level of detail is this geometry, and what is this object thematically? A design model knows it contains a wall; it does not know whether the wall is part of a building whose LoD2 envelope should include it. Getting that mapping wrong produces a file that validates against the schema and is useless for analysis, which is the worst of both outcomes.

Prerequisites

• Python 3.9+.
• lxml>=4.9 — the parser used throughout, for its XPath and namespace support. The standard library’s ElementTree can read these files but its XPath subset makes real CityGML traversal painful.
• shapely>=2.0 for the geometry representation on the Python side.
• pyproj>=3.5 where the source and target coordinate reference systems differ.
• Familiarity with XML namespaces. Every element in a CityGML document is namespaced, and this is the single most common reason a first attempt returns nothing.

# lxml>=4.9  shapely>=2.0  pyproj>=3.5
pip install "lxml>=4.9" "shapely>=2.0" "pyproj>=3.5"

Architectural Overview

A CityGML document is a GML feature collection. GML supplies the geometry primitives and the feature machinery; CityGML supplies the thematic modules — Building, Transportation, Vegetation, Relief and others — that classify what those geometries represent.

Three properties of that structure govern how it is parsed.

Everything is namespaced, and the namespaces are version-specific. The CityGML core, the Building module and GML itself each have their own URI, and those URIs change between CityGML 1.0, 2.0 and 3.0. A namespace map is therefore not boilerplate; it is version detection.

Geometry is expressed as GML primitives, not as coordinate arrays. A surface is a gml:Polygon whose exterior is a gml:LinearRing whose points arrive as a whitespace-separated gml:posList, with an srsDimension attribute saying whether the values come in twos or threes. There is no shortcut to a coordinate list; the nesting is the format.

The level of detail is declared per geometry property, not per file. A single building can carry an LoD0 footprint and an LoD2 solid simultaneously, in different properties. A parser that takes the first geometry it finds is picking a generalisation level at random.

The mapping into a Python pipeline therefore has three jobs: select a namespace map from the file’s version, walk from city objects down to the geometry property at the level of detail you want, and decode the coordinate strings into something Shapely can hold.

Step-by-Step Implementation

1. Detect the version and build the namespace map

# lxml>=4.9
from lxml import etree

GML = "http://www.opengis.net/gml"

NS_BY_VERSION = {
    "http://www.opengis.net/citygml/2.0": {
        "core": "http://www.opengis.net/citygml/2.0",
        "bldg": "http://www.opengis.net/citygml/building/2.0",
        "gml": GML,
    },
    "http://www.opengis.net/citygml/1.0": {
        "core": "http://www.opengis.net/citygml/1.0",
        "bldg": "http://www.opengis.net/citygml/building/1.0",
        "gml": GML,
    },
}

def open_citygml(path: str):
    tree = etree.parse(path)
    root = tree.getroot()
    uri = etree.QName(root).namespace
    ns = NS_BY_VERSION.get(uri)
    if ns is None:
        raise ValueError(f"unmapped CityGML namespace {uri!r} — refusing to guess")
    return tree, ns

Refusing an unknown namespace rather than falling back is deliberate. A fallback map produces a parse that finds no buildings, which reads exactly like an empty file.

2. Walk from city objects to the level of detail you want

# lxml>=4.9
LOD_PROPERTIES = ["lod2Solid", "lod2MultiSurface", "lod1Solid", "lod0FootPrint"]

def buildings(tree, ns):
    return tree.findall(".//bldg:Building", namespaces=ns)

def preferred_geometry(building, ns):
    """Return (lod_name, element) for the most detailed geometry present."""
    for prop in LOD_PROPERTIES:
        found = building.find(f"bldg:{prop}", namespaces=ns)
        if found is not None:
            return prop, found
    return None, None

The order of LOD_PROPERTIES is a policy decision and belongs in configuration, not in a parser. A visualisation pipeline wants the most detailed geometry available; a footprint pipeline wants LoD0 even when LoD2 exists, because extracting a footprint from a solid is work it does not need to do.

3. Decode GML coordinate strings

# lxml>=4.9, shapely>=2.0
from shapely.geometry import Polygon

def ring_coords(ring_el, ns):
    """A gml:LinearRing -> a list of (x, y) tuples."""
    pos_list = ring_el.find(".//gml:posList", namespaces=ns)
    if pos_list is None:
        raise ValueError("LinearRing without a posList")
    dim = int(pos_list.get("srsDimension", "3"))
    values = [float(v) for v in pos_list.text.split()]
    if len(values) % dim:
        raise ValueError(f"posList length {len(values)} is not a multiple of {dim}")
    pts = [tuple(values[i:i + dim]) for i in range(0, len(values), dim)]
    return [(p[0], p[1]) for p in pts]          # drop Z for a planar footprint

def polygon_from_surface(surface_el, ns) -> Polygon:
    exterior = surface_el.find(".//gml:exterior//gml:LinearRing", namespaces=ns)
    interiors = surface_el.findall(".//gml:interior//gml:LinearRing", namespaces=ns)
    return Polygon(ring_coords(exterior, ns),
                   [ring_coords(i, ns) for i in interiors])

4. Handle axis order before anything else touches the coordinates

The srsName attribute names the coordinate reference system, and GML follows that system’s declared axis order. For a geographic CRS that is latitude then longitude — the reverse of what Polygon above assumes.

# pyproj>=3.5
from pyproj import CRS

def needs_axis_swap(srs_name: str) -> bool:
    """True when the declared CRS puts latitude/northing first."""
    crs = CRS.from_user_input(srs_name)
    first = crs.axis_info[0].abbrev.lower()
    return first in {"lat", "n"}

Read srsName from the outermost element that declares it and apply the decision once, rather than testing per ring.

5. Preserve identity on the way out

# lxml>=4.9
def city_object_id(el, ns) -> str:
    gml_id = el.get(f"{{{ns['gml']}}}id")
    if not gml_id:
        raise ValueError("city object without a gml:id — output cannot be reconciled")
    return gml_id

An output whose features cannot be matched back to their source is a dead end: a discrepancy found later cannot be attributed to an object. Treat the identifier as required rather than optional, and carry it into whatever the pipeline writes.

Edge Cases and Gotchas

Unqualified XPath matches nothing. The most common first failure. .//Building finds no elements in a namespaced document; .//bldg:Building with a namespace map finds them all. The symptom is an empty result rather than an error, which sends people looking at the file instead of the query.

srsDimension is not always three. A posList may be two-dimensional, and a parser that hard-codes a stride of three reads coordinates that are interleaved nonsense — still numeric, still parseable, geometrically meaningless. Read the attribute.

Solids include vertical faces. Projecting an LoD2 solid to a footprint by dropping Z produces the wall faces as degenerate zero-area polygons alongside the roof and floor. Filter degenerate results before the union, exactly as when extracting IFC wall geometries to Shapely.

Ring orientation is not guaranteed. CityGML surfaces carry an orientation through their boundary composition rather than through winding, so rings arrive in either direction. Normalise winding before treating a ring as a polygon boundary, or hole detection becomes unreliable.

Generic attributes hide in a separate module. Thematic attributes defined by a producer rather than by the standard live in gen:stringAttribute and its siblings, not as elements on the feature. A pipeline reading only the standard attributes silently drops everything the producer added — which, in practice, is usually the part the client cares about.

CityGML 3.0 restructures level of detail. LoD4 no longer exists; interior structure is modelled through spaces instead. Code written against 2.0’s LoD4 property finds nothing in a 3.0 file and reports the building as having no interior.

Validation and Testing

The assertions worth writing are about identity and generalisation, because geometry errors in this format tend to be visible and schema errors do not.

# pytest, lxml>=4.9, shapely>=2.0
def test_every_building_has_an_id_and_geometry(citygml_path):
    tree, ns = open_citygml(citygml_path)
    bs = buildings(tree, ns)
    assert bs, "no buildings parsed — check the namespace map"
    for b in bs:
        assert city_object_id(b, ns)
        lod, geom = preferred_geometry(b, ns)
        assert geom is not None, f"{city_object_id(b, ns)} has no geometry at any LoD"

def test_footprints_are_valid_and_plausible(citygml_path):
    tree, ns = open_citygml(citygml_path)
    for b in buildings(tree, ns):
        _, geom_el = preferred_geometry(b, ns)
        for surface in geom_el.findall(".//gml:Polygon", namespaces=ns):
            poly = polygon_from_surface(surface, ns)
            assert poly.is_valid, city_object_id(b, ns)
            assert 1.0 < poly.area < 1e6, "implausible footprint area"

The area bound is doing real work: it catches an axis swap, an unapplied unit and a coordinate stride error at once, because all three produce areas that are absurd rather than merely wrong.

Performance and Scale

City-scale CityGML files reach hundreds of megabytes, and a full DOM parse of one holds the entire tree in memory at several times the file size. For anything beyond a district, parse incrementally and release each feature as it is consumed:

# lxml>=4.9
def iter_buildings(path: str, bldg_ns: str):
    context = etree.iterparse(path, events=("end",), tag=f"{{{bldg_ns}}}Building")
    for _, element in context:
        yield element
        element.clear()                       # release this subtree
        while element.getprevious() is not None:
            del element.getparent()[0]        # release the already-processed siblings

The sibling deletion is the part that is easy to omit and the part that matters: without it the parent accumulates every element already processed, and memory grows exactly as it would have with a full parse.

Beyond parsing, the expensive stage is geometry assembly rather than XML traversal — polygon construction and union dominate once features are in Shapely. Where only footprints are needed, prefer the LoD0 property when it exists rather than projecting a solid, and where the whole city is being processed, partition by tile and process tiles in separate processes; the work is embarrassingly parallel because city objects are independent.

FAQ

Related Pages

• Core Format Fundamentals & Schema Mapping — the parent section on how each format encodes geometry and meaning
• IFC4x3 Schema Mapping — the building-scale schema a city model generalises from
• Metadata Extraction Strategies — the attribute side of the same mapping problem
• DXF vs IFC for GIS Ingestion — choosing the interchange representation before this mapping is written
• GeoPackage vs PostGIS for CAD Output — where the parsed city objects are usually stored