[FEATURE] Add method to calculate pole of inaccessibility for polygons

Implements a method in QgsGeometry and a processing algorithm to
calculate the pole of inaccessibility for a surface, which is the
most distant internal point from the boundary of the surface. This function
uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative
approach guaranteed to find the true pole of inaccessibility within a specified
tolerance. More precise tolerances require more iterations and will take longer
to calculate.

python/core/geometry/qgsgeometry.sip

@@ -510,15 +510,37 @@ class QgsGeometry
     /** Returns a simplified version of this geometry using a specified tolerance value */
     QgsGeometry simplify( double tolerance )  const;
 
-    /** Returns the center of mass of a geometry
+    /**
+     * Returns the center of mass of a geometry.
      * @note for line based geometries, the center point of the line is returned,
      * and for point based geometries, the point itself is returned
+     * @see pointOnSurface()
+     * @see poleOfInaccessibility()
      */
     QgsGeometry centroid() const;
 
-    /** Returns a point within a geometry */
+    /**
+     * Returns a point guaranteed to lie on the surface of a geometry. While the centroid()
+     * of a geometry may be located outside of the geometry itself (eg for concave shapes),
+     * the point on surface will always be inside the geometry.
+     * @see centroid()
+     * @see poleOfInaccessibility()
+     */
     QgsGeometry pointOnSurface() const;
 
+    /**
+     * Calculates the approximate pole of inaccessibility for a surface, which is the
+     * most distant internal point from the boundary of the surface. This function
+     * uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative
+     * approach guaranteed to find the true pole of inaccessibility within a specified
+     * tolerance. More precise tolerances require more iterations and will take longer
+     * to calculate.
+     * @see centroid()
+     * @see pointOnSurface()
+     * @note added in QGIS 3.0
+     */
+    QgsGeometry poleOfInaccessibility( double precision ) const;
+
     /** Returns the smallest convex polygon that contains all the points in the geometry. */
     QgsGeometry convexHull() const;
 

python/plugins/processing/algs/help/qgis.yaml

@@ -337,6 +337,9 @@ qgis:polarplot: >
 
   Two fields must be entered as parameters: one that define the category each feature two (to group features) and another one with the variable to plot (this has to be a numeric one)
 
+qgis:poleofinaccessibility: >
+  This algorithm calculates the pole of inaccessibility for a polygon layer, which is the most distant internal point from the boundary of the surface. This algorithm uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative approach guaranteed to find the true pole of inaccessibility within a specified tolerance (in layer units). More precise tolerances require more iterations and will take longer to calculate.
+
 qgis:polygoncentroids: >
   This algorithm creates a new point layer, with points representing the centroid of polygons of an input layer.
 

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

@@ -45,6 +45,9 @@ class PointOnSurface(GeoAlgorithm):
     INPUT_LAYER = 'INPUT_LAYER'
     OUTPUT_LAYER = 'OUTPUT_LAYER'
 
+    def getIcon(self):
+        return QIcon(os.path.join(pluginPath, 'images', 'ftools', 'centroids.png'))
+
     def defineCharacteristics(self):
         self.name, self.i18n_name = self.trAlgorithm('Point on surface')
         self.group, self.i18n_group = self.trAlgorithm('Vector geometry tools')

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

@@ -0,0 +1,91 @@
+# -*- coding: utf-8 -*-
+
+"""
+***************************************************************************
+    PoleOfInaccessibility.py
+    ------------------------
+    Date                 : November 2016
+    Copyright            : (C) 2016 by Nyall Dawson
+    Email                : nyall dot dawson at gmail dot com
+***************************************************************************
+*                                                                         *
+*   This program is free software; you can redistribute it and/or modify  *
+*   it under the terms of the GNU General Public License as published by  *
+*   the Free Software Foundation; either version 2 of the License, or     *
+*   (at your option) any later version.                                   *
+*                                                                         *
+***************************************************************************
+"""
+
+__author__ = 'Nyall Dawson'
+__date__ = 'November 2016'
+__copyright__ = '(C) 2016, Nyall Dawson'
+
+# This will get replaced with a git SHA1 when you do a git archive323
+
+__revision__ = '$Format:%H$'
+
+import os
+
+from qgis.core import QgsGeometry, QgsWkbTypes
+
+from qgis.PyQt.QtGui import QIcon
+
+from processing.core.GeoAlgorithm import GeoAlgorithm
+from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException
+from processing.core.parameters import ParameterVector, ParameterNumber
+from processing.core.outputs import OutputVector
+from processing.tools import dataobjects, vector
+
+pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]
+
+
+class PoleOfInaccessibility(GeoAlgorithm):
+
+    INPUT_LAYER = 'INPUT_LAYER'
+    TOLERANCE = 'TOLERANCE'
+    OUTPUT_LAYER = 'OUTPUT_LAYER'
+
+    def getIcon(self):
+        return QIcon(os.path.join(pluginPath, 'images', 'ftools', 'centroids.png'))
+
+    def defineCharacteristics(self):
+        self.name, self.i18n_name = self.trAlgorithm('Pole of Inaccessibility')
+        self.group, self.i18n_group = self.trAlgorithm('Vector geometry tools')
+
+        self.addParameter(ParameterVector(self.INPUT_LAYER,
+                                          self.tr('Input layer'),
+                                          [dataobjects.TYPE_VECTOR_POLYGON]))
+        self.addParameter(ParameterNumber(self.TOLERANCE,
+                                          self.tr('Tolerance (layer units)'), default=1.0, minValue=0.0))
+        self.addOutput(OutputVector(self.OUTPUT_LAYER, self.tr('Point'), datatype=[dataobjects.TYPE_VECTOR_POINT]))
+
+    def processAlgorithm(self, progress):
+        layer = dataobjects.getObjectFromUri(
+            self.getParameterValue(self.INPUT_LAYER))
+        tolerance = self.getParameterValue(self.TOLERANCE)
+
+        writer = self.getOutputFromName(
+            self.OUTPUT_LAYER).getVectorWriter(
+                layer.fields(),
+                QgsWkbTypes.Point,
+                layer.crs())
+
+        features = vector.features(layer)
+        total = 100.0 / len(features)
+
+        for current, input_feature in enumerate(features):
+            output_feature = input_feature
+            input_geometry = input_feature.geometry()
+            if input_geometry:
+                output_geometry = input_geometry.poleOfInaccessibility(tolerance)
+                if not output_geometry:
+                    raise GeoAlgorithmExecutionException(
+                        self.tr('Error calculating pole of inaccessibility'))
+
+                output_feature.setGeometry(output_geometry)
+
+            writer.addFeature(output_feature)
+            progress.setPercentage(int(current * total))
+
+        del writer

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

@@ -176,6 +176,7 @@
 from .ExtractSpecificNodes import ExtractSpecificNodes
 from .GeometryByExpression import GeometryByExpression
 from .SnapGeometries import SnapGeometriesToLayer
+from .PoleOfInaccessibility import PoleOfInaccessibility
 
 pluginPath = os.path.normpath(os.path.join(
     os.path.split(os.path.dirname(__file__))[0], os.pardir))
@@ -238,7 +239,8 @@ def __init__(self):
                         IdwInterpolationZValue(), IdwInterpolationAttribute(),
                         TinInterpolationZValue(), TinInterpolationAttribute(),
                         RemoveNullGeometry(), ExtractByExpression(), ExtendLines(),
-                        ExtractSpecificNodes(), GeometryByExpression(), SnapGeometriesToLayer()
+                        ExtractSpecificNodes(), GeometryByExpression(), SnapGeometriesToLayer(),
+                        PoleOfInaccessibility()
                         ]
 
         if hasMatplotlib:

python/plugins/processing/tests/testdata/expected/pole_inaccessibility_polys.gfs

@@ -0,0 +1,32 @@
+<GMLFeatureClassList>
+  <GMLFeatureClass>
+    <Name>pole_inaccessibility_polys</Name>
+    <ElementPath>pole_inaccessibility_polys</ElementPath>
+    <!--POINT-->
+    <GeometryType>1</GeometryType>
+    <SRSName>EPSG:4326</SRSName>
+    <DatasetSpecificInfo>
+      <FeatureCount>6</FeatureCount>
+      <ExtentXMin>0.50000</ExtentXMin>
+      <ExtentXMax>6.58578</ExtentXMax>
+      <ExtentYMin>-2.41422</ExtentYMin>
+      <ExtentYMax>5.50000</ExtentYMax>
+    </DatasetSpecificInfo>
+    <PropertyDefn>
+      <Name>name</Name>
+      <ElementPath>name</ElementPath>
+      <Type>String</Type>
+      <Width>5</Width>
+    </PropertyDefn>
+    <PropertyDefn>
+      <Name>intval</Name>
+      <ElementPath>intval</ElementPath>
+      <Type>Integer</Type>
+    </PropertyDefn>
+    <PropertyDefn>
+      <Name>floatval</Name>
+      <ElementPath>floatval</ElementPath>
+      <Type>Real</Type>
+    </PropertyDefn>
+  </GMLFeatureClass>
+</GMLFeatureClassList>

python/plugins/processing/tests/testdata/expected/pole_inaccessibility_polys.gml

@@ -0,0 +1,58 @@
+<?xml version="1.0" encoding="utf-8" ?>
+<ogr:FeatureCollection
+     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+     xsi:schemaLocation=""
+     xmlns:ogr="http://ogr.maptools.org/"
+     xmlns:gml="http://www.opengis.net/gml">
+  <gml:boundedBy>
+    <gml:Box>
+      <gml:coord><gml:X>0.5</gml:X><gml:Y>-2.414215087890625</gml:Y></gml:coord>
+      <gml:coord><gml:X>6.585784912109375</gml:X><gml:Y>5.5</gml:Y></gml:coord>
+    </gml:Box>
+  </gml:boundedBy>
+
+  <gml:featureMember>
+    <ogr:pole_inaccessibility_polys fid="polys.0">
+      <ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>0.5,0.5</gml:coordinates></gml:Point></ogr:geometryProperty>
+      <ogr:name>aaaaa</ogr:name>
+      <ogr:intval>33</ogr:intval>
+      <ogr:floatval>44.123456</ogr:floatval>
+    </ogr:pole_inaccessibility_polys>
+  </gml:featureMember>
+  <gml:featureMember>
+    <ogr:pole_inaccessibility_polys fid="polys.1">
+      <ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>4.999996185302734,4.414211273193359</gml:coordinates></gml:Point></ogr:geometryProperty>
+      <ogr:name>Aaaaa</ogr:name>
+      <ogr:intval>-33</ogr:intval>
+      <ogr:floatval>0</ogr:floatval>
+    </ogr:pole_inaccessibility_polys>
+  </gml:featureMember>
+  <gml:featureMember>
+    <ogr:pole_inaccessibility_polys fid="polys.2">
+      <ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>2.5,5.5</gml:coordinates></gml:Point></ogr:geometryProperty>
+      <ogr:name>bbaaa</ogr:name>
+      <ogr:floatval>0.123</ogr:floatval>
+    </ogr:pole_inaccessibility_polys>
+  </gml:featureMember>
+  <gml:featureMember>
+    <ogr:pole_inaccessibility_polys fid="polys.3">
+      <ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>6.585784912109375,-2.414215087890625</gml:coordinates></gml:Point></ogr:geometryProperty>
+      <ogr:name>ASDF</ogr:name>
+      <ogr:intval>0</ogr:intval>
+    </ogr:pole_inaccessibility_polys>
+  </gml:featureMember>
+  <gml:featureMember>
+    <ogr:pole_inaccessibility_polys fid="polys.4">
+      <ogr:intval>120</ogr:intval>
+      <ogr:floatval>-100291.43213</ogr:floatval>
+    </ogr:pole_inaccessibility_polys>
+  </gml:featureMember>
+  <gml:featureMember>
+    <ogr:pole_inaccessibility_polys fid="polys.5">
+      <ogr:geometryProperty><gml:Point srsName="EPSG:4326"><gml:coordinates>4.289703369140625,-0.232696533203125</gml:coordinates></gml:Point></ogr:geometryProperty>
+      <ogr:name>elim</ogr:name>
+      <ogr:intval>2</ogr:intval>
+      <ogr:floatval>3.33</ogr:floatval>
+    </ogr:pole_inaccessibility_polys>
+  </gml:featureMember>
+</ogr:FeatureCollection>

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

@@ -1464,3 +1464,15 @@ tests:
       OUTPUT:
         name: expected/snap_point_to_lines_prefer_closest.gml
         type: vector
+
+  - algorithm: qgis:poleofinaccessibility
+    name: Pole of inaccessibility (polygons)
+    params:
+      INPUT_LAYER:
+        name: polys.gml
+        type: vector
+      TOLERANCE: 1.0e-05
+    results:
+      OUTPUT_LAYER:
+        name: expected/pole_inaccessibility_polys.gml
+        type: vector

src/core/geometry/qgsgeometry.cpp

@@ -1480,6 +1480,13 @@ QgsGeometry QgsGeometry::pointOnSurface() const
   return QgsGeometry( pt.clone() );
 }
 
+QgsGeometry QgsGeometry::poleOfInaccessibility( double precision ) const
+{
+  QgsInternalGeometryEngine engine( *this );
+
+  return engine.poleOfInaccessibility( precision );
+}
+
 QgsGeometry QgsGeometry::convexHull() const
 {
   if ( !d->geometry )

src/core/geometry/qgsgeometry.h

@@ -563,15 +563,37 @@ class CORE_EXPORT QgsGeometry
     //! Returns a simplified version of this geometry using a specified tolerance value
     QgsGeometry simplify( double tolerance ) const;
 
-    /** Returns the center of mass of a geometry
+    /**
+     * Returns the center of mass of a geometry.
      * @note for line based geometries, the center point of the line is returned,
      * and for point based geometries, the point itself is returned
+     * @see pointOnSurface()
+     * @see poleOfInaccessibility()
      */
     QgsGeometry centroid() const;
 
-    //! Returns a point within a geometry
+    /**
+     * Returns a point guaranteed to lie on the surface of a geometry. While the centroid()
+     * of a geometry may be located outside of the geometry itself (eg for concave shapes),
+     * the point on surface will always be inside the geometry.
+     * @see centroid()
+     * @see poleOfInaccessibility()
+     */
     QgsGeometry pointOnSurface() const;
 
+    /**
+     * Calculates the approximate pole of inaccessibility for a surface, which is the
+     * most distant internal point from the boundary of the surface. This function
+     * uses the 'polylabel' algorithm (Vladimir Agafonkin, 2016), which is an iterative
+     * approach guaranteed to find the true pole of inaccessibility within a specified
+     * tolerance. More precise tolerances require more iterations and will take longer
+     * to calculate.
+     * @see centroid()
+     * @see pointOnSurface()
+     * @note added in QGIS 3.0
+     */
+    QgsGeometry poleOfInaccessibility( double precision ) const;
+
     //! Returns the smallest convex polygon that contains all the points in the geometry.
     QgsGeometry convexHull() const;