Allow specifying a destination CRS in QgsFeatureRequest

If set, all geometries will be reprojected from their original
coordinate reference system to the destination CRS while
iterating over features.

If a CRS has been set as the destination CRS, then the filterRect
parameter should be specified in the same CRS as this destination
CRS.

Additionally, a callback function can be specified on the request
to be called if a transform exception is encountered while
iterating over features.

This is designed to make it easier for scripts and plugins to
correctly reproject layers in an efficient and robust way, instead
of having to implement lots of repeated code themselves and
potentially missing some of the important considerations which
come with reprojecting geometries & bounding boxes.

Now, if a script wants the features from a layer in a specific
CRS, they can call:

    crs = QgsCoordinateReferenceSystem('epsg:4326')
    request = QgsFeatureRequest().setDestinationCrs(crs)
    for f in layer.getFeatures(reqeuest):
        print('geometry in 4326 is {}.format(f.geometry().exportToWkt()))

Author: nyalldawson

python/core/qgsfeatureiterator.sip

@@ -102,6 +102,28 @@ end of iterating: free the resources / lock
  :rtype: bool
 %End
 
+    void transformFeatureGeometry( QgsFeature &feature, const QgsCoordinateTransform &transform ) const;
+%Docstring
+ Transforms ``feature``'s geometry according to the specified coordinate ``transform``.
+ If ``feature`` has no geometry or ``transform`` is invalid then calling this method
+ has no effect and will be shortcut.
+ Iterators should call this method before returning features to ensure that any
+ QgsFeatureRequest.destinationCrs() set on the request is respected.
+.. versionadded:: 3.0
+%End
+
+
+    QgsRectangle transformedFilterRect( const QgsCoordinateTransform &transform ) const;
+%Docstring
+ Returns a rectangle representing the original request's QgsFeatureRequest.filterRect().
+ If ``transform`` is a valid coordinate transform, the return rectangle will represent
+ the requested filterRect() transformed to the source's coordinate reference system.
+ Iterators should call this method and use the returned rectangle for filtering
+ features to ensure that any QgsFeatureRequest.destinationCrs() set on the request is respected.
+.. versionadded:: 3.0
+ :rtype: QgsRectangle
+%End
+
 
 
 

python/core/qgsfeaturerequest.sip

@@ -267,10 +267,16 @@ construct a request with feature ID filter
 %Docstring
 construct a request with feature ID filter
 %End
-    explicit QgsFeatureRequest( const QgsRectangle &rect );
+
+    explicit QgsFeatureRequest( const QgsRectangle &rectangle );
 %Docstring
-construct a request with rectangle filter
+ Construct a request with ``rectangle`` bounding box filter.
+
+ When a destination CRS is set using setDestinationCrs(), ``rectangle``
+ is expected to be in the same CRS as the destinationCrs(). Otherwise, ``rectangle``
+ should use the same CRS as the source layer/provider.
 %End
+
     explicit QgsFeatureRequest( const QgsExpression &expr, const QgsExpressionContext &context = QgsExpressionContext() );
 %Docstring
 construct a request with a filter expression
@@ -288,16 +294,28 @@ copy constructor
  :rtype: FilterType
 %End
 
-    QgsFeatureRequest &setFilterRect( const QgsRectangle &rect );
+    QgsFeatureRequest &setFilterRect( const QgsRectangle &rectangle );
 %Docstring
- Set rectangle from which features will be taken. Empty rectangle removes the filter.
+ Sets the ``rectangle`` from which features will be taken. An empty rectangle removes the filter.
+
+ When a destination CRS is set using setDestinationCrs(), ``rectangle``
+ is expected to be in the same CRS as the destinationCrs(). Otherwise, ``rectangle``
+ should use the same CRS as the source layer/provider.
+
+.. seealso:: filterRect()
  :rtype: QgsFeatureRequest
 %End
 
     const QgsRectangle &filterRect() const;
 %Docstring
- Get the rectangle from which features will be taken. If the returned
+ Returns the rectangle from which features will be taken. If the returned
  rectangle is null, then no filter rectangle is set.
+
+ When a destination CRS is set using setDestinationCrs(), the rectangle
+ will be in the same CRS as the destinationCrs(). Otherwise, the rectangle
+ will use the same CRS as the source layer/provider.
+
+.. seealso:: setFilterRect()
  :rtype: QgsRectangle
 %End
 
@@ -528,6 +546,60 @@ Set a subset of attributes by names that will be fetched
  :rtype: QgsSimplifyMethod
 %End
 
+    QgsCoordinateReferenceSystem destinationCrs() const;
+%Docstring
+ Returns the destination coordinate reference system for feature's geometries,
+ or an invalid QgsCoordinateReferenceSystem if no reprojection will be done
+ and all features will be left with their original geometry.
+.. seealso:: setDestinationCrs()
+.. versionadded:: 3.0
+ :rtype: QgsCoordinateReferenceSystem
+%End
+
+    QgsFeatureRequest &setDestinationCrs( const QgsCoordinateReferenceSystem &crs );
+%Docstring
+ Sets the destination ``crs`` for feature's geometries. If set, all
+ geometries will be reprojected from their original coordinate reference
+ system to this desired reference system. If ``crs`` is an invalid
+ QgsCoordinateReferenceSystem then no reprojection will be done
+ and all features will be left with their original geometry.
+
+ When a ``crs`` is set using setDestinationCrs(), then any filterRect()
+ set on the request is expected to be in the same CRS as the destination
+ CRS.
+
+.. seealso:: destinationCrs()
+.. versionadded:: 3.0
+ :rtype: QgsFeatureRequest
+%End
+
+    QgsFeatureRequest &setTransformErrorCallback( SIP_PYCALLABLE / AllowNone / );
+%Docstring
+ Sets a callback function to use when encountering a transform error when iterating
+ features and a destinationCrs() is set. This function will be
+ called using the feature which encountered the transform error as a parameter.
+.. versionadded:: 3.0
+.. seealso:: transformErrorCallback()
+.. seealso:: setDestinationCrs()
+ :rtype: QgsFeatureRequest
+%End
+%MethodCode
+    Py_BEGIN_ALLOW_THREADS
+
+    sipCpp->setTransformErrorCallback( [a0]( const QgsFeature &arg )
+    {
+      SIP_BLOCK_THREADS
+      Py_XDECREF( sipCallMethod( NULL, a0, "D", &arg, sipType_QgsFeature, NULL ) );
+      SIP_UNBLOCK_THREADS
+    } );
+
+    sipRes = sipCpp;
+
+    Py_END_ALLOW_THREADS
+%End
+
+
+
     bool acceptFeature( const QgsFeature &feature );
 %Docstring
  Check if a feature is accepted by this requests filter

src/core/providers/memory/qgsmemoryfeatureiterator.cpp

@@ -27,25 +27,31 @@
 QgsMemoryFeatureIterator::QgsMemoryFeatureIterator( QgsMemoryFeatureSource *source, bool ownSource, const QgsFeatureRequest &request )
   : QgsAbstractFeatureIteratorFromSource<QgsMemoryFeatureSource>( source, ownSource, request )
 {
+  if ( mRequest.destinationCrs().isValid() && mRequest.destinationCrs() != mSource->mCrs )
+  {
+    mTransform = QgsCoordinateTransform( mSource->mCrs, mRequest.destinationCrs() );
+  }
+  mFilterRect = transformedFilterRect( mTransform );
+
   if ( !mSource->mSubsetString.isEmpty() )
   {
     mSubsetExpression = new QgsExpression( mSource->mSubsetString );
     mSubsetExpression->prepare( &mSource->mExpressionContext );
   }
 
-  if ( !mRequest.filterRect().isNull() && mRequest.flags() & QgsFeatureRequest::ExactIntersect )
+  if ( !mFilterRect.isNull() && mRequest.flags() & QgsFeatureRequest::ExactIntersect )
   {
-    mSelectRectGeom = QgsGeometry::fromRect( request.filterRect() );
+    mSelectRectGeom = QgsGeometry::fromRect( mFilterRect );
     mSelectRectEngine.reset( QgsGeometry::createGeometryEngine( mSelectRectGeom.geometry() ) );
     mSelectRectEngine->prepareGeometry();
   }
 
   // if there's spatial index, use it!
   // (but don't use it when selection rect is not specified)
-  if ( !mRequest.filterRect().isNull() && mSource->mSpatialIndex )
+  if ( !mFilterRect.isNull() && mSource->mSpatialIndex )
   {
     mUsingFeatureIdList = true;
-    mFeatureIdList = mSource->mSpatialIndex->intersects( mRequest.filterRect() );
+    mFeatureIdList = mSource->mSpatialIndex->intersects( mFilterRect );
     QgsDebugMsg( "Features returned by spatial index: " + QString::number( mFeatureIdList.count() ) );
   }
   else if ( mRequest.filterType() == QgsFeatureRequest::FilterFid )
@@ -92,7 +98,7 @@ bool QgsMemoryFeatureIterator::nextFeatureUsingList( QgsFeature &feature )
   // option 1: we have a list of features to traverse
   while ( mFeatureIdListIterator != mFeatureIdList.constEnd() )
   {
-    if ( !mRequest.filterRect().isNull() && mRequest.flags() & QgsFeatureRequest::ExactIntersect )
+    if ( !mFilterRect.isNull() && mRequest.flags() & QgsFeatureRequest::ExactIntersect )
     {
       // do exact check in case we're doing intersection
       if ( mSource->mFeatures.value( *mFeatureIdListIterator ).hasGeometry() && mSelectRectEngine->intersects( *mSource->mFeatures.value( *mFeatureIdListIterator ).geometry().geometry() ) )
@@ -124,7 +130,10 @@ bool QgsMemoryFeatureIterator::nextFeatureUsingList( QgsFeature &feature )
     close();
 
   if ( hasFeature )
+  {
     feature.setFields( mSource->mFields ); // allow name-based attribute lookups
+    transformFeatureGeometry( feature, mTransform );
+  }
 
   return hasFeature;
 }
@@ -137,7 +146,7 @@ bool QgsMemoryFeatureIterator::nextFeatureTraverseAll( QgsFeature &feature )
   // option 2: traversing the whole layer
   while ( mSelectIterator != mSource->mFeatures.constEnd() )
   {
-    if ( mRequest.filterRect().isNull() )
+    if ( mFilterRect.isNull() )
     {
       // selection rect empty => using all features
       hasFeature = true;
@@ -153,7 +162,7 @@ bool QgsMemoryFeatureIterator::nextFeatureTraverseAll( QgsFeature &feature )
       else
       {
         // check just bounding box against rect when not using intersection
-        if ( mSelectIterator->hasGeometry() && mSelectIterator->geometry().boundingBox().intersects( mRequest.filterRect() ) )
+        if ( mSelectIterator->hasGeometry() && mSelectIterator->geometry().boundingBox().intersects( mFilterRect ) )
           hasFeature = true;
       }
     }
@@ -178,6 +187,7 @@ bool QgsMemoryFeatureIterator::nextFeatureTraverseAll( QgsFeature &feature )
     ++mSelectIterator;
     feature.setValid( true );
     feature.setFields( mSource->mFields ); // allow name-based attribute lookups
+    transformFeatureGeometry( feature, mTransform );
   }
   else
     close();
@@ -216,6 +226,7 @@ QgsMemoryFeatureSource::QgsMemoryFeatureSource( const QgsMemoryProvider *p )
   , mFeatures( p->mFeatures )
   , mSpatialIndex( p->mSpatialIndex ? new QgsSpatialIndex( *p->mSpatialIndex ) : nullptr )  // just shallow copy
   , mSubsetString( p->mSubsetString )
+  , mCrs( p->mCrs )
 {
   mExpressionContext << QgsExpressionContextUtils::globalScope()
                      << QgsExpressionContextUtils::projectScope( QgsProject::instance() );

src/core/providers/memory/qgsmemoryfeatureiterator.h

@@ -42,6 +42,7 @@ class QgsMemoryFeatureSource : public QgsAbstractFeatureSource
     std::unique_ptr< QgsSpatialIndex > mSpatialIndex;
     QString mSubsetString;
     QgsExpressionContext mExpressionContext;
+    QgsCoordinateReferenceSystem mCrs;
 
     friend class QgsMemoryFeatureIterator;
 };
@@ -67,11 +68,13 @@ class QgsMemoryFeatureIterator : public QgsAbstractFeatureIteratorFromSource<Qgs
 
     QgsGeometry mSelectRectGeom;
     std::unique_ptr< QgsGeometryEngine > mSelectRectEngine;
+    QgsRectangle mFilterRect;
     QgsFeatureMap::const_iterator mSelectIterator;
     bool mUsingFeatureIdList = false;
     QList<QgsFeatureId> mFeatureIdList;
     QList<QgsFeatureId>::const_iterator mFeatureIdListIterator;
     QgsExpression *mSubsetExpression = nullptr;
+    QgsCoordinateTransform mTransform;
 
 };
 

src/core/qgsfeatureiterator.cpp

@@ -16,7 +16,7 @@
 #include "qgslogger.h"
 
 #include "qgssimplifymethod.h"
-
+#include "qgscsexception.h"
 #include "qgsexpressionsorter.h"
 
 QgsAbstractFeatureIterator::QgsAbstractFeatureIterator( const QgsFeatureRequest &request )
@@ -98,6 +98,43 @@ bool QgsAbstractFeatureIterator::nextFeatureFilterFids( QgsFeature &f )
   return false;
 }
 
+void QgsAbstractFeatureIterator::transformFeatureGeometry( QgsFeature &feature, const QgsCoordinateTransform &transform ) const
+{
+  if ( transform.isValid() && feature.hasGeometry() )
+  {
+    try
+    {
+      QgsGeometry g = feature.geometry();
+      g.transform( transform );
+      feature.setGeometry( g );
+    }
+    catch ( QgsCsException & )
+    {
+      // transform error
+      if ( mRequest.transformErrorCallback() )
+      {
+        mRequest.transformErrorCallback()( feature );
+      }
+    }
+  }
+}
+
+QgsRectangle QgsAbstractFeatureIterator::transformedFilterRect( const QgsCoordinateTransform &transform ) const
+{
+  if ( mRequest.filterRect().isNull() )
+    return QgsRectangle();
+
+  try
+  {
+    return transform.transformBoundingBox( mRequest.filterRect(), QgsCoordinateTransform::ReverseTransform );
+  }
+  catch ( QgsCsException & )
+  {
+    // can't reproject mFilterRect
+    return mRequest.filterRect();
+  }
+}
+
 void QgsAbstractFeatureIterator::ref()
 {
   // Prepare if required the simplification of geometries to fetch:

src/core/qgsfeatureiterator.h

@@ -114,6 +114,27 @@ class CORE_EXPORT QgsAbstractFeatureIterator
      */
     virtual bool nextFeatureFilterFids( QgsFeature &f );
 
+    /**
+     * Transforms \a feature's geometry according to the specified coordinate \a transform.
+     * If \a feature has no geometry or \a transform is invalid then calling this method
+     * has no effect and will be shortcut.
+     * Iterators should call this method before returning features to ensure that any
+     * QgsFeatureRequest::destinationCrs() set on the request is respected.
+     * \since QGIS 3.0
+     */
+    void transformFeatureGeometry( QgsFeature &feature, const QgsCoordinateTransform &transform ) const;
+
+
+    /**
+     * Returns a rectangle representing the original request's QgsFeatureRequest::filterRect().
+     * If \a transform is a valid coordinate transform, the return rectangle will represent
+     * the requested filterRect() transformed to the source's coordinate reference system.
+     * Iterators should call this method and use the returned rectangle for filtering
+     * features to ensure that any QgsFeatureRequest::destinationCrs() set on the request is respected.
+     * \since QGIS 3.0
+     */
+    QgsRectangle transformedFilterRect( const QgsCoordinateTransform &transform ) const;
+
     //! A copy of the feature request.
     QgsFeatureRequest mRequest;
 

src/core/qgsfeaturerequest.cpp

@@ -83,6 +83,8 @@ QgsFeatureRequest &QgsFeatureRequest::operator=( const QgsFeatureRequest &rh )
   mSimplifyMethod = rh.mSimplifyMethod;
   mLimit = rh.mLimit;
   mOrderBy = rh.mOrderBy;
+  mCrs = rh.mCrs;
+  mTransformErrorCallback = rh.mTransformErrorCallback;
   return *this;
 }
 
@@ -234,6 +236,24 @@ QgsFeatureRequest &QgsFeatureRequest::setSimplifyMethod( const QgsSimplifyMethod
   return *this;
 }
 
+
+QgsCoordinateReferenceSystem QgsFeatureRequest::destinationCrs() const
+{
+  return mCrs;
+}
+
+QgsFeatureRequest &QgsFeatureRequest::setDestinationCrs( const QgsCoordinateReferenceSystem &crs )
+{
+  mCrs = crs;
+  return *this;
+}
+
+QgsFeatureRequest &QgsFeatureRequest::setTransformErrorCallback( std::function<void ( const QgsFeature & )> callback )
+{
+  mTransformErrorCallback = callback;
+  return *this;
+}
+
 bool QgsFeatureRequest::acceptFeature( const QgsFeature &feature )
 {
   if ( !mFilterRect.isNull() )