Much faster save as vector layer when using restricted extent

Use a provider filter rect to avoid looping through every single feature in this case, and also use a prepared extent geometry to speed up the remaining intersection checks. (cherry-picked from d2f3eb1)
qgis · Apr 10, 2017 · bf51886 · bf51886
1 parent 248d8c9
commit bf51886
Show file tree

Hide file tree

Showing 2 changed files with 90 additions and 6 deletions.
diff --git a/src/core/qgsvectorfilewriter.cpp b/src/core/qgsvectorfilewriter.cpp
@@ -28,6 +28,8 @@
 #include "qgssymbollayerv2.h"
 #include "qgsvectordataprovider.h"
 #include "qgslocalec.h"
+#include "qgscsexception.h"
+#include "qgsgeometryengine.h"
 
 #include <QFile>
 #include <QSettings>
@@ -2443,6 +2445,34 @@ QgsVectorFileWriter::writeAsVectorFormat( QgsVectorLayer* layer,
   req.setSubsetOfAttributes( attributes );
   if ( options.onlySelectedFeatures )
     req.setFilterFids( layer->selectedFeaturesIds() );
+
+  QScopedPointer< QgsGeometry > filterRectGeometry;
+  QScopedPointer< QgsGeometryEngine  > filterRectEngine;
+  if ( !options.filterExtent.isNull() )
+  {
+    QgsRectangle filterRect = options.filterExtent;
+    bool useFilterRect = true;
+    if ( shallTransform )
+    {
+      try
+      {
+        // map filter rect back from destination CRS to layer CRS
+        filterRect = options.ct->transformBoundingBox( filterRect, QgsCoordinateTransform::ReverseTransform );
+      }
+      catch ( QgsCsException & )
+      {
+        useFilterRect = false;
+      }
+    }
+    if ( useFilterRect )
+    {
+      req.setFilterRect( filterRect );
+    }
+    filterRectGeometry.reset( QgsGeometry::fromRect( options.filterExtent ) );
+    filterRectEngine.reset( QgsGeometry::createGeometryEngine( filterRectGeometry->geometry() ) );
+    filterRectEngine->prepareGeometry();
+  }
+
   QgsFeatureIterator fit = layer->getFeatures( req );
 
   //create symbol table if needed
@@ -2494,7 +2524,7 @@ QgsVectorFileWriter::writeAsVectorFormat( QgsVectorLayer* layer,
     {
       try
       {
-        if ( fet.geometry() )
+        if ( fet.constGeometry() )
         {
           fet.geometry()->transform( *( options.ct ) );
         }
@@ -2513,7 +2543,7 @@ QgsVectorFileWriter::writeAsVectorFormat( QgsVectorLayer* layer,
       }
     }
 
-    if ( fet.constGeometry() && !options.filterExtent.isNull() && !fet.constGeometry()->intersects( options.filterExtent ) )
+    if ( fet.constGeometry() && filterRectEngine && !filterRectEngine->intersects( *fet.constGeometry()->geometry() ) )
       continue;
 
     if ( attributes.size() < 1 && options.skipAttributeCreation )

diff --git a/tests/src/python/test_qgsvectorfilewriter.py b/tests/src/python/test_qgsvectorfilewriter.py
@@ -24,16 +24,19 @@
                        QgsCoordinateReferenceSystem,
                        QgsVectorFileWriter,
                        QgsFeatureRequest,
-                       QgsWKBTypes
+                       QgsWKBTypes,
+                       QgsRectangle,
+                       QgsCoordinateTransform
                        )
 from qgis.PyQt.QtCore import QDate, QTime, QDateTime, QVariant, QDir
 import os
 import osgeo.gdal
 from osgeo import gdal, ogr
 import platform
 from qgis.testing import start_app, unittest
-from utilities import writeShape, compareWkt
+from utilities import writeShape, compareWkt, unitTestDataPath
 
+TEST_DATA_DIR = unitTestDataPath()
 start_app()
 
 
@@ -66,8 +69,7 @@ def convert(self, idx, value):
         return 'unexpected_idx'
 
 
-class TestQgsVectorLayer(unittest.TestCase):
-
+class TestQgsVectorFileWriter(unittest.TestCase):
     mMemoryLayer = None
 
     def testWrite(self):
@@ -145,6 +147,58 @@ def testDateTimeWriteShapefile(self):
         self.assertIsInstance(f.attributes()[datetime_idx], str)
         self.assertEqual(f.attributes()[datetime_idx], QDateTime(QDate(2014, 3, 5), QTime(13, 45, 22)).toString("yyyy/MM/dd hh:mm:ss.zzz"))
 
+    def testWriterWithExtent(self):
+        """Check writing using extent filter."""
+        source_file = os.path.join(TEST_DATA_DIR, 'points.shp')
+        source_layer = QgsVectorLayer(source_file, 'Points', 'ogr')
+        self.assertTrue(source_layer.isValid())
+
+        options = QgsVectorFileWriter.SaveVectorOptions()
+        options.driverName = 'ESRI Shapefile'
+        options.filterExtent = QgsRectangle(-111, 26, -96, 38)
+
+        dest_file_name = os.path.join(str(QDir.tempPath()), 'extent_no_transform.shp')
+        write_result = QgsVectorFileWriter.writeAsVectorFormat(
+            source_layer,
+            dest_file_name,
+            options)
+        self.assertEqual(write_result, QgsVectorFileWriter.NoError)
+
+        # Open result and check
+        created_layer = QgsVectorLayer('{}|layerid=0'.format(dest_file_name), 'test', 'ogr')
+        features = [f for f in created_layer.getFeatures()]
+        self.assertEqual(len(features), 5)
+        for f in features:
+            self.assertTrue(f.geometry().intersects(options.filterExtent))
+
+    def testWriterWithExtentAndReprojection(self):
+        """Check writing using extent filter with reprojection."""
+        source_file = os.path.join(TEST_DATA_DIR, 'points.shp')
+        source_layer = QgsVectorLayer(source_file, 'Points', 'ogr')
+        self.assertTrue(source_layer.isValid())
+
+        options = QgsVectorFileWriter.SaveVectorOptions()
+        options.driverName = 'ESRI Shapefile'
+        options.filterExtent = QgsRectangle(-12511460, 3045157, -10646621, 4683497)
+        crs = QgsCoordinateReferenceSystem()
+        self.assertTrue(crs.createFromOgcWmsCrs('EPSG:3785'))
+        ct = QgsCoordinateTransform(source_layer.crs(), crs)
+        options.ct = ct
+
+        dest_file_name = os.path.join(str(QDir.tempPath()), 'extent_transform.shp')
+        write_result = QgsVectorFileWriter.writeAsVectorFormat(
+            source_layer,
+            dest_file_name,
+            options)
+        self.assertEqual(write_result, QgsVectorFileWriter.NoError)
+
+        # Open result and check
+        created_layer = QgsVectorLayer('{}|layerid=0'.format(dest_file_name), 'test', 'ogr')
+        features = [f for f in created_layer.getFeatures()]
+        self.assertEqual(len(features), 5)
+        for f in features:
+            self.assertTrue(f.geometry().intersects(options.filterExtent))
+
     def testDateTimeWriteTabfile(self):
         """Check writing date and time fields to an MapInfo tabfile."""
         ml = QgsVectorLayer(