[WFS provider] Fix handling of LatLongBoundingBox in WFS 1.0

According to the specification, the values of LatLongBoundingBox
are supposed to be in the SRS, not necessarily in WGS84.

But some servers do not follow the spec and return values in WGS84,
so let's try to accomodate for this too.

Fix #14876

Author: rouault

src/providers/wfs/qgswfscapabilities.cpp

@@ -20,6 +20,8 @@
 #include "qgslogger.h"
 #include "qgsmessagelog.h"
 #include "qgsogcutils.h"
+#include "qgscrscache.h"
+
 #include <QDomDocument>
 #include <QSettings>
 #include <QStringList>
@@ -312,31 +314,72 @@ void QgsWFSCapabilities::capabilitiesReplyFinished()
     QDomElement latLongBB = featureTypeElem.firstChildElement( "LatLongBoundingBox" );
     if ( latLongBB.hasAttributes() )
     {
-      featureType.bboxLongLat = QgsRectangle(
-                                  latLongBB.attribute( "minx" ).toDouble(),
-                                  latLongBB.attribute( "miny" ).toDouble(),
-                                  latLongBB.attribute( "maxx" ).toDouble(),
-                                  latLongBB.attribute( "maxy" ).toDouble() );
+      // Despite the name LatLongBoundingBox, the coordinates are supposed to
+      // be expressed in <SRS>. From the WFS schema;
+      // <!-- The LatLongBoundingBox element is used to indicate the edges of
+      // an enclosing rectangle in the SRS of the associated feature type.
+      featureType.bbox = QgsRectangle(
+                           latLongBB.attribute( "minx" ).toDouble(),
+                           latLongBB.attribute( "miny" ).toDouble(),
+                           latLongBB.attribute( "maxx" ).toDouble(),
+                           latLongBB.attribute( "maxy" ).toDouble() );
+      featureType.bboxSRSIsWGS84 = false;
+
+      // But some servers do not honour this and systematically reproject to WGS84
+      // such as GeoServer. See http://osgeo-org.1560.x6.nabble.com/WFS-LatLongBoundingBox-td3813810.html
+      // This is also true of TinyOWS
+      if ( !featureType.crslist.isEmpty() &&
+           featureType.bbox.xMinimum() >= -180 && featureType.bbox.yMinimum() >= -90 &&
+           featureType.bbox.xMaximum() <= 180 && featureType.bbox.yMaximum() < 90 )
+      {
+        QgsCoordinateReferenceSystem crs = QgsCRSCache::instance()->crsByOgcWmsCrs( featureType.crslist[0] );
+        if ( !crs.geographicFlag() )
+        {
+          // If the CRS is projected then check that projecting the corner of the bbox, assumed to be in WGS84,
+          // into the CRS, and then back to WGS84, works (check that we are in the validity area)
+          QgsCoordinateReferenceSystem crsWGS84 = QgsCRSCache::instance()->crsByOgcWmsCrs( "CRS:84" );
+          QgsCoordinateTransform ct( crsWGS84, crs );
+
+          QgsPoint ptMin( featureType.bbox.xMinimum(), featureType.bbox.yMinimum() );
+          QgsPoint ptMinBack( ct.transform( ct.transform( ptMin, QgsCoordinateTransform::ForwardTransform ), QgsCoordinateTransform::ReverseTransform ) );
+          QgsPoint ptMax( featureType.bbox.xMaximum(), featureType.bbox.yMaximum() );
+          QgsPoint ptMaxBack( ct.transform( ct.transform( ptMax, QgsCoordinateTransform::ForwardTransform ), QgsCoordinateTransform::ReverseTransform ) );
+
+          QgsDebugMsg( featureType.bbox.toString() );
+          QgsDebugMsg( ptMinBack.toString() );
+          QgsDebugMsg( ptMaxBack.toString() );
+
+          if ( fabs( featureType.bbox.xMinimum() - ptMinBack.x() ) < 1e-5 &&
+               fabs( featureType.bbox.yMinimum() - ptMinBack.y() ) < 1e-5 &&
+               fabs( featureType.bbox.xMaximum() - ptMaxBack.x() ) < 1e-5 &&
+               fabs( featureType.bbox.yMaximum() - ptMaxBack.y() ) < 1e-5 )
+          {
+            QgsDebugMsg( "Values of LatLongBoundingBox are consistent with WGS84 long/lat bounds, so as the CRS is projected, assume they are indeed in WGS84 and not in the CRS units" );
+            featureType.bboxSRSIsWGS84 = true;
+          }
+        }
+      }
     }
     else
     {
       // WFS 1.1 way
-      latLongBB = featureTypeElem.firstChildElement( "WGS84BoundingBox" );
-      if ( !latLongBB.isNull() )
+      QDomElement WGS84BoundingBox = featureTypeElem.firstChildElement( "WGS84BoundingBox" );
+      if ( !WGS84BoundingBox.isNull() )
       {
-        QDomElement lowerCorner = latLongBB.firstChildElement( "LowerCorner" );
-        QDomElement upperCorner = latLongBB.firstChildElement( "UpperCorner" );
+        QDomElement lowerCorner = WGS84BoundingBox.firstChildElement( "LowerCorner" );
+        QDomElement upperCorner = WGS84BoundingBox.firstChildElement( "UpperCorner" );
         if ( !lowerCorner.isNull() && !upperCorner.isNull() )
         {
           QStringList lowerCornerList = lowerCorner.text().split( " ", QString::SkipEmptyParts );
           QStringList upperCornerList = upperCorner.text().split( " ", QString::SkipEmptyParts );
           if ( lowerCornerList.size() == 2 && upperCornerList.size() == 2 )
           {
-            featureType.bboxLongLat = QgsRectangle(
-                                        lowerCornerList[0].toDouble(),
-                                        lowerCornerList[1].toDouble(),
-                                        upperCornerList[0].toDouble(),
-                                        upperCornerList[1].toDouble() );
+            featureType.bbox = QgsRectangle(
+                                 lowerCornerList[0].toDouble(),
+                                 lowerCornerList[1].toDouble(),
+                                 upperCornerList[0].toDouble(),
+                                 upperCornerList[1].toDouble() );
+            featureType.bboxSRSIsWGS84 = true;
           }
         }
       }

src/providers/wfs/qgswfscapabilities.h

@@ -36,13 +36,14 @@ class QgsWFSCapabilities : public QgsWFSRequest
     struct FeatureType
     {
       //! Default constructor
-      FeatureType() : insertCap( false ), updateCap( false ), deleteCap( false ) {}
+      FeatureType() : bboxSRSIsWGS84( false ), insertCap( false ), updateCap( false ), deleteCap( false ) {}
 
       QString name;
       QString title;
       QString abstract;
       QList<QString> crslist; // first is default
-      QgsRectangle bboxLongLat;
+      QgsRectangle bbox;
+      bool bboxSRSIsWGS84; // if false, the bbox is expressed in crslist[0] CRS
       bool insertCap;
       bool updateCap;
       bool deleteCap;

src/providers/wfs/qgswfsprovider.cpp

@@ -1457,21 +1457,28 @@ bool QgsWFSProvider::getCapabilities()
   {
     if ( thisLayerName == mShared->mCaps.featureTypes[i].name )
     {
-      const QgsRectangle& r = mShared->mCaps.featureTypes[i].bboxLongLat;
+      const QgsRectangle& r = mShared->mCaps.featureTypes[i].bbox;
       if ( mShared->mSourceCRS.authid().isEmpty() && mShared->mCaps.featureTypes[i].crslist.size() != 0 )
       {
         mShared->mSourceCRS = QgsCRSCache::instance()->crsByOgcWmsCrs( mShared->mCaps.featureTypes[i].crslist[0] );
       }
       if ( !r.isNull() )
       {
-        QgsCoordinateReferenceSystem src = QgsCRSCache::instance()->crsByOgcWmsCrs( "CRS:84" );
-        QgsCoordinateTransform ct( src, mShared->mSourceCRS );
+        if ( mShared->mCaps.featureTypes[i].bboxSRSIsWGS84 )
+        {
+          QgsCoordinateReferenceSystem src = QgsCRSCache::instance()->crsByOgcWmsCrs( "CRS:84" );
+          QgsCoordinateTransform ct( src, mShared->mSourceCRS );
 
-        QgsDebugMsg( "latlon ext:" + r.toString() );
-        QgsDebugMsg( "src:" + src.authid() );
-        QgsDebugMsg( "dst:" + mShared->mSourceCRS.authid() );
+          QgsDebugMsg( "latlon ext:" + r.toString() );
+          QgsDebugMsg( "src:" + src.authid() );
+          QgsDebugMsg( "dst:" + mShared->mSourceCRS.authid() );
 
-        mShared->mCapabilityExtent = ct.transformBoundingBox( r, QgsCoordinateTransform::ForwardTransform );
+          mShared->mCapabilityExtent = ct.transformBoundingBox( r, QgsCoordinateTransform::ForwardTransform );
+        }
+        else
+        {
+          mShared->mCapabilityExtent = r;
+        }
 
         QgsDebugMsg( "layer ext:" + mShared->mCapabilityExtent.toString() );
       }

tests/src/python/test_provider_wfs.py

@@ -369,8 +369,9 @@ def testWFS10(self):
       <Name>my:typename</Name>
       <Title>Title</Title>
       <Abstract>Abstract</Abstract>
-      <SRS>EPSG:4326</SRS>
-      <LatLongBoundingBox minx="-71.123" miny="66.33" maxx="-65.32" maxy="78.3"/>
+      <SRS>EPSG:32631</SRS>
+      <!-- in WFS 1.0, LatLongBoundingBox is in SRS units, not necessarily lat/long... -->
+      <LatLongBoundingBox minx="400000" miny="5400000" maxx="450000" maxy="5500000"/>
     </FeatureType>
   </FeatureTypeList>
 </WFS_Capabilities>""".encode('UTF-8'))
@@ -402,11 +403,11 @@ def testWFS10(self):
         self.assertEqual(vl.wkbType(), QgsWKBTypes.Point)
         self.assertEqual(len(vl.fields()), 5)
         self.assertEqual(vl.featureCount(), 0)
-        reference = QgsGeometry.fromRect(QgsRectangle(-71.123, 66.33, -65.32, 78.3))
+        reference = QgsGeometry.fromRect(QgsRectangle(400000.0, 5400000.0, 450000.0, 5500000.0))
         vl_extent = QgsGeometry.fromRect(vl.extent())
         assert QgsGeometry.compare(vl_extent.asPolygon()[0], reference.asPolygon()[0], 0.00001), 'Expected {}, got {}'.format(reference.exportToWkt(), vl_extent.exportToWkt())
 
-        with open(sanitize(endpoint, '?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.0.0&TYPENAME=my:typename&SRSNAME=EPSG:4326'), 'wb') as f:
+        with open(sanitize(endpoint, '?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.0.0&TYPENAME=my:typename&SRSNAME=EPSG:32631'), 'wb') as f:
             f.write("""
 <wfs:FeatureCollection
                        xmlns:wfs="http://www.opengis.net/wfs"
@@ -416,7 +417,7 @@ def testWFS10(self):
   <gml:featureMember>
     <my:typename fid="typename.0">
       <my:geometryProperty>
-          <gml:Point srsName="http://www.opengis.net/gml/srs/epsg.xml#4326"><gml:coordinates decimal="." cs="," ts=" ">2,49</gml:coordinates></gml:Point></my:geometryProperty>
+          <gml:Point srsName="http://www.opengis.net/gml/srs/epsg.xml#4326"><gml:coordinates decimal="." cs="," ts=" ">426858,5427937</gml:coordinates></gml:Point></my:geometryProperty>
       <my:INTFIELD>1</my:INTFIELD>
       <my:GEOMETRY>2</my:GEOMETRY>
       <my:longfield>1234567890123</my:longfield>
@@ -448,7 +449,7 @@ def testWFS10(self):
 
         got_f = [f for f in vl.getFeatures()]
         got = got_f[0].geometry().geometry()
-        self.assertEqual((got.x(), got.y()), (2.0, 49.0))
+        self.assertEqual((got.x(), got.y()), (426858.0, 5427937.0))
 
         self.assertEqual(vl.featureCount(), 1)
 
@@ -459,6 +460,51 @@ def testWFS10(self):
 
         assert not vl.dataProvider().deleteFeatures([0])
 
+    def testWFS10_latlongboundingbox_in_WGS84(self):
+        """Test WFS 1.0 with non conformatn LatLongBoundingBox"""
+
+        endpoint = self.__class__.basetestpath + '/fake_qgis_http_endpoint_WFS1.0_latlongboundingbox_in_WGS84'
+
+        with open(sanitize(endpoint, '?SERVICE=WFS?REQUEST=GetCapabilities?VERSION=1.0.0'), 'wb') as f:
+            f.write("""
+<WFS_Capabilities version="1.0.0" xmlns="http://www.opengis.net/wfs" xmlns:ogc="http://www.opengis.net/ogc">
+  <FeatureTypeList>
+    <FeatureType>
+      <Name>my:typename</Name>
+      <Title>Title</Title>
+      <Abstract>Abstract</Abstract>
+      <SRS>EPSG:32631</SRS>
+      <!-- in WFS 1.0, LatLongBoundingBox are supposed to be in SRS units, not necessarily lat/long...
+        But some servers do not honour this, so let's try to be robust -->
+      <LatLongBoundingBox minx="1.63972075372399" miny="48.7449841112119" maxx="2.30733562794991" maxy="49.6504711179582"/>
+    </FeatureType>
+  </FeatureTypeList>
+</WFS_Capabilities>""".encode('UTF-8'))
+
+        with open(sanitize(endpoint, '?SERVICE=WFS&REQUEST=DescribeFeatureType&VERSION=1.0.0&TYPENAME=my:typename'), 'wb') as f:
+            f.write("""
+<xsd:schema xmlns:my="http://my" xmlns:gml="http://www.opengis.net/gml" xmlns:xsd="http://www.w3.org/2001/XMLSchema" elementFormDefault="qualified" targetNamespace="http://my">
+  <xsd:import namespace="http://www.opengis.net/gml"/>
+  <xsd:complexType name="typenameType">
+    <xsd:complexContent>
+      <xsd:extension base="gml:AbstractFeatureType">
+        <xsd:sequence>
+          <xsd:element maxOccurs="1" minOccurs="0" name="geometryProperty" nillable="true" type="gml:PointPropertyType"/>
+        </xsd:sequence>
+      </xsd:extension>
+    </xsd:complexContent>
+  </xsd:complexType>
+  <xsd:element name="typename" substitutionGroup="gml:_Feature" type="my:typenameType"/>
+</xsd:schema>
+""".encode('UTF-8'))
+
+        vl = QgsVectorLayer(u"url='http://" + endpoint + u"' typename='my:typename' version='1.0.0'", u'test', u'WFS')
+        assert vl.isValid()
+
+        reference = QgsGeometry.fromRect(QgsRectangle(399999.9999999680439942,5399338.9090830031782389,449999.9999999987776391,5500658.0448500607162714))
+        vl_extent = QgsGeometry.fromRect(vl.extent())
+        assert QgsGeometry.compare(vl_extent.asPolygon()[0], reference.asPolygon()[0], 0.00001), 'Expected {}, got {}'.format(reference.exportToWkt(), vl_extent.exportToWkt())
+
     def testWFST10(self):
         """Test WFS-T 1.0 (read-write)"""
 
@@ -1959,14 +2005,19 @@ def testGeomedia(self):
         assert vl.isValid()
         self.assertEqual(vl.wkbType(), QgsWKBTypes.MultiPolygon)
 
+        # Extent before downloading features
+        reference = QgsGeometry.fromRect(QgsRectangle(243900.3520259926444851,4427769.1559739429503679,1525592.3040170343592763,5607994.6020106188952923))
+        vl_extent = QgsGeometry.fromRect(vl.extent())
+        assert QgsGeometry.compare(vl_extent.asPolygon()[0], reference.asPolygon()[0], 0.00001), 'Expected {}, got {}'.format(reference.exportToWkt(), vl_extent.exportToWkt())
+
         # Download all features
         features = [f for f in vl.getFeatures()]
         self.assertEqual(len(features), 2)
 
         reference = QgsGeometry.fromRect(QgsRectangle(500000, 4500000, 510000, 4510000))
         vl_extent = QgsGeometry.fromRect(vl.extent())
-        self.assertEqual(features[0]['intfield'], 1)
         assert QgsGeometry.compare(vl_extent.asPolygon()[0], reference.asPolygon()[0], 0.00001), 'Expected {}, got {}'.format(reference.exportToWkt(), vl_extent.exportToWkt())
+        self.assertEqual(features[0]['intfield'], 1)
         self.assertEqual(features[1]['intfield'], 2)