GML and GeoJSON Payload Handling

Every OGC feature service speaks one of two payload dialects, and a production client eventually has to speak both. WFS returns GML by default; OGC API - Features returns GeoJSON. They disagree about axis order, about how a coordinate list is written, about whether a feature must have an identifier, and about how much of a reference system a payload is allowed to carry. This guide covers the structural differences that actually break code, how to stream both without exhausting memory, and how to convert between them without silently losing the parts GeoJSON cannot express.

Prerequisites & Architecture Context

You need lxml for GML — the standard library’s ElementTree will parse it, but iterparse with element clearing is what keeps memory flat on a large response, and lxml implements it more efficiently. For GeoJSON you need ijson if payloads may exceed a few tens of megabytes, and shapely on both paths as the common geometry representation. pyproj handles reference system conversion where a service returns coordinates in something other than CRS84.

The architectural decision to make early is where normalisation happens. Two viable designs exist: convert at the edge, so everything downstream sees one internal representation; or carry the native encoding through and convert at each consumer. The first is almost always right — a single normalisation point is where axis order can be fixed once, and the SRS and Coordinate Reference System Handling guide explains why leaving that decision to individual consumers guarantees it will be made inconsistently.

Specification Deep-Dive: Where the Encodings Diverge

The differences below are not stylistic. Each one has a corresponding class of bug that appears only against real data.

Coordinates. GML 3.2 writes a geometry’s coordinates as a whitespace-separated gml:posList, with the number of ordinates per position given by a srsDimension attribute that defaults to 2. A three-dimensional geometry therefore looks identical to a two-dimensional one apart from that attribute, and a parser that assumes pairs will silently read every third value as a latitude. GeoJSON nests arrays and states dimensionality by array length, which removes the ambiguity but makes streaming harder.

Axis order. This is the difference that costs the most engineering time. GML coordinate order follows the axis definition of the srsName on the geometry, so urn:ogc:def:crs:EPSG::4326 means latitude first. RFC 7946 GeoJSON has no axis question: coordinates are always longitude then latitude, and the only permitted reference system is CRS84. A converter that copies ordinates across without consulting srsName produces geometries in the Indian Ocean, which is why the reprojection guide treats always_xy as non-negotiable.

Feature identity. gml:id is mandatory on every GML feature and unique within the document — it is what a WFS-T Update or Delete targets. GeoJSON’s id member is optional, and many producers omit it, putting the identifier in properties instead or not at all. A pipeline that reads GeoJSON and later needs to write back through WFS-T must preserve that identifier deliberately.

Schema. A GML feature type validates against an XSD that the service publishes through DescribeFeatureType, which means the attribute names, types and cardinalities are discoverable before a single feature is fetched. GeoJSON validates against one generic schema that says nothing about which properties a given collection carries, so the equivalent knowledge has to come from the OGC API - Features collection metadata or from inspecting the data.

Python Implementation: A Streaming Ingest Path

The pattern below reads either encoding without materialising the whole response, and converges on a single internal feature type.

from __future__ import annotations

from dataclasses import dataclass, field
from typing import Any, Iterator

import requests
from lxml import etree
from shapely.geometry import shape as shapely_shape

GML_NS = "http://www.opengis.net/gml/3.2"
WFS_NS = "http://www.opengis.net/wfs/2.0"

# Reference systems whose authority axis order is latitude first. GML honours
# the authority definition; GeoJSON never does.
LAT_FIRST = {
    "urn:ogc:def:crs:EPSG::4326",
    "http://www.opengis.net/def/crs/EPSG/0/4326",
}

@dataclass
class Feature:
    """One normalised feature: geometry always longitude-first."""
    fid: str | None
    geometry: Any
    properties: dict[str, Any] = field(default_factory=dict)
    srs: str = "CRS84"

def stream_gml(url: str, params: dict[str, str], chunk: int = 65536) -> Iterator[Feature]:
    """Yield features from a WFS GML response without buffering the document.

    `iterparse` fires on the closing tag of each member, so the tree for a
    single feature is complete when we see it — and clearing it there keeps
    peak memory proportional to one feature, not to the response.
    """
    with requests.get(url, params=params, stream=True, timeout=120) as resp:
        resp.raise_for_status()
        resp.raw.decode_content = True
        context = etree.iterparse(
            resp.raw, events=("end",), tag=f"{{{WFS_NS}}}member", huge_tree=False)
        for _, member in context:
            for child in member:
                yield _gml_feature(child)
            member.clear()
            while member.getprevious() is not None:
                del member.getparent()[0]

def _gml_feature(node: etree._Element) -> Feature:
    fid = node.get(f"{{{GML_NS}}}id")
    geom_node = next((el for el in node.iter()
                      if el.tag.startswith(f"{{{GML_NS}}}")
                      and el.tag.rsplit("}", 1)[1] in _GEOMETRY_TAGS), None)
    srs = (geom_node.get("srsName") if geom_node is not None else None) or "CRS84"
    geometry = _parse_gml_geometry(geom_node, srs) if geom_node is not None else None
    props = {
        el.tag.rsplit("}", 1)[1]: (el.text or "").strip()
        for el in node
        if not el.tag.startswith(f"{{{GML_NS}}}") and len(el) == 0
    }
    return Feature(fid=fid, geometry=geometry, properties=props, srs=srs)

_GEOMETRY_TAGS = {"Point", "LineString", "Polygon", "MultiSurface",
                  "MultiCurve", "MultiPoint", "Surface", "Curve"}

def _positions(text: str, dimension: int, lat_first: bool) -> list[tuple[float, ...]]:
    """Split a posList into positions, fixing axis order as we go."""
    flat = [float(v) for v in text.split()]
    if len(flat) % dimension:
        raise ValueError(f"posList holds {len(flat)} ordinates, not a multiple "
                         f"of srsDimension={dimension}")
    out = []
    for i in range(0, len(flat), dimension):
        pos = tuple(flat[i:i + dimension])
        out.append((pos[1], pos[0], *pos[2:]) if lat_first else pos)
    return out

def _parse_gml_geometry(node: etree._Element, srs: str):
    """Handle the common cases; delegate anything exotic to a full GML reader."""
    from shapely.geometry import LineString, Point, Polygon

    lat_first = srs in LAT_FIRST
    local = node.tag.rsplit("}", 1)[1]
    dim = int(node.get("srsDimension", "2"))

    if local == "Point":
        pos = node.find(f"{{{GML_NS}}}pos")
        return Point(_positions(pos.text, dim, lat_first)[0])
    if local == "LineString":
        pl = node.find(f"{{{GML_NS}}}posList")
        return LineString(_positions(pl.text, dim, lat_first))
    if local == "Polygon":
        ext = node.find(f".//{{{GML_NS}}}exterior//{{{GML_NS}}}posList")
        holes = node.findall(f".//{{{GML_NS}}}interior//{{{GML_NS}}}posList")
        return Polygon(_positions(ext.text, dim, lat_first),
                       [_positions(h.text, dim, lat_first) for h in holes])
    raise NotImplementedError(f"geometry type {local} needs a full GML reader")

def stream_geojson(url: str, params: dict[str, str]) -> Iterator[Feature]:
    """Yield features from a GeoJSON FeatureCollection, one at a time."""
    import ijson

    with requests.get(url, params=params, stream=True, timeout=120) as resp:
        resp.raise_for_status()
        for item in ijson.items(resp.raw, "features.item"):
            geom = item.get("geometry")
            yield Feature(
                fid=str(item["id"]) if "id" in item else None,
                geometry=shapely_shape(geom) if geom else None,
                properties=dict(item.get("properties") or {}),
                srs="CRS84",
            )

The two readers return the same Feature type, which is the whole point: every consumer downstream — a validator, a database writer, a tile renderer — is written once against a representation whose axis order is already settled.

Error Handling & Edge Cases

srsDimension defaults to 2 and is often wrong. A service publishing elevation data may emit three-dimensional positions without setting the attribute on every geometry, because it is permitted to set it once on an enclosing element. The _positions helper raises when the ordinate count is not a multiple of the assumed dimension, which converts a silent misparse into a loud failure — the alternative is coordinates shifted by one position for the rest of the geometry.

Empty geometries are legal. GML permits a feature with no geometry property at all, and GeoJSON permits "geometry": null. Both are common in metadata-only records. The Feature type above models this as None rather than an empty geometry, because an empty Polygon and a missing Polygon mean different things to a spatial index.

iterparse leaks without explicit clearing. Clearing the member element is not enough — lxml keeps the preceding siblings alive on the parent, so a long response still accumulates the whole document. Deleting the earlier siblings, as the loop above does, is what makes memory genuinely flat.

Namespace versions differ. GML 3.1.1 uses http://www.opengis.net/gml while GML 3.2 uses http://www.opengis.net/gml/3.2, and WFS 1.1.0 responses wrap features in gml:featureMember where WFS 2.0 uses wfs:member. A reader hard-coded to one pair silently yields zero features against the other — the same class of version trap covered in WFS 2.0 vs 1.1.0 Breaking Changes.

Testing & Compliance Verification

Test the axis-order fix explicitly, because it is the one that produces plausible-looking wrong answers rather than exceptions.

import pytest
from shapely.geometry import Point

GML_LAT_FIRST = b'''<wfs:FeatureCollection
    xmlns:wfs="http://www.opengis.net/wfs/2.0"
    xmlns:gml="http://www.opengis.net/gml/3.2">
  <wfs:member>
    <app:station gml:id="s.1" xmlns:app="http://example.org/app">
      <app:location>
        <gml:Point srsName="urn:ogc:def:crs:EPSG::4326">
          <gml:pos>47.3769 8.5417</gml:pos>
        </gml:Point>
      </app:location>
      <app:name>Zurich</app:name>
    </app:station>
  </wfs:member>
</wfs:FeatureCollection>'''

def test_lat_first_srs_is_swapped(tmp_path):
    path = tmp_path / "fc.gml"
    path.write_bytes(GML_LAT_FIRST)
    features = list(_stream_local(path))
    assert features[0].geometry.equals_exact(Point(8.5417, 47.3769), 1e-9)
    assert features[0].fid == "s.1"
    assert features[0].properties["name"] == "Zurich"

Assert on longitude explicitly rather than on “the first ordinate”: a test written as assert coords[0] == 47.3769 passes against the bug it was written to catch. For formal conformance, the OGC CITE WFS 2.0 suite exercises GML output against the published schemas, and the OGC API - Features suite checks GeoJSON structure — neither validates axis order against your expectations, only against the specification’s.

Performance & Scaling Notes

Stream both encodings. A GetFeature over a national dataset returns hundreds of megabytes of GML. etree.fromstring on that response is a multi-gigabyte resident set; iterparse with clearing is a few megabytes regardless of length. The same asymmetry applies to json.loads versus ijson on the GeoJSON side.

Ask for less rather than parsing faster. Both protocols support restricting the returned properties — propertyName in WFS, and increasingly properties in OGC API - Features. Parsing cost is dominated by the number of elements, so dropping twenty unused attributes is a larger win than any parser optimisation.

Prefer GeoJSON when you have the choice. For a client that only needs geometry and flat attributes, GeoJSON parses several times faster than the equivalent GML and carries no axis ambiguity. Reach for GML when you need gml:id round-tripping for transactions, a non-CRS84 reference system, or the schema information DescribeFeatureType provides.

Push paging, not filtering, to the client. A single enormous response is worse than many bounded ones even when total bytes are identical, because failure in the middle of a stream costs the whole run. The paging discipline in Paginating OGC API - Features Collections in Python applies equally to WFS with count and startIndex.

Frequently Asked Questions

---

Back to OGC Standards Architecture & Service Fundamentals

Related

• WFS Transactional Operations Deep Dive — the service that returns GML by default
• OGC API - Features and the REST Transition — the service that returns GeoJSON by default
• Parsing GML 3.2 Geometries Into Shapely Objects — the geometry reader in full
• Streaming Large GeoJSON Feature Collections in Python — constant-memory ingest on the JSON side
• SRS and Coordinate Reference System Handling — why axis order must be settled at one boundary