Stream the GML members, normalise each geometry into shapely with longitude first, reproject anything that is not already CRS84, flatten nested properties with a documented separator, and emit newline-delimited GeoJSON. The conversion is lossy by specification — the job is to decide deliberately where each lost construct goes rather than discovering the loss downstream.
GML is a schema-driven, reference-system-aware, arbitrarily nested encoding. RFC 7946 GeoJSON is a fixed-CRS format with a flat properties object and an optional identifier. The gap is not a matter of parser effort; three constructs have no conformant destination.
The gml:id attribute is mandatory on every GML feature and is what a WFS-T Update or Delete addresses. GeoJSON’s id member is optional and, in practice, ignored by a great many consumers — so writing it in both places is the only reliable way to keep round-tripping possible.
The srsName on a GML geometry can name any authority definition. RFC 7946 removed the crs member entirely and defines all coordinates as CRS84. A conformant conversion therefore has to reproject, and the original identifier has to be recorded somewhere outside the geometry if it will ever be needed again.
Nested properties are the quiet one. GML permits a property that contains a complex type with its own children; GeoJSON’s properties is a flat object. Flattening is unavoidable, and the only question is whether the flattening is deterministic and documented or accidental.
from __future__ import annotations
import json
from typing import Any, Iterator, TextIO
import requests
from lxml import etree
from pyproj import CRS, Transformer
from shapely.geometry import mapping
from shapely.ops import transform as shapely_transform
WFS20 = "http://www.opengis.net/wfs/2.0"
GML = "http://www.opengis.net/gml/3.2"
CRS84 = CRS.from_epsg(4326)
FLATTEN_SEPARATOR = "."
class ConversionReport:
"""What the conversion had to change, so it is never a surprise later."""
def __init__(self) -> None:
self.features = 0
self.reprojected_from: set[str] = set()
self.flattened_keys: set[str] = set()
self.missing_ids = 0
def as_dict(self) -> dict[str, Any]:
return {
"features": self.features,
"reprojected_from": sorted(self.reprojected_from),
"flattened_keys": sorted(self.flattened_keys),
"features_without_gml_id": self.missing_ids,
}
def _transformer_for(srs: str) -> Transformer | None:
"""None when the source is already CRS84 — skip the work entirely."""
source = CRS.from_user_input(srs)
if source.equals(CRS84):
return None
return Transformer.from_crs(source, CRS84, always_xy=True)
def _flatten(node: etree._Element, prefix: str, out: dict[str, Any],
report: ConversionReport) -> None:
"""Flatten a GML property subtree into dotted keys.
A collision means two different GML paths produced the same flat key,
which would silently discard one of them — so it is an error, not a
last-write-wins.
"""
for child in node:
name = etree.QName(child).localname
key = f"{prefix}{FLATTEN_SEPARATOR}{name}" if prefix else name
if len(child):
report.flattened_keys.add(key)
_flatten(child, key, out, report)
else:
if key in out:
raise ValueError(f"flatten collision on {key!r}")
out[key] = (child.text or "").strip() or None
def convert(url: str, params: dict[str, str], sink: TextIO,
parse_geometry, report: ConversionReport | None = None) -> ConversionReport:
"""Stream a WFS GML response into newline-delimited GeoJSON.
`parse_geometry` is the reader from the GML geometry guide: it returns a
shapely object with axis order already corrected.
"""
report = report or ConversionReport()
with requests.get(url, params=params, stream=True, timeout=300) as resp:
resp.raise_for_status()
resp.raw.decode_content = True
context = etree.iterparse(resp.raw, events=("end",), tag=f"{{{WFS20}}}member")
for _, member in context:
for feature_node in member:
sink.write(json.dumps(
_to_geojson(feature_node, parse_geometry, report)) + "\n")
report.features += 1
member.clear()
while member.getprevious() is not None:
del member.getparent()[0]
return report
def _to_geojson(node: etree._Element, parse_geometry,
report: ConversionReport) -> dict[str, Any]:
fid = node.get(f"{{{GML}}}id")
if not fid:
report.missing_ids += 1
geometry = None
geom_node = _first_geometry(node)
if geom_node is not None:
geom = parse_geometry(geom_node)
srs = geom_node.get("srsName") or "EPSG:4326"
transformer = _transformer_for(srs)
if transformer is not None:
report.reprojected_from.add(srs)
geom = shapely_transform(transformer.transform, geom)
geometry = mapping(geom)
properties: dict[str, Any] = {}
for child in node:
if etree.QName(child).namespace == GML or child is geom_node:
continue
if _first_geometry(child) is not None:
continue # the geometry property itself
_flatten(child, "", properties, report)
# gml:id lives in both places: the id member for conformant consumers,
# and properties for the many that ignore it.
if fid:
properties.setdefault("gml_id", fid)
out: dict[str, Any] = {"type": "Feature", "geometry": geometry,
"properties": properties}
if fid:
out["id"] = fid
return out
def _first_geometry(node: etree._Element) -> etree._Element | None:
for el in node.iter():
q = etree.QName(el)
if q.namespace == GML and q.localname in {
"Point", "LineString", "Polygon", "MultiSurface", "MultiCurve"}:
return el
return None
if __name__ == "__main__":
from gml_geometry import parse_geometry # from the geometry guide
with open("parcels.ndjson", "w", encoding="utf-8") as fh:
summary = convert(
"https://example.org/geoserver/wfs",
{"SERVICE": "WFS", "VERSION": "2.0.0", "REQUEST": "GetFeature",
"typeNames": "cadastre:parcels", "count": "50000"},
fh, parse_geometry)
print(json.dumps(summary.as_dict(), indent=2))
The report is the deliverable, not a log line. A conversion that quietly reprojects from EPSG:25832 and flattens four nested properties has changed the data in ways a consumer three systems away will eventually care about. Returning a structured summary — which reference systems were reprojected from, which keys were flattened, how many features arrived without an identifier — makes that visible at the point of conversion rather than at the point of confusion.
Reprojection is skipped when it is a no-op. _transformer_for returns None when the source system already equals CRS84, so a service already publishing in WGS 84 pays nothing. Building a Transformer costs tens of milliseconds, which is irrelevant once but significant across a million features, which is why the transformer is built per distinct source system rather than per geometry.
Flatten collisions are errors. Two GML paths that flatten to the same dotted key mean one value would overwrite the other. Silently keeping the last is the behaviour of most quick conversion scripts and it loses data invisibly. Raising names the key and stops the run, which is recoverable; discovering it later, from a consumer, is not.
Newline-delimited output is append-only. Writing a single JSON array means a run that dies at ninety per cent leaves a file with no closing bracket, which no parser will read. One feature per line means a partial file is a valid partial result, the run can be resumed from the last line, and downstream consumers can start before the conversion finishes.
Convert a small extract and check the three lossy constructs explicitly:
python convert_gml.py && head -n 1 parcels.ndjson | python -m json.tool
{
"type": "Feature",
"id": "parcels.4471",
"geometry": {
"type": "Polygon",
"coordinates": [[[8.5401, 47.3752], [8.5419, 47.3752], [8.5419, 47.3768]]]
},
"properties": {
"gml_id": "parcels.4471",
"owner.name": "Kanton Zurich",
"owner.type": "public",
"area_m2": "1284.5"
}
}
The summary reports what changed:
{
"features": 50000,
"reprojected_from": ["urn:ogc:def:crs:EPSG::2056"],
"flattened_keys": ["owner"],
"features_without_gml_id": 0
}
Assert on longitude, not on the first ordinate — a test written as coordinates[0][0][0] == 47.37 passes against exactly the axis bug it should catch.
Property values arrive as strings. GML carries type information in the schema, not in the instance document, so a flattened property is text until something casts it. Reading DescribeFeatureType and applying the declared types during conversion is the correct fix; the shortcut of guessing from the value produces a field that is an integer in one file and a string in the next.
Reprojecting changes coordinate precision. A transform from a projected system into degrees produces long decimal expansions that inflate the output substantially. Rounding to a documented precision — six decimal places is roughly ten centimetres — is worth doing deliberately, because the alternative is a file that is larger than it needs to be and implies accuracy the data does not have.
Not every GML feature has a geometry. Metadata-only records are legal and common. Emitting "geometry": null is the conformant representation and is different from an empty geometry, which is what a converter that constructs an empty Polygon produces.
Check what the service can do before converting anything. Many WFS deployments offer outputFormat=application/json directly, which makes the whole conversion unnecessary. Reading the advertised output formats from GetCapabilities first — as OGC API - Features and the REST Transition discusses — is the cheapest possible optimisation.
Yes, to a documented precision. Reprojecting into degrees produces fifteen significant figures that imply nanometre accuracy no survey supports, and the extra digits are pure payload. Six decimal places is about ten centimetres at the equator and is a sensible default; state the choice in the conversion report so a consumer knows the precision is deliberate rather than accidental.
Because consumers disagree about which they read. The top-level id member is the conformant home and is what a well-behaved client uses, but a large number of tools — including several popular GIS desktop importers — drop it and surface only properties. Writing both costs a few bytes per feature and is the difference between a round trip that works and one that cannot address the original records.
Only if collisions are treated as errors, as the code above does. A dotted key is also worth checking against the consumer: some document stores interpret dots as path separators, which turns a flat key back into nesting on write. Where that is a risk, use a double underscore instead and document it in the same place the separator is defined.
Yes, but it needs information GeoJSON does not carry. Producing valid GML requires a target feature type with a schema, a namespace, an element name per property and a feature identifier — none of which are in the GeoJSON. In practice the source of truth for a GeoJSON-to-GML conversion is the DescribeFeatureType response of the destination service, not the GeoJSON itself.
Back to GML and GeoJSON Payload Handling
Related