1

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?

  • could you query the first feature or the bounding box as the seems to be where it is declared in the GML – Ian Turton Feb 5 '15 at 10:43
  • Hi Ian. The concept of bounding box does not seem to exist in OGR Python. As for the Feature class, there are no CRS related methods available. – Luís de Sousa Feb 5 '15 at 10:47
1

After some digging I find 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'

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.