1# -*- coding: utf-8 -*- 2 3""" 4*************************************************************************** 5 VoronoiPolygons.py 6 --------------------- 7 Date : August 2012 8 Copyright : (C) 2012 by Victor Olaya 9 Email : volayaf at gmail dot com 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__authors__ = 'Victor Olaya, Håvard Tveite' 21__date__ = 'August 2012' 22__copyright__ = '(C) 2012, Victor Olaya' 23 24import os 25 26from qgis.PyQt.QtGui import QIcon 27 28from qgis.core import (QgsApplication, 29 QgsFeatureRequest, 30 QgsFeatureSink, 31 QgsFeature, 32 QgsGeometry, 33 QgsPointXY, 34 QgsWkbTypes, 35 QgsProcessing, 36 QgsProcessingException, 37 QgsProcessingParameterFeatureSource, 38 QgsProcessingParameterFeatureSink, 39 QgsProcessingParameterNumber) 40 41from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm 42 43from . import voronoi 44 45pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0] 46 47 48class VoronoiPolygons(QgisAlgorithm): 49 INPUT = 'INPUT' 50 BUFFER = 'BUFFER' 51 OUTPUT = 'OUTPUT' 52 53 def icon(self): 54 return QgsApplication.getThemeIcon("/algorithms/mAlgorithmVoronoi.svg") 55 56 def svgIconPath(self): 57 return QgsApplication.iconPath("/algorithms/mAlgorithmVoronoi.svg") 58 59 def group(self): 60 return self.tr('Vector geometry') 61 62 def groupId(self): 63 return 'vectorgeometry' 64 65 def __init__(self): 66 super().__init__() 67 68 def initAlgorithm(self, config=None): 69 self.addParameter( 70 QgsProcessingParameterFeatureSource( 71 self.INPUT, self.tr('Input layer'), 72 [QgsProcessing.TypeVectorPoint])) 73 self.addParameter( 74 QgsProcessingParameterNumber( 75 self.BUFFER, self.tr('Buffer region (% of extent)'), 76 minValue=0.0, defaultValue=0.0)) 77 self.addParameter( 78 QgsProcessingParameterFeatureSink( 79 self.OUTPUT, self.tr('Voronoi polygons'), 80 type=QgsProcessing.TypeVectorPolygon)) 81 82 def name(self): 83 return 'voronoipolygons' 84 85 def displayName(self): 86 return self.tr('Voronoi polygons') 87 88 def tags(self): 89 return self.tr('thiessen').split(',') 90 91 def processAlgorithm(self, parameters, context, feedback): 92 source = self.parameterAsSource(parameters, self.INPUT, context) 93 if source is None: 94 raise QgsProcessingException( 95 self.invalidSourceError(parameters, self.INPUT)) 96 97 buf = self.parameterAsDouble(parameters, self.BUFFER, context) 98 (sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, 99 context, source.fields(), 100 QgsWkbTypes.Polygon, 101 source.sourceCrs()) 102 if sink is None: 103 raise QgsProcessingException( 104 self.invalidSinkError(parameters, self.OUTPUT)) 105 outFeat = QgsFeature() 106 extent = source.sourceExtent() 107 extraX = extent.width() * (buf / 100.0) 108 # Adjust the extent 109 extent.setXMinimum(extent.xMinimum() - extraX) 110 extent.setXMaximum(extent.xMaximum() + extraX) 111 extraY = extent.height() * (buf / 100.0) 112 extent.setYMinimum(extent.yMinimum() - extraY) 113 extent.setYMaximum(extent.yMaximum() + extraY) 114 height = extent.height() 115 width = extent.width() 116 c = voronoi.Context() 117 pts = [] 118 ptDict = {} 119 ptNdx = -1 120 # Find the minimum and maximum x and y for the input points 121 xmin = width 122 xmax = 0 123 ymin = height 124 ymax = 0 125 features = source.getFeatures() 126 total = 100.0 / source.featureCount() if source.featureCount() else 0 127 for current, inFeat in enumerate(features): 128 if feedback.isCanceled(): 129 break 130 geom = inFeat.geometry() 131 point = geom.asPoint() 132 x = point.x() - extent.xMinimum() 133 y = point.y() - extent.yMinimum() 134 pts.append((x, y)) 135 ptNdx += 1 136 ptDict[ptNdx] = inFeat.id() 137 if x < xmin: 138 xmin = x 139 if y < ymin: 140 ymin = y 141 if x > xmax: 142 xmax = x 143 if y > ymax: 144 ymax = y 145 feedback.setProgress(int(current * total)) 146 if xmin == xmax or ymin == ymax: 147 raise QgsProcessingException('The extent of the input points is ' 148 'not a polygon (all the points are ' 149 'on a vertical or horizontal line) ' 150 '- cannot make a Voronoi diagram!') 151 xyminmax = [xmin, ymin, xmax, ymax] 152 if len(pts) < 3: 153 raise QgsProcessingException( 154 self.tr('Input file should contain at least 3 points. Choose ' 155 'another file and try again.')) 156 # Eliminate duplicate points 157 uniqueSet = set(item for item in pts) 158 ids = [pts.index(item) for item in uniqueSet] 159 sl = voronoi.SiteList([voronoi.Site(i[0], i[1], sitenum=j) 160 for (j, i) in enumerate(uniqueSet)]) 161 voronoi.voronoi(sl, c) 162 if len(c.polygons) == 0: 163 raise QgsProcessingException( 164 self.tr('There were no polygons created.')) 165 166 inFeat = QgsFeature() 167 current = 0 168 total = 100.0 / len(c.polygons) 169 # Clip each of the generated "polygons" 170 for (site, edges) in list(c.polygons.items()): 171 if feedback.isCanceled(): 172 break 173 request = QgsFeatureRequest().setFilterFid(ptDict[ids[site]]) 174 inFeat = next(source.getFeatures(request)) 175 boundarypoints = self.clip_voronoi(edges, c, width, 176 height, extent, 177 inFeat.geometry().asPoint(), 178 xyminmax) 179 ptgeom = QgsGeometry.fromMultiPointXY(boundarypoints) 180 geom = QgsGeometry(ptgeom.convexHull()) 181 outFeat.setGeometry(geom) 182 outFeat.setAttributes(inFeat.attributes()) 183 sink.addFeature(outFeat, QgsFeatureSink.FastInsert) 184 current += 1 185 feedback.setProgress(int(current * total)) 186 return {self.OUTPUT: dest_id} 187 188 def clip_voronoi(self, edges, c, width, height, extent, point, xyminmax): 189 """Clip voronoi function based on code written for Inkscape. 190 Copyright (C) 2010 Alvin Penner, penner@vaxxine.com 191 Clips one Thiessen polygon (convex polygon) to extent 192 """ 193 194 pt_x = point.x() - extent.xMinimum() 195 pt_y = point.y() - extent.yMinimum() 196 (xmin, ymin, xmax, ymax) = xyminmax 197 198 def isclose(a, b, rel_tol=1e-9, abs_tol=0.0): 199 return abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol) 200 201 def clip_line(x1, y1, x2, y2, w, h): 202 if x1 < 0 and x2 < 0: 203 # Completely to the left 204 return [0, 0, 0, 0] 205 if x1 > w and x2 > w: 206 # Completely to the right 207 return [0, 0, 0, 0] 208 if y1 < 0 and y2 < 0: 209 # Completely below 210 return [0, 0, 0, 0] 211 if y1 > h and y2 > h: 212 # Completely above 213 return [0, 0, 0, 0] 214 # Clip on the left envelope boundary 215 if x1 < 0: 216 # First point to the left 217 y1 = (y1 * x2 - y2 * x1) / (x2 - x1) 218 x1 = 0 219 if x2 < 0: 220 # Last point to the left 221 y2 = (y1 * x2 - y2 * x1) / (x2 - x1) 222 x2 = 0 223 # Clip on the right envelope boundary 224 if x1 > w: 225 # First point to the right 226 y1 = y1 + (w - x1) * (y2 - y1) / (x2 - x1) 227 x1 = w 228 if x2 > w: 229 # Last point to the right 230 y2 = y1 + (w - x1) * (y2 - y1) / (x2 - x1) 231 x2 = w 232 if isclose(x1, x2) and isclose(y1, y2): 233 return [0, 0, 0, 0] 234 # Clip on the bottom envelope boundary 235 if y1 < 0: 236 # First point below 237 x1 = (x1 * y2 - x2 * y1) / (y2 - y1) 238 y1 = 0 239 if y2 < 0: 240 # Second point below 241 x2 = (x1 * y2 - x2 * y1) / (y2 - y1) 242 y2 = 0 243 # Clip on the top envelope boundary 244 if y1 > h: 245 # First point above 246 x1 = x1 + (h - y1) * (x2 - x1) / (y2 - y1) 247 y1 = h 248 if y2 > h: 249 # Second point above 250 x2 = x1 + (h - y1) * (x2 - x1) / (y2 - y1) 251 y2 = h 252 if isclose(x1, x2) and isclose(y1, y2): 253 return [0, 0, 0, 0] 254 return [x1, y1, x2, y2] 255 256 bndpoints = [] 257 hasXMin = False 258 hasYMin = False 259 hasXMax = False 260 hasYMax = False 261 XMinNumber = 0 262 XMaxNumber = 0 263 YMinNumber = 0 264 YMaxNumber = 0 265 # The same line may appear twice for collinear input points, 266 # so have to remember which lines have contributed 267 XMinLine = -1 268 XMaxLine = -1 269 YMinLine = -1 270 YMaxLine = -1 271 for edge in edges: 272 if edge[1] >= 0 and edge[2] >= 0: 273 # Two vertices 274 [x1, y1, x2, y2] = clip_line( 275 c.vertices[edge[1]][0], 276 c.vertices[edge[1]][1], 277 c.vertices[edge[2]][0], 278 c.vertices[edge[2]][1], 279 width, 280 height 281 ) 282 elif edge[1] >= 0: 283 # Only one (left) vertex 284 if c.lines[edge[0]][1] == 0: 285 # Vertical line 286 # xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0] 287 xtemp = c.vertices[edge[1]][0] 288 ytemp = 0 - 1 289 # if c.vertices[edge[1]][1] > height / 2: 290 # ytemp = height 291 # else: 292 # ytemp = 0 293 else: 294 # Create an end of the line at the right edge - OK 295 xtemp = width 296 ytemp = (c.lines[edge[0]][2] - width * 297 c.lines[edge[0]][0]) / c.lines[edge[0]][1] 298 [x1, y1, x2, y2] = clip_line( 299 c.vertices[edge[1]][0], 300 c.vertices[edge[1]][1], 301 xtemp, 302 ytemp, 303 width, 304 height 305 ) 306 elif edge[2] >= 0: 307 # Only one (right) vertex 308 if c.lines[edge[0]][1] == 0: 309 # Vertical line 310 # xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0] 311 xtemp = c.vertices[edge[2]][0] 312 ytemp = height + 1 313 # if c.vertices[edge[2]][1] > height / 2: 314 # ytemp = height 315 # else: 316 # ytemp = 0.0 317 else: 318 # End the line at the left edge - OK 319 xtemp = 0.0 320 ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1] 321 [x1, y1, x2, y2] = clip_line( 322 xtemp, 323 ytemp, 324 c.vertices[edge[2]][0], 325 c.vertices[edge[2]][1], 326 width, 327 height, 328 ) 329 else: 330 # No vertex, only a line 331 if c.lines[edge[0]][1] == 0: 332 # Vertical line - should not happen 333 xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0] 334 ytemp = 0.0 335 xend = xtemp 336 yend = height 337 else: 338 # End the line at both edges - ??? 339 xtemp = 0.0 340 ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1] 341 xend = width 342 yend = (c.lines[edge[0]][2] - width * 343 c.lines[edge[0]][0]) / c.lines[edge[0]][1] 344 [x1, y1, x2, y2] = clip_line( 345 xtemp, 346 ytemp, 347 xend, 348 yend, 349 width, 350 height, 351 ) 352 if x1 or x2 or y1 or y2: 353 bndpoints.append(QgsPointXY(x1 + extent.xMinimum(), 354 y1 + extent.yMinimum())) 355 bndpoints.append(QgsPointXY(x2 + extent.xMinimum(), 356 y2 + extent.yMinimum())) 357 if 0 in (x1, x2): 358 hasXMin = True 359 if XMinLine != edge[0]: 360 XMinNumber = XMinNumber + 1 361 XMinLine = edge[0] 362 if 0 in (y1, y2): 363 hasYMin = True 364 if YMinLine != edge[0]: 365 YMinNumber = YMinNumber + 1 366 YMinLine = edge[0] 367 if height in (y1, y2): 368 hasYMax = True 369 if YMaxLine != edge[0]: 370 YMaxNumber = YMaxNumber + 1 371 YMaxLine = edge[0] 372 if width in (x1, x2): 373 hasXMax = True 374 if XMaxLine != edge[0]: 375 XMaxNumber = XMaxNumber + 1 376 XMaxLine = edge[0] 377 378 # Add auxiliary points for corner cases, if necessary (duplicate 379 # points is not a problem - will be ignored later). 380 # a) Extreme input points (lowest, leftmost, rightmost, highest) 381 # A point can be extreme on both axis 382 if pt_x == xmin: # leftmost point 383 if XMinNumber == 0: 384 bndpoints.append(QgsPointXY(extent.xMinimum(), 385 extent.yMinimum())) 386 bndpoints.append(QgsPointXY(extent.xMinimum(), 387 height + extent.yMinimum())) 388 elif XMinNumber == 1: 389 if hasYMin: 390 bndpoints.append(QgsPointXY(extent.xMinimum(), 391 extent.yMinimum())) 392 elif hasYMax: 393 bndpoints.append(QgsPointXY(extent.xMinimum(), 394 height + extent.yMinimum())) 395 elif pt_x == xmax: # rightmost point 396 if XMaxNumber == 0: 397 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 398 extent.yMinimum())) 399 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 400 height + extent.yMinimum())) 401 elif XMaxNumber == 1: 402 if hasYMin: 403 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 404 extent.yMinimum())) 405 elif hasYMax: 406 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 407 height + extent.yMinimum())) 408 if pt_y == ymin: # lowest point 409 if YMinNumber == 0: 410 bndpoints.append(QgsPointXY(extent.xMinimum(), 411 extent.yMinimum())) 412 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 413 extent.yMinimum())) 414 elif YMinNumber == 1: 415 if hasXMin: 416 bndpoints.append(QgsPointXY(extent.xMinimum(), 417 extent.yMinimum())) 418 elif hasXMax: 419 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 420 extent.yMinimum())) 421 elif pt_y == ymax: # highest point 422 if YMaxNumber == 0: 423 bndpoints.append(QgsPointXY(extent.xMinimum(), 424 height + extent.yMinimum())) 425 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 426 height + extent.yMinimum())) 427 elif YMaxNumber == 1: 428 if hasXMin: 429 bndpoints.append(QgsPointXY(extent.xMinimum(), 430 height + extent.yMinimum())) 431 elif hasXMax: 432 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 433 height + extent.yMinimum())) 434 # b) Polygon that covers the x or the y extent: 435 if hasYMin and hasYMax: 436 if YMaxNumber > 1: 437 if hasXMin: 438 bndpoints.append(QgsPointXY(extent.xMinimum(), 439 extent.yMinimum())) 440 elif hasXMax: 441 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 442 extent.yMinimum())) 443 elif YMinNumber > 1: 444 if hasXMin: 445 bndpoints.append(QgsPointXY(extent.xMinimum(), 446 height + extent.yMinimum())) 447 elif hasXMax: 448 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 449 height + extent.yMinimum())) 450 elif hasXMin and hasXMax: 451 if XMaxNumber > 1: 452 if hasYMin: 453 bndpoints.append(QgsPointXY(extent.xMinimum(), 454 extent.yMinimum())) 455 elif hasYMax: 456 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 457 extent.yMinimum())) 458 elif XMinNumber > 1: 459 if hasYMin: 460 bndpoints.append(QgsPointXY(extent.xMinimum(), 461 height + extent.yMinimum())) 462 elif hasYMax: 463 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 464 height + extent.yMinimum())) 465 # c) Simple corners: 466 if XMinNumber == 1 and YMinNumber == 1 and not hasXMax and not hasYMax: 467 bndpoints.append(QgsPointXY(extent.xMinimum(), 468 extent.yMinimum())) 469 if XMinNumber == 1 and YMaxNumber == 1 and not hasXMax and not hasYMin: 470 bndpoints.append(QgsPointXY(extent.xMinimum(), 471 height + extent.yMinimum())) 472 if XMaxNumber == 1 and YMinNumber == 1 and not hasXMin and not hasYMax: 473 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 474 extent.yMinimum())) 475 if XMaxNumber == 1 and YMaxNumber == 1 and not hasXMin and not hasYMin: 476 bndpoints.append(QgsPointXY(width + extent.xMinimum(), 477 height + extent.yMinimum())) 478 return bndpoints 479