Locate the geometry element by local name, read srsName from it or the nearest ancestor that declares one, read srsDimension and assert the ordinate count is a multiple of it, chunk the posList into positions, and swap the first two ordinates when the reference system’s authority definition is latitude-first. Everything else about GML geometry parsing is bookkeeping.
A gml:posList is a whitespace-separated list of numbers. Nothing in it says how many ordinates make up a position, what order they are in, or what they mean. Those three facts live in attributes — srsDimension and srsName — that may sit on the geometry element, on an enclosing element, or nowhere at all, in which case defaults apply. Getting any of the three wrong produces a geometry that is structurally valid and geographically wrong.
The mapping onto shapely is straightforward for everything except curved segments. GML can express an Arc through three points or a Circle by centre and radius; shapely has no curve primitive, so those must be densified into a LineString at a chosen tolerance before they enter the pipeline. Silently dropping them — which is what a parser that matches only on known tags does — loses real geometry.
from __future__ import annotations
import math
from typing import Iterable
from lxml import etree
from shapely.geometry import (LineString, MultiLineString, MultiPolygon,
Point, Polygon)
GML = "http://www.opengis.net/gml/3.2"
GML31 = "http://www.opengis.net/gml"
# Authority definitions whose axis order is latitude, longitude. The short
# "EPSG:4326" form is conventionally longitude-first and is absent here.
LAT_FIRST = {
"urn:ogc:def:crs:EPSG::4326",
"urn:x-ogc:def:crs:EPSG:4326",
"http://www.opengis.net/def/crs/EPSG/0/4326",
"urn:ogc:def:crs:EPSG::4258",
}
GEOMETRY_TAGS = {"Point", "LineString", "LinearRing", "Polygon", "Surface",
"MultiSurface", "MultiCurve", "MultiPoint", "Curve"}
def _local(el: etree._Element) -> str:
return etree.QName(el).localname
def _inherited(el: etree._Element, attribute: str, default: str | None) -> str | None:
"""Walk up until an ancestor declares the attribute — GML allows this."""
node = el
while node is not None:
value = node.get(attribute)
if value:
return value
node = node.getparent()
return default
def find_geometry(feature: etree._Element) -> etree._Element | None:
"""The first GML geometry element at any depth inside a feature."""
for el in feature.iter():
qname = etree.QName(el)
if qname.namespace in (GML, GML31) and qname.localname in GEOMETRY_TAGS:
return el
return None
def positions(text: str, dimension: int, swap: bool) -> list[tuple[float, ...]]:
"""Chunk a coordinate string into positions, fixing axis order.
The multiple assertion is what turns a 3D geometry read as 2D into a
loud failure instead of coordinates silently shifted by one ordinate.
"""
flat = [float(v) for v in text.split()]
if dimension < 2 or len(flat) % dimension:
raise ValueError(
f"{len(flat)} ordinates is not a multiple of srsDimension={dimension}")
out: list[tuple[float, ...]] = []
for i in range(0, len(flat), dimension):
pos = tuple(flat[i:i + dimension])
out.append((pos[1], pos[0], *pos[2:]) if swap else pos)
return out
def _coords(el: etree._Element, dimension: int, swap: bool):
"""Read whichever coordinate carrier this element uses."""
for tag in ("posList", "pos", "coordinates"):
for ns in (GML, GML31):
node = el.find(f".//{{{ns}}}{tag}")
if node is not None and node.text:
text = node.text.replace(",", " ") if tag == "coordinates" else node.text
return positions(text, dimension, swap)
raise ValueError(f"no coordinate carrier inside <{_local(el)}>")
def _rings(polygon: etree._Element, dimension: int, swap: bool):
ns = etree.QName(polygon).namespace
exterior = polygon.find(f"{{{ns}}}exterior")
interiors = polygon.findall(f"{{{ns}}}interior")
if exterior is None: # GML 3.1.1 spelling
exterior = polygon.find(f"{{{ns}}}outerBoundaryIs")
interiors = polygon.findall(f"{{{ns}}}innerBoundaryIs")
shell = _coords(exterior, dimension, swap)
holes = [_coords(node, dimension, swap) for node in interiors]
return shell, holes
def parse_geometry(el: etree._Element, densify_arcs: float = 1.0):
"""Return a shapely geometry, with axis order already corrected."""
srs = _inherited(el, "srsName", None) or ""
swap = srs in LAT_FIRST
dimension = int(_inherited(el, "srsDimension", "2"))
ns = etree.QName(el).namespace
kind = _local(el)
if kind == "Point":
return Point(_coords(el, dimension, swap)[0][:2])
if kind in ("LineString", "LinearRing", "Curve"):
return LineString([p[:2] for p in _coords(el, dimension, swap)])
if kind in ("Polygon", "Surface"):
target = el if kind == "Polygon" else el.find(f".//{{{ns}}}Polygon")
shell, holes = _rings(target, dimension, swap)
return Polygon([p[:2] for p in shell], [[p[:2] for p in h] for h in holes])
if kind == "MultiSurface":
return MultiPolygon([
parse_geometry(poly, densify_arcs)
for poly in el.findall(f".//{{{ns}}}Polygon")])
if kind == "MultiCurve":
return MultiLineString([
parse_geometry(line, densify_arcs)
for line in el.findall(f".//{{{ns}}}LineString")])
raise NotImplementedError(f"unsupported GML geometry <{kind}>")
def densify_arc(start, middle, end, step: float = 1.0) -> list[tuple[float, float]]:
"""Approximate a three-point gml:Arc as a polyline, `step` degrees apart.
shapely has no curve primitive, so an arc that is not densified is an
arc that is silently dropped.
"""
(x1, y1), (x2, y2), (x3, y3) = start, middle, end
d = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
if abs(d) < 1e-12:
return [start, end] # collinear: it is a segment
ux = ((x1 ** 2 + y1 ** 2) * (y2 - y3) + (x2 ** 2 + y2 ** 2) * (y3 - y1)
+ (x3 ** 2 + y3 ** 2) * (y1 - y2)) / d
uy = ((x1 ** 2 + y1 ** 2) * (x3 - x2) + (x2 ** 2 + y2 ** 2) * (x1 - x3)
+ (x3 ** 2 + y3 ** 2) * (x2 - x1)) / d
radius = math.hypot(x1 - ux, y1 - uy)
a0, a1 = math.atan2(y1 - uy, x1 - ux), math.atan2(y3 - uy, x3 - ux)
if a1 < a0:
a1 += 2 * math.pi
steps = max(2, int((a1 - a0) / math.radians(step)))
return [(ux + radius * math.cos(a0 + (a1 - a0) * i / steps),
uy + radius * math.sin(a0 + (a1 - a0) * i / steps))
for i in range(steps + 1)]
_inherited exists because GML allows attribute inheritance. A service is permitted to declare srsName once on the enclosing FeatureCollection and omit it on every geometry inside. A parser that reads the attribute only from the geometry element sees None, falls back to a default, and treats a latitude-first document as longitude-first. Walking up the tree costs nothing and removes the whole class of bug.
The multiple assertion is the single most valuable line. Three-dimensional geometry appears in elevation and bathymetry services more often than expected, and it does not announce itself: a posList of ninety numbers is a valid thirty-position 3D line and an equally valid forty-five-position 2D line. Asserting that the count divides evenly by the declared dimension converts a wrong answer into an exception. Without it, the parser reads the elevation of point one as the latitude of point two and every subsequent coordinate is nonsense that still plots on a map.
Coordinate carriers vary by GML version. GML 3.2 uses posList for sequences and pos for single positions. GML 3.1.1 additionally permits coordinates, a comma-separated form deprecated but still emitted by older MapServer builds. Trying all three in order, and normalising commas to spaces, handles every real response without a version switch.
Exterior and interior rings must be matched on their parent element, not their position. A Polygon writes its shell inside gml:exterior and each hole inside its own gml:interior; the underlying LinearRing elements are identical. Searching for LinearRing at any depth and taking the first as the shell works until a polygon has holes, at which point the holes become separate outlines — the “holes are outlines” branch above.
Assert on longitude by name, never on ordinate position, or the test passes against the bug it exists to catch.
from lxml import etree
from shapely.geometry import Point
DOC = b'''<gml:Point xmlns:gml="http://www.opengis.net/gml/3.2"
srsName="urn:ogc:def:crs:EPSG::4326">
<gml:pos>47.3769 8.5417</gml:pos>
</gml:Point>'''
def test_authority_form_is_latitude_first():
geom = parse_geometry(etree.fromstring(DOC))
assert geom.equals_exact(Point(8.5417, 47.3769), 1e-9)
assert geom.x == 8.5417, "x must be longitude after the swap"
def test_three_dimensional_poslist_is_rejected_when_declared_2d():
doc = etree.fromstring(b'''<gml:LineString
xmlns:gml="http://www.opengis.net/gml/3.2" srsDimension="2">
<gml:posList>0 0 5 1 1 5 2 2 5</gml:posList>
</gml:LineString>''')
# 9 ordinates is not a multiple of 2 — the parser must refuse it.
with pytest.raises(ValueError, match="not a multiple"):
parse_geometry(doc)
Running these against a service you control confirms the two failures that matter most before any of the data reaches a database or a renderer.
The short and long CRS forms mean opposite things. EPSG:4326 is conventionally longitude-first; urn:ogc:def:crs:EPSG::4326 is defined by the authority as latitude-first. Both appear in production GML from mainstream servers. Keeping the latitude-first identifiers in an explicit set — rather than pattern-matching on “4326” — is what keeps the two apart.
Shapely drops the third ordinate on most operations. The parser above keeps only the first two ordinates when constructing geometry, which is deliberate: shapely stores Z but silently discards it through most predicates and transformations, so carrying it further gives a false impression that elevation survives the pipeline. If Z matters, keep it alongside the geometry rather than inside it.
LinearRing must close. GML requires the first and last positions of a ring to be identical, and shapely will close an open ring for you without comment. A ring that arrives open is a sign the source is generating geometry incorrectly, and it is worth asserting rather than accepting.
MultiSurface members can themselves be surfaces. The recursive findall for Polygon above flattens the common case. A MultiSurface whose members are Surface elements containing PolygonPatch needs one more level of unwrapping — rare outside INSPIRE-conformant services, but exactly where it appears, per the ISO 19115 metadata ecosystem.
Both are excellent and both are the right answer when you can take the dependency. GDAL’s GML driver in particular handles curved geometry, schema-driven attribute typing and version differences that hand-written code will not. The reason to write the reader yourself is the streaming case: when a response is large enough that memory matters, controlling the parse loop and the element clearing directly is what keeps the resident set flat.
Fail loudly. Dropping it produces a polygon with a straight edge where the data has a curve, which is a silent geometric error that will propagate into every downstream area calculation and intersection test. Either densify at a documented tolerance or refuse the feature and record why.
No, and that is the trap. It is permitted on the geometry, on a coordinate element, or on an enclosing element, and it defaults to 2 when absent everywhere. Reading it through an inherited lookup and asserting the ordinate count divides evenly is the only combination that is safe against all three placements.
GML permits it, and find_geometry above takes the first, which is usually wrong for such features. When a feature type genuinely carries several geometries — a parcel with a centroid and a boundary, say — read DescribeFeatureType to learn the property names and select by name rather than by document order.
Back to GML and GeoJSON Payload Handling
Related