|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | + |
| 3 | +""" |
| 4 | +*************************************************************************** |
| 5 | + ConcaveHull.py |
| 6 | + --------------------- |
| 7 | + Date : May 2014 |
| 8 | + Copyright : (C) 2012 by Piotr Pociask |
| 9 | + Email : piotr dot pociask at gis-support dot pl |
| 10 | +*************************************************************************** |
| 11 | +* * |
| 12 | +* This program is free software; you can redistribute it and/or modify * |
| 13 | +* it under the terms of the GNU General Public License as published by * |
| 14 | +* the Free Software Foundation; either version 2 of the License, or * |
| 15 | +* (at your option) any later version. * |
| 16 | +* * |
| 17 | +*************************************************************************** |
| 18 | +""" |
| 19 | + |
| 20 | +__author__ = 'Piotr Pociask' |
| 21 | +__date__ = 'May 2014' |
| 22 | +__copyright__ = '(C) 2014, Piotr Pociask' |
| 23 | + |
| 24 | +# This will get replaced with a git SHA1 when you do a git archive |
| 25 | + |
| 26 | +__revision__ = '$Format:%H$' |
| 27 | + |
| 28 | +from PyQt4.QtCore import * |
| 29 | +from PyQt4.QtGui import * |
| 30 | +from qgis.core import * |
| 31 | +from processing.core.GeoAlgorithm import GeoAlgorithm |
| 32 | +from processing.parameters.ParameterVector import ParameterVector |
| 33 | +from processing.parameters.ParameterNumber import ParameterNumber |
| 34 | +from processing.parameters.ParameterBoolean import ParameterBoolean |
| 35 | +from processing.outputs.OutputVector import OutputVector |
| 36 | +from processing.tools import dataobjects |
| 37 | +import processing |
| 38 | +from math import sqrt |
| 39 | + |
| 40 | +class ConcaveHull(GeoAlgorithm): |
| 41 | + |
| 42 | + INPUT = 'INPUT' |
| 43 | + ALPHA = 'ALPHA' |
| 44 | + HOLES = 'HOLES' |
| 45 | + NO_MULTIGEOMETRY = 'NO_MULTIGEOMETRY' |
| 46 | + OUTPUT = 'OUTPUT' |
| 47 | + |
| 48 | + def defineCharacteristics(self): |
| 49 | + self.name = 'Concave hull' |
| 50 | + self.group = 'Vector geometry tools' |
| 51 | + self.addParameter(ParameterVector(ConcaveHull.INPUT, 'Input point layer', |
| 52 | + [ParameterVector.VECTOR_TYPE_POINT])) |
| 53 | + self.addParameter(ParameterNumber(self.ALPHA, |
| 54 | + 'Threshold (0-1, where 1 is equivalent with Convex Hull)', 0, 1, 0.3)) |
| 55 | + self.addParameter(ParameterBoolean(self.HOLES, 'Allow holes', True)) |
| 56 | + self.addParameter(ParameterBoolean(self.NO_MULTIGEOMETRY, |
| 57 | + 'Split multipart geometry into singleparts geometries', False)) |
| 58 | + self.addOutput(OutputVector(ConcaveHull.OUTPUT, 'Convex hull')) |
| 59 | + |
| 60 | + def processAlgorithm(self, progress): |
| 61 | + #get parameters |
| 62 | + layer = dataobjects.getObjectFromUri(self.getParameterValue(ConcaveHull.INPUT)) |
| 63 | + alpha = self.getParameterValue(self.ALPHA) |
| 64 | + holes = self.getParameterValue(self.HOLES) |
| 65 | + no_multigeom = self.getParameterValue(self.NO_MULTIGEOMETRY) |
| 66 | + #Delaunay triangulation from input point layer |
| 67 | + progress.setText('Creating Delaunay triangles ...') |
| 68 | + delone_triangles = processing.runalg("qgis:delaunaytriangulation", layer, None)['OUTPUT'] |
| 69 | + delaunay_layer = processing.getObject(delone_triangles) |
| 70 | + #get max edge length from Delaunay triangles |
| 71 | + progress.setText('Computing edges max length ...') |
| 72 | + features = delaunay_layer.getFeatures() |
| 73 | + counter = 50./delaunay_layer.featureCount() |
| 74 | + lengths = [] |
| 75 | + edges = {} |
| 76 | + for feat in features: |
| 77 | + line = feat.geometry().asPolygon()[0] |
| 78 | + for i in range(len(line)-1): |
| 79 | + lengths.append(sqrt(line[i].sqrDist(line[i+1]))) |
| 80 | + edges[feat.id()] = max(lengths[-3:]) |
| 81 | + progress.setPercentage(feat.id()*counter) |
| 82 | + max_length = max(lengths) |
| 83 | + #get features with longest edge longer than alpha*max_length |
| 84 | + progress.setText('Removing features ...') |
| 85 | + counter = 50./len(edges) |
| 86 | + i = 0 |
| 87 | + ids = [] |
| 88 | + for id, max_len in edges.iteritems(): |
| 89 | + if max_len > alpha*max_length: |
| 90 | + ids.append(id) |
| 91 | + progress.setPercentage(50+i*counter) |
| 92 | + i += 1 |
| 93 | + #remove features |
| 94 | + delaunay_layer.setSelectedFeatures(ids) |
| 95 | + delaunay_layer.startEditing() |
| 96 | + delaunay_layer.deleteSelectedFeatures() |
| 97 | + delaunay_layer.commitChanges() |
| 98 | + #dissolve all Delaunay triangles |
| 99 | + progress.setText('Dissolving Delaunay triangles ...') |
| 100 | + dissolved = processing.runalg("qgis:dissolve", delaunay_layer, |
| 101 | + True, '', None)['OUTPUT'] |
| 102 | + dissolved_layer = processing.getObject(dissolved) |
| 103 | + #save result |
| 104 | + progress.setText('Saving data ...') |
| 105 | + feat = QgsFeature() |
| 106 | + dissolved_layer.getFeatures(QgsFeatureRequest().setFilterFid(0)).nextFeature(feat) |
| 107 | + writer = self.getOutputFromName( |
| 108 | + self.OUTPUT).getVectorWriter(layer.pendingFields().toList(), |
| 109 | + QGis.WKBPolygon, layer.crs()) |
| 110 | + geom = feat.geometry() |
| 111 | + if no_multigeom and geom.isMultipart(): |
| 112 | + #only singlepart geometries are allowed |
| 113 | + geom_list = geom.asMultiPolygon() |
| 114 | + for single_geom_list in geom_list: |
| 115 | + single_feature = QgsFeature() |
| 116 | + single_geom = QgsGeometry.fromPolygon(single_geom_list) |
| 117 | + if not holes: |
| 118 | + #delete holes |
| 119 | + deleted = True |
| 120 | + while deleted: |
| 121 | + deleted = single_geom.deleteRing(1) |
| 122 | + single_feature.setGeometry(single_geom) |
| 123 | + writer.addFeature(single_feature) |
| 124 | + else: |
| 125 | + #multipart geometries are allowed |
| 126 | + if not holes: |
| 127 | + #delete holes |
| 128 | + deleted = True |
| 129 | + while deleted: |
| 130 | + deleted = geom.deleteRing(1) |
| 131 | + writer.addFeature(feat) |
| 132 | + del writer |
0 commit comments