Parsing CityGML with lxml and Shapely

To parse CityGML, select a namespace map from the file’s root element, iterate features with iterparse while clearing processed subtrees, decode each gml:posList using its declared srsDimension, and assemble the rings into validated Shapely polygons. The two mistakes that account for most first attempts are unqualified XPath, which matches nothing, and a hard-coded coordinate stride, which reads interleaved nonsense. This page sits under CityGML and GML Interchange.

How the Format Resists a Naive Read

Three properties of the format shape the parser.

Namespaces are mandatory and version-specific. Every element belongs to a namespace, and the URIs differ between CityGML 1.0, 2.0 and 3.0. An XPath without a namespace map matches nothing, and a map for the wrong version matches nothing either — in both cases returning an empty result rather than an error, which sends people to inspect the file.

Coordinates arrive as text with a declared dimension. A gml:posList is whitespace-separated numbers, and the srsDimension attribute says whether to group them in twos or threes. Assuming three on a two-dimensional list produces coordinates built from ordinates of successive points: numeric, parseable, meaningless.

Axis order follows the declared CRS. GML honours the authority axis order, so a geographic system yields latitude first. Shapely, GeoJSON and every downstream consumer expect the opposite. The swap has to happen once, driven by the declared srsName, rather than being applied by feel.

Production-Ready Script

# lxml>=4.9, shapely>=2.0, pyproj>=3.5, Python 3.9+
from __future__ import annotations

from dataclasses import dataclass
from lxml import etree
from pyproj import CRS
from shapely.geometry import Polygon
from shapely.validation import make_valid

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


@dataclass(frozen=True)
class CityFeature:
    gml_id: str
    polygons: list[Polygon]


def namespaces(path: str) -> dict:
    """Version detection: read only the root element."""
    for _, root in etree.iterparse(path, events=("start",)):
        uri = etree.QName(root).namespace
        ns = NS_BY_CORE.get(uri)
        if ns is None:
            raise ValueError(f"unmapped CityGML namespace {uri!r}")
        return {**ns, "gml": GML}
    raise ValueError("empty document")


def axis_swap_needed(srs_name: str | None) -> bool:
    if not srs_name:
        return False
    first = CRS.from_user_input(srs_name).axis_info[0].abbrev.lower()
    return first in {"lat", "n"}


def decode_pos_list(el, swap: bool) -> list[tuple[float, float]]:
    dim = int(el.get("srsDimension", "3"))
    vals = [float(v) for v in (el.text or "").split()]
    if not vals or len(vals) % dim:
        raise ValueError(f"posList of {len(vals)} values is not a multiple of {dim}")
    pts = [(vals[i], vals[i + 1]) for i in range(0, len(vals), dim)]
    return [(y, x) for x, y in pts] if swap else pts


def polygons_of(element, ns: dict, swap: bool) -> list[Polygon]:
    out = []
    for poly_el in element.findall(f".//{{{GML}}}Polygon"):
        ext = poly_el.find(f".//{{{GML}}}exterior//{{{GML}}}posList")
        if ext is None:
            continue
        holes = [decode_pos_list(h, swap) for h in
                 poly_el.findall(f".//{{{GML}}}interior//{{{GML}}}posList")]
        poly = Polygon(decode_pos_list(ext, swap), holes)
        if not poly.is_valid:
            poly = make_valid(poly)
        if not poly.is_empty:
            out.append(poly)
    return out


def iter_buildings(path: str):
    """Incremental, namespace-aware, memory-bounded."""
    ns = namespaces(path)
    swap = None
    tag = f"{{{ns['bldg']}}}Building"
    for _, el in etree.iterparse(path, events=("end",), tag=tag):
        if swap is None:
            srs = el.get("srsName") or _inherited_srs(el)
            swap = axis_swap_needed(srs)
        gml_id = el.get(f"{{{GML}}}id")
        if not gml_id:
            raise ValueError("city object without a gml:id — output cannot be reconciled")
        yield CityFeature(gml_id, polygons_of(el, ns, swap))
        el.clear()
        while el.getprevious() is not None:
            del el.getparent()[0]              # release processed siblings


def _inherited_srs(el) -> str | None:
    node = el
    while node is not None:
        srs = node.get("srsName")
        if srs:
            return srs
        node = node.getparent()
    return None


if __name__ == "__main__":
    for feature in iter_buildings("city.gml"):
        total = sum(p.area for p in feature.polygons)
        print(f"{feature.gml_id}: {len(feature.polygons)} surface(s), {total:.1f} m2")

Key implementation notes:

• The sibling deletion inside iter_buildings is what makes iterparse bounded. Clearing the element alone is not enough — the parent keeps every processed sibling alive.
• srsName is looked up on the feature and then inherited from ancestors, because it is normally declared once on an enclosing envelope rather than per feature.
• The axis decision is made once and reused. Testing per ring would be correct and would spend the whole parse in pyproj.
• Missing gml:id raises. A feature that cannot be reconciled with its source is not usable output, however good its geometry.

Compatibility Matrix

Fallback Strategies

1. No features found. Almost always the namespace map. Print the root namespace and compare it against NS_BY_CORE; a 3.0 file against a 2.0 map yields exactly this symptom.

2. Coordinates in the wrong hemisphere. Axis order. Confirm what srsName declares rather than swapping until the map looks right, because a swap applied to a projected system breaks a file that was previously correct.

3. srsDimension absent. The default in the code above is three, which matches most CityGML. Where a producer omits the attribute on two-dimensional data, the stride is wrong for the whole file — assert that the resulting coordinate count is plausible for the feature.

4. Deprecated gml:coordinates. Older exports use it, with a comma between ordinates. Add a branch that splits on the declared separators; the rest of the pipeline is unchanged.

5. Enormous single features. A terrain surface can be one feature with millions of vertices, and iterparse bounds memory per feature rather than within one. Where a single feature does not fit, the file needs tiling before parsing rather than a better parser.

FAQ

Related Pages

• CityGML and GML Interchange — parent reference on levels of detail and the GML geometry model
• Converting IFC Buildings to CityGML LoD1 with Python — the writing counterpart to this reading guide
• Mapping GML Geometry to PostGIS with OGR — where the parsed features usually go next