Reading GML features with shapely in Python?

Question

Is it possible to read GML features directly with shapely in Python?

We are able to read a GML string with GDAL to shapely via below code, but we prefer to avoid GDAL and use shapely directly:

from osgeo import ogr
import json
from shapely.geometry import shape

geom_str = '<gml:Polygon srsDimension="3" srsName="urn:ogc:def:crs:EPSG::28992"><gml:exterior><gml:LinearRing><gml:posList count="10">207819.658 373249.741 0.0 207817.964 373252.252 0.0 207811.797 373248.147 0.0 207813.477 373245.626 0.0 207814.879 373243.523 0.0 207816.619 373240.965 0.0 207818.809 373242.423 0.0 207826.753 373247.712 0.0 207823.637 373252.39 0.0 207819.658 373249.741 0.0</gml:posList></gml:LinearRing></gml:exterior></gml:Polygon>'
geom_shape = ogr.CreateGeometryFromGML(geom_str)
shape(json.loads(geom_shape.ExportToJson()))

Answer 1

You can use any Python library that parse xml files (see Using Python to parse an XML containing GML tags or How to extract data from GML file for example)

With pygml that provides the Geo Interface protocol

import pygml
geom = pygml.parse("""<gml:Polygon srsDimension="3" srsName="urn:ogc:def:crs:EPSG::28992" xmlns:gml="http://www.opengis.net/gml/3.2"><gml:exterior><gml:LinearRing><gml:posList count="10">207819.658 373249.741 0.0 207817.964 373252.252 0.0 207811.797 373248.147 0.0 207813.477 373245.626 0.0 207814.879 373243.523 0.0 207816.619 373240.965 0.0 207818.809 373242.423 0.0 207826.753 373247.712 0.0 207823.637 373252.39 0.0 207819.658 373249.741 0.0</gml:posList></gml:LinearRing></gml:exterior></gml:Polygon>""")
print(geom.__geo_interface__)
{'type': 'Polygon', 'coordinates': [[(207819.658, 373249.741), (0.0, 207817.964), (373252.252, 0.0), (207811.797, 373248.147), (0.0, 207813.477), (373245.626, 0.0), (207814.879, 373243.523), (0.0, 207816.619), (373240.965, 0.0), (207818.809, 373242.423), (0.0, 207826.753), (373247.712, 0.0), (207823.637, 373252.39), (0.0, 207819.658), (373249.741, 0.0)]], 'crs': {'type': 'name', 'properties': {'name': 'urn:ogc:def:crs:EPSG::28992'}}}


from shapely.geometry import shape
geom_shape=  shape(geom.__geo_interface__)
print(type(geom_shape))
<class 'shapely.geometry.polygon.Polygon'>
print(geom_shape)
POLYGON ((207819.658 373249.741, 0 207817.964, 373252.252 0, 207811.797 373248.147, 0 207813.477, 373245.626 0, 207814.879 373243.523, 0 207816.619, 373240.965 0, 207818.809 373242.423, 0 207826.753, 373247.712 0, 207823.637 373252.39, 0 207819.658, 373249.741 0, 207819.658 373249.741))