QGIS/python/plugins/processing/algs/qgis/ConcaveHull.py

# -*- coding: utf-8 -*-

"""
***************************************************************************
    ConcaveHull.py
    ---------------------
    Date                 : May 2014
    Copyright            : (C) 2012 by Piotr Pociask
    Email                : piotr dot pociask at gis-support dot pl
***************************************************************************
*                                                                         *
*   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.                                   *
*                                                                         *
***************************************************************************
"""
from builtins import range

__author__ = 'Piotr Pociask'
__date__ = 'May 2014'
__copyright__ = '(C) 2014, Piotr Pociask'

# This will get replaced with a git SHA1 when you do a git archive

__revision__ = '$Format:%H$'

from math import sqrt

from qgis.core import (QgsFeature,
                       QgsFeatureSink,
                       QgsWkbTypes,
                       QgsProcessing,
                       QgsProcessingException,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterNumber,
                       QgsProcessingParameterBoolean,
                       QgsProcessingParameterFeatureSink)
import processing
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm


class ConcaveHull(QgisAlgorithm):

    INPUT = 'INPUT'
    ALPHA = 'ALPHA'
    HOLES = 'HOLES'
    NO_MULTIGEOMETRY = 'NO_MULTIGEOMETRY'
    OUTPUT = 'OUTPUT'

    def group(self):
        return self.tr('Vector geometry')

    def __init__(self):
        super().__init__()

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT, self.tr('Input point layer'), [QgsProcessing.TypeVectorPoint]))
        self.addParameter(QgsProcessingParameterNumber(self.ALPHA,
                                                       self.tr('Threshold (0-1, where 1 is equivalent with Convex Hull)'),
                                                       minValue=0, maxValue=1, defaultValue=0.3, type=QgsProcessingParameterNumber.Double))

        self.addParameter(QgsProcessingParameterBoolean(self.HOLES,
                                                        self.tr('Allow holes'), defaultValue=True))
        self.addParameter(QgsProcessingParameterBoolean(self.NO_MULTIGEOMETRY,
                                                        self.tr('Split multipart geometry into singleparts geometries'), defaultValue=False))

        self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr('Concave hull'), type=QgsProcessing.TypeVectorPolygon))

    def name(self):
        return 'concavehull'

    def displayName(self):
        return self.tr('Concave hull')

    def processAlgorithm(self, parameters, context, feedback):
        layer = self.parameterAsSource(parameters, ConcaveHull.INPUT, context)
        alpha = self.parameterAsDouble(parameters, self.ALPHA, context)
        holes = self.parameterAsBool(parameters, self.HOLES, context)
        no_multigeom = self.parameterAsBool(parameters, self.NO_MULTIGEOMETRY, context)

        # Delaunay triangulation from input point layer
        feedback.setProgressText(self.tr('Creating Delaunay triangles...'))
        delaunay_layer = processing.run("qgis:delaunaytriangulation", {'INPUT': parameters[ConcaveHull.INPUT], 'OUTPUT': 'memory:'}, feedback=feedback, context=context)['OUTPUT']

        # Get max edge length from Delaunay triangles
        feedback.setProgressText(self.tr('Computing edges max length...'))

        features = delaunay_layer.getFeatures()
        count = delaunay_layer.featureCount()
        if count == 0:
            raise QgsProcessingException(self.tr('No Delaunay triangles created.'))

        counter = 50. / count
        lengths = []
        edges = {}
        for feat in features:
            if feedback.isCanceled():
                break

            line = feat.geometry().asPolygon()[0]
            for i in range(len(line) - 1):
                lengths.append(sqrt(line[i].sqrDist(line[i + 1])))
            edges[feat.id()] = max(lengths[-3:])
            feedback.setProgress(feat.id() * counter)
        max_length = max(lengths)

        # Get features with longest edge longer than alpha*max_length
        feedback.setProgressText(self.tr('Removing features...'))
        counter = 50. / len(edges)
        i = 0
        ids = []
        for id, max_len in list(edges.items()):
            if feedback.isCanceled():
                break

            if max_len > alpha * max_length:
                ids.append(id)
            feedback.setProgress(50 + i * counter)
            i += 1

        # Remove features
        delaunay_layer.dataProvider().deleteFeatures(ids)

        # Dissolve all Delaunay triangles
        feedback.setProgressText(self.tr('Dissolving Delaunay triangles...'))
        dissolved_layer = processing.run("native:dissolve", {'INPUT': delaunay_layer, 'OUTPUT': 'memory:'}, feedback=feedback, context=context)['OUTPUT']

        # Save result
        feedback.setProgressText(self.tr('Saving data...'))
        feat = QgsFeature()
        dissolved_layer.getFeatures().nextFeature(feat)

        # Not needed anymore, free up some resources
        del delaunay_layer
        del dissolved_layer

        (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
                                               layer.fields(), QgsWkbTypes.Polygon, layer.sourceCrs())

        geom = feat.geometry()
        if no_multigeom and geom.isMultipart():
            # Only singlepart geometries are allowed
            geom_list = geom.asGeometryCollection()
            for single_geom in geom_list:
                if feedback.isCanceled():
                    break

                single_feature = QgsFeature()
                if not holes:
                    # Delete holes
                    single_geom = single_geom.removeInteriorRings()
                single_feature.setGeometry(single_geom)
                sink.addFeature(single_feature, QgsFeatureSink.FastInsert)
        else:
            # Multipart geometries are allowed
            if not holes:
                # Delete holes
                geom = geom.removeInteriorRings()
                feat.setGeometry(geom)
            sink.addFeature(feat, QgsFeatureSink.FastInsert)

        return {self.OUTPUT: dest_id}