Merge pull request #4732 from nyalldawson/alg

[FEATURE] Subdivide algorithm for geometries

Author: nyalldawson

python/core/geometry/qgsgeometry.sip

@@ -782,6 +782,23 @@ Returns the smallest convex polygon that contains all the points in the geometry
  If ``edgesOnly`` is true than line string boundary geometries will be returned
  instead of polygons.
  An empty geometry will be returned if the diagram could not be calculated.
+.. versionadded:: 3.0
+ :rtype: QgsGeometry
+%End
+
+    QgsGeometry subdivide( int maxNodes = 256 ) const;
+%Docstring
+ Subdivides the geometry. The returned geometry will be a collection containing subdivided parts
+ from the original geometry, where no part has more then the specified maximum number of nodes (``maxNodes``).
+
+ This is useful for dividing a complex geometry into less complex parts, which are better able to be spatially
+ indexed and faster to perform further operations such as intersects on. The returned geometry parts may
+ not be valid and may contain self-intersections.
+
+ The minimum allowed value for ``maxNodes`` is 8.
+
+ Curved geometries will be segmentized before subdivision.
+
 .. versionadded:: 3.0
  :rtype: QgsGeometry
 %End
@@ -828,6 +845,16 @@ Returns a geometry representing the points shared by this geometry and other.
  :rtype: QgsGeometry
 %End
 
+    QgsGeometry clipped( const QgsRectangle &rectangle );
+%Docstring
+ Clips the geometry using the specified ``rectangle``.
+
+ Performs a fast, non-robust intersection between the geometry and
+ a ``rectangle``. The returned geometry may be invalid.
+.. versionadded:: 3.0
+ :rtype: QgsGeometry
+%End
+
     QgsGeometry combine( const QgsGeometry &geometry ) const;
 %Docstring
  Returns a geometry representing all the points in this geometry and other (a
@@ -1129,7 +1156,6 @@ Ring 0 is outer ring and can't be deleted.
 .. versionadded:: 2.10
 %End
 
-
     void draw( QPainter &p ) const;
 %Docstring
  Draws the geometry onto a QPainter

python/plugins/processing/algs/qgis/MultipartToSingleparts.py

@@ -29,7 +29,11 @@
 
 from qgis.PyQt.QtGui import QIcon
 
-from qgis.core import QgsWkbTypes, QgsProcessingUtils
+from qgis.core import (QgsWkbTypes,
+                       QgsProcessingUtils,
+                       QgsProcessingParameterFeatureSource,
+                       QgsProcessingParameterFeatureSink,
+                       QgsProcessingOutputVectorLayer)
 
 from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
 from processing.core.parameters import ParameterVector
@@ -51,8 +55,11 @@ def group(self):
 
     def __init__(self):
         super().__init__()
-        self.addParameter(ParameterVector(self.INPUT, self.tr('Input layer')))
-        self.addOutput(OutputVector(self.OUTPUT, self.tr('Single parts')))
+        self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT,
+                                                              self.tr('Input layer')))
+
+        self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr('Single parts')))
+        self.addOutput(QgsProcessingOutputVectorLayer(self.OUTPUT, self.tr('Single parts')))
 
     def name(self):
         return 'multiparttosingleparts'
@@ -61,28 +68,32 @@ def displayName(self):
         return self.tr('Multipart to singleparts')
 
     def processAlgorithm(self, parameters, context, feedback):
-        layer = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.INPUT), context)
-        geomType = QgsWkbTypes.singleType(layer.wkbType())
+        source = self.parameterAsSource(parameters, self.INPUT, context)
+        geomType = QgsWkbTypes.singleType(source.wkbType())
+
+        (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
+                                               source.fields(), geomType, source.sourceCrs())
 
-        writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(layer.fields(), geomType, layer.crs(),
-                                                                     context)
+        features = source.getFeatures()
+        total = 100.0 / source.featureCount()
 
-        features = QgsProcessingUtils.getFeatures(layer, context)
-        total = 100.0 / QgsProcessingUtils.featureCount(layer, context)
         for current, f in enumerate(features):
+            if feedback.isCanceled():
+                break
+
             input_geometry = f.geometry()
             if input_geometry:
                 if input_geometry.isMultipart():
                     for g in input_geometry.asGeometryCollection():
                         output_feature = f
                         output_feature.setGeometry(g)
-                        writer.addFeature(output_feature)
+                        sink.addFeature(output_feature)
                 else:
-                    writer.addFeature(f)
+                    sink.addFeature(f)
             else:
                 #input feature with null geometry
-                writer.addFeature(f)
+                sink.addFeature(f)
 
             feedback.setProgress(int(current * total))
 
-        del writer
+        return {self.OUTPUT: dest_id}

python/plugins/processing/algs/qgis/QGISAlgorithmProvider.py

@@ -66,7 +66,7 @@
 # from .Delaunay import Delaunay
 # from .VoronoiPolygons import VoronoiPolygons
 # from .DensifyGeometries import DensifyGeometries
-# from .MultipartToSingleparts import MultipartToSingleparts
+from .MultipartToSingleparts import MultipartToSingleparts
 # from .SimplifyGeometries import SimplifyGeometries
 # from .LinesToPolygons import LinesToPolygons
 # from .PolygonsToLines import PolygonsToLines
@@ -206,7 +206,7 @@ def getAlgs(self):
         #         ReprojectLayer(), ExportGeometryInfo(), Centroids(),
         #         Delaunay(), VoronoiPolygons(), SimplifyGeometries(),
         #         DensifyGeometries(), DensifyGeometriesInterval(),
-        #         MultipartToSingleparts(), SinglePartsToMultiparts(),
+        #         , SinglePartsToMultiparts(),
         #         PolygonsToLines(), LinesToPolygons(), ExtractNodes(),
         #         ConvexHull(), FixedDistanceBuffer(),
         #         VariableDistanceBuffer(), Dissolve(), Difference(),
@@ -271,7 +271,8 @@ def getAlgs(self):
                 ExtentFromLayer(),
                 ExtractByExpression(),
                 GridPolygon(),
-                Merge()
+                Merge(),
+                MultipartToSingleparts()
                 ]
 
         if hasPlotly:

python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml

@@ -389,17 +389,17 @@ tests:
 #        name: expected/merge_lines.gml
 #        type: vector
 #
-#  - algorithm: qgis:multiparttosingleparts
-#    name: Multiparts to singleparts
-#    params:
-#      INPUT:
-#        name: multilines.gml
-#        type: vector
-#    results:
-#      OUTPUT:
-#        name: expected/multi_to_single.gml
-#        type: vector
-#
+  - algorithm: qgis:multiparttosingleparts
+    name: Multiparts to singleparts
+    params:
+      INPUT:
+        name: multilines.gml
+        type: vector
+    results:
+      OUTPUT:
+        name: expected/multi_to_single.gml
+        type: vector
+
   - algorithm: qgis:boundingboxes
     name: Bounding boxes for lines
     params:

src/core/geometry/qgsgeometry.cpp

@@ -1655,6 +1655,30 @@ QgsGeometry QgsGeometry::delaunayTriangulation( double tolerance, bool edgesOnly
   return geos.delaunayTriangulation( tolerance, edgesOnly );
 }
 
+QgsGeometry QgsGeometry::subdivide( int maxNodes ) const
+{
+  if ( !d->geometry )
+  {
+    return QgsGeometry();
+  }
+
+  const QgsAbstractGeometry *geom = d->geometry;
+  std::unique_ptr< QgsAbstractGeometry > segmentizedCopy;
+  if ( QgsWkbTypes::isCurvedType( d->geometry->wkbType() ) )
+  {
+    segmentizedCopy.reset( d->geometry->segmentize() );
+    geom = segmentizedCopy.get();
+  }
+
+  QgsGeos geos( geom );
+  QgsAbstractGeometry *result = geos.subdivide( maxNodes );
+  if ( !result )
+  {
+    return QgsGeometry();
+  }
+  return QgsGeometry( result );
+}
+
 QgsGeometry QgsGeometry::interpolate( double distance ) const
 {
   if ( !d->geometry )
@@ -2089,17 +2113,17 @@ void QgsGeometry::mapToPixel( const QgsMapToPixel &mtp )
   }
 }
 
-#if 0
-void QgsGeometry::clip( const QgsRectangle &rect )
+QgsGeometry QgsGeometry::clipped( const QgsRectangle &rectangle )
 {
-  if ( d->geometry )
+  if ( !d->geometry || rectangle.isNull() || rectangle.isEmpty() )
   {
-    detach();
-    d->geometry->clip( rect );
-    removeWkbGeos();
+    return QgsGeometry();
   }
+
+  QgsGeos geos( d->geometry );
+  QgsAbstractGeometry *resultGeom = geos.clip( rectangle );
+  return QgsGeometry( resultGeom );
 }
-#endif
 
 void QgsGeometry::draw( QPainter &p ) const
 {

src/core/geometry/qgsgeometry.h

@@ -733,6 +733,22 @@ class CORE_EXPORT QgsGeometry
      */
     QgsGeometry delaunayTriangulation( double tolerance = 0.0, bool edgesOnly = false ) const;
 
+    /**
+     * Subdivides the geometry. The returned geometry will be a collection containing subdivided parts
+     * from the original geometry, where no part has more then the specified maximum number of nodes (\a maxNodes).
+     *
+     * This is useful for dividing a complex geometry into less complex parts, which are better able to be spatially
+     * indexed and faster to perform further operations such as intersects on. The returned geometry parts may
+     * not be valid and may contain self-intersections.
+     *
+     * The minimum allowed value for \a maxNodes is 8.
+     *
+     * Curved geometries will be segmentized before subdivision.
+     *
+     * \since QGIS 3.0
+     */
+    QgsGeometry subdivide( int maxNodes = 256 ) const;
+
     /**
      * Return interpolated point on line at distance
      * \since QGIS 1.9
@@ -765,6 +781,15 @@ class CORE_EXPORT QgsGeometry
     //! Returns a geometry representing the points shared by this geometry and other.
     QgsGeometry intersection( const QgsGeometry &geometry ) const;
 
+    /**
+     * Clips the geometry using the specified \a rectangle.
+     *
+     * Performs a fast, non-robust intersection between the geometry and
+     * a \a rectangle. The returned geometry may be invalid.
+     * \since QGIS 3.0
+     */
+    QgsGeometry clipped( const QgsRectangle &rectangle );
+
     /** Returns a geometry representing all the points in this geometry and other (a
      * union geometry operation).
      * \note this operation is not called union since its a reserved word in C++.
@@ -999,13 +1024,6 @@ class CORE_EXPORT QgsGeometry
      */
     void mapToPixel( const QgsMapToPixel &mtp );
 
-    // not implemented for 2.10
-    /* Clips the geometry using the specified rectangle
-     * \param rect clip rectangle
-     * \since QGIS 2.10
-     */
-    // void clip( const QgsRectangle& rect );
-
     /** Draws the geometry onto a QPainter
      * \param p destination QPainter
      * \since QGIS 2.10

src/core/geometry/qgsgeometryfactory.cpp

@@ -264,3 +264,39 @@ QgsAbstractGeometry *QgsGeometryFactory::geomFromWkbType( QgsWkbTypes::Type t )
       return nullptr;
   }
 }
+
+std::unique_ptr<QgsGeometryCollection> QgsGeometryFactory::createCollectionOfType( QgsWkbTypes::Type t )
+{
+  QgsWkbTypes::Type type = QgsWkbTypes::flatType( QgsWkbTypes::multiType( t ) );
+  std::unique_ptr< QgsGeometryCollection > collect;
+  switch ( type )
+  {
+    case QgsWkbTypes::MultiPoint:
+      collect.reset( new QgsMultiPointV2() );
+      break;
+    case QgsWkbTypes::MultiLineString:
+      collect.reset( new QgsMultiLineString() );
+      break;
+    case QgsWkbTypes::MultiCurve:
+      collect.reset( new QgsMultiCurve() );
+      break;
+    case QgsWkbTypes::MultiPolygon:
+      collect.reset( new QgsMultiPolygonV2() );
+      break;
+    case QgsWkbTypes::MultiSurface:
+      collect.reset( new QgsMultiSurface() );
+      break;
+    case QgsWkbTypes::GeometryCollection:
+      collect.reset( new QgsGeometryCollection() );
+      break;
+    default:
+      // should not be possible
+      return nullptr;
+  }
+  if ( QgsWkbTypes::hasM( t ) )
+    collect->addMValue();
+  if ( QgsWkbTypes::hasZ( t ) )
+    collect->addZValue();
+
+  return collect;
+}

src/core/geometry/qgsgeometryfactory.h

@@ -20,11 +20,13 @@
 
 #include "qgis_core.h"
 #include <QString>
+#include <memory>
 
 class QgsAbstractGeometry;
 class QgsLineString;
 class QgsConstWkbPtr;
 class QgsRectangle;
+class QgsGeometryCollection;
 
 //compatibility with old classes
 #include "qgspointxy.h"
@@ -70,6 +72,12 @@ class CORE_EXPORT QgsGeometryFactory
     //! Return empty geometry from wkb type
     static QgsAbstractGeometry *geomFromWkbType( QgsWkbTypes::Type t );
 
+    /**
+     * Returns a new geometry collection matching a specified WKB \a type. For instance, if
+     * type is PolygonM the returned geometry will be a QgsMultiPolygonV2 with M values.
+     */
+    static std::unique_ptr< QgsGeometryCollection > createCollectionOfType( QgsWkbTypes::Type type );
+
   private:
     static QgsLineString *linestringFromPolyline( const QgsPolyline &polyline );
 };