[GIS] How to get GML layer CRS with OGR Python

gmlogrpythonwfs

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?

Best Answer

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'
Related Question