Use prepared geometries to check for intersection between

pixel rectangles and polygons when doing exact statistics calculations in zonal stats Sppeds up zonal stats calculation by ~10x when pixel size is large (compared to polygon sizes)
qgis · Jun 8, 2018 · ff9e7ba · ff9e7ba
1 parent 43d8567
commit ff9e7ba
Showing 1 changed file with 11 additions and 3 deletions.
diff --git a/src/analysis/vector/qgszonalstatistics.cpp b/src/analysis/vector/qgszonalstatistics.cpp
@@ -444,6 +444,13 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &p
   QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
   QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );
 
+  std::unique_ptr< QgsGeometryEngine > polyEngine( QgsGeometry::createGeometryEngine( poly.constGet( ) ) );
+  if ( !polyEngine )
+  {
+    return;
+  }
+  polyEngine->prepareGeometry();
+
   std::unique_ptr< QgsRasterBlock > block( mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY ) );
   for ( int i = 0; i < nCellsY; ++i )
   {
@@ -457,11 +464,13 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &p
       }
 
       pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );
-      if ( !pixelRectGeometry.isNull() )
+      // GEOS intersects tests on prepared geometry is MAGNITUDES faster than calculating the intersection itself,
+      // so we first test to see if there IS an intersection before doing the actual calculation
+      if ( !pixelRectGeometry.isNull() && polyEngine->intersects( pixelRectGeometry.constGet() ) )
       {
         //intersection
         QgsGeometry intersectGeometry = pixelRectGeometry.intersection( poly );
-        if ( !intersectGeometry.isNull() )
+        if ( !intersectGeometry.isEmpty() )
         {
           double intersectionArea = intersectGeometry.area();
           if ( intersectionArea >= 0.0 )
@@ -470,7 +479,6 @@ void QgsZonalStatistics::statisticsFromPreciseIntersection( const QgsGeometry &p
             stats.addValue( pixelValue, weight );
           }
         }
-        pixelRectGeometry = QgsGeometry();
       }
       currentX += cellSizeX;
     }