How to get GML layer CRS with OGR Python?

Question

I have a GML layer obtained from a WFS that looks like this:

<wfs:FeatureCollection xmlns:ms="http://mapserver.gis.umn.edu/mapserver" xmlns:wfs="http://www.opengis.net/wfs" xmlns:gml="http://www.opengis.net/gml" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.0.0/WFS-basic.xsd http://mapserver.gis.umn.edu/mapserver http://maps.iguess.list.lu/cgi-bin/mapserv?map=/srv/mapserv/MapFiles/LUX_samples.map&SERVICE=WFS&VERSION=1.0.0&REQUEST=DescribeFeatureType&TYPENAME=PointsOfInterest2169&OUTPUTFORMAT=XMLSCHEMA">
<gml:boundedBy>
<gml:Box srsName="EPSG:2169">
<gml:coordinates>
60252.601105,62852.917507 88920.374823,102462.589590
</gml:coordinates>
</gml:Box>
</gml:boundedBy>
<gml:featureMember>
<ms:PointsOfInterest2169>
<gml:boundedBy>
<gml:Box srsName="EPSG:2169">
<gml:coordinates>
71194.499471,102462.589590 71194.499471,102462.589590
</gml:coordinates>
</gml:Box>
</gml:boundedBy>
<ms:msGeometry>
<gml:Point srsName="EPSG:2169">
<gml:coordinates>71194.499471,102462.589590</gml:coordinates>
</gml:Point>
</ms:msGeometry>
<ms:fid>LuxPOIs.0</ms:fid>
<ms:id>7</ms:id>
</ms:PointsOfInterest2169>
</gml:featureMember>
<gml:featureMember>
<ms:PointsOfInterest2169>
<gml:boundedBy>
<gml:Box srsName="EPSG:2169">
<gml:coordinates>
60252.601105,92833.719028 60252.601105,92833.719028
</gml:coordinates>
</gml:Box>
</gml:boundedBy>
<ms:msGeometry>
<gml:Point srsName="EPSG:2169">
<gml:coordinates>60252.601105,92833.719028</gml:coordinates>
</gml:Point>
</ms:msGeometry>
<ms:fid>LuxPOIs.1</ms:fid>
<ms:id>6</ms:id>
</ms:PointsOfInterest2169>
</gml:featureMember>
</wfs:FeatureCollection>

I need to obtain the CRS of this layer in a PyWPS process, but I can not figure how. I saved this GML bit into a file (test.gml) and made a few tests:

$ ogrinfo test.gml
Had to open data source read-only.
INFO: Open of `test.gml'
      using driver `GML' successful.
1: PointsOfInterest2169

$ python
Python 2.7.6 (default, Mar 22 2014, 22:59:56) 
[GCC 4.8.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from osgeo import ogr
>>> inSource = ogr.Open("test.gml")
>>> inLayer = inSource.GetLayer()
>>> inLayer.GetSpatialRef()
>>> 

The CRS is always reported as null, even though the information is there in the GML. What would be the correct way of obtaining it with OGR Python?

Answer 1

After some digging I found out a GetSpatialReference method for the Geometry class (note the different name to the GetSpatialRef method in the Layer class). So the code should follow something like:

$ python
Python 2.7.6 (default, Mar 22 2014, 22:59:56) 
[GCC 4.8.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from osgeo import ogr
>>> inSource = ogr.Open("test.gml")
>>> inLayer = inSource.GetLayer()
>>> f = inLayer.GetNextFeature()
>>> g = f.GetGeometryRef()
>>> r = g.GetSpatialReference()
>>> str(r.GetAuthorityCode("PROJCS"))
'2169'