[gdal-dev] ExportToGML with compound CRS

[Joep Keijsers]

Greetings,

I'm trying to export a geometry object to GML using the Python API. The geometry has a compound CRS (EPSG:5973). When exporting to GML, the vertical component of the compound CRS is dropped, although the Z-values are still there. For example:

from osgeo import ogr, osr

# Create geometry from WKT string
srs = osr.SpatialReference()
srs.ImportFromEPSG(5973)
wkt = "POINT (30 10 5)"
geometry = ogr.CreateGeometryFromWkt(wkt, srs)

# Export to GML
xml = geometry.ExportToGML(options=['FORMAT=GML32', 'SRSDIMENSION_LOC=GEOMETRY', 'SRSNAME_FORMAT=SHORT', 'GMLID=TestPoint-1'])
print(xml)

... which returns EPSG:11023 instead of EPSG:5973:

<gml:Point srsName="EPSG:11023" srsDimension="3" gml:id="TestPoint-1"><gml:pos>30 10 5</gml:pos></gml:Point>

I believe this EPSG code change happens in this line: https://github.com/OSGeo/gdal/blob/30e7918f245022983d84f93f1a53641fd50337ed/ogr/ogr2gmlgeometry.cpp#L562.

Is there a way to keep the compound CRS in the GML output?

Version info: gdal 3.12.3 

Thanks!

Joep