QGIS/python/plugins/processing/algs/qgis/VoronoiPolygons.py
Nyall Dawson 8cad2a6e16 [processing] Fix random selection count parameter
Remove a bunch of manual "max" values for numeric parameters
where the maximum just represents a 'large number' and not a real
constraint, and let the default parameter max value handling kick in instead.

In the case of random selection the max value exceeded the possible
range for integers in spin boxes and broke the widget.

Fixes #20015
2018-10-05 10:48:17 +10:00

481 lines
20 KiB
Python
Executable File

# -*- coding: utf-8 -*-
"""
***************************************************************************
VoronoiPolygons.py
---------------------
Date : August 2012
Copyright : (C) 2012 by Victor Olaya
Email : volayaf 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. *
* *
***************************************************************************
"""
__authors__ = 'Victor Olaya, Håvard Tveite'
__date__ = 'August 2012'
__copyright__ = '(C) 2012, Victor Olaya'
# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'
import os
from qgis.PyQt.QtGui import QIcon
from qgis.core import (QgsApplication,
QgsFeatureRequest,
QgsFeatureSink,
QgsFeature,
QgsGeometry,
QgsPointXY,
QgsWkbTypes,
QgsProcessing,
QgsProcessingException,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterFeatureSink,
QgsProcessingParameterNumber)
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
from . import voronoi
pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]
class VoronoiPolygons(QgisAlgorithm):
INPUT = 'INPUT'
BUFFER = 'BUFFER'
OUTPUT = 'OUTPUT'
def icon(self):
return QgsApplication.getThemeIcon("/algorithms/mAlgorithmVoronoi.svg")
def svgIconPath(self):
return QgsApplication.iconPath("/algorithms/mAlgorithmVoronoi.svg")
def group(self):
return self.tr('Vector geometry')
def groupId(self):
return 'vectorgeometry'
def __init__(self):
super().__init__()
def initAlgorithm(self, config=None):
self.addParameter(
QgsProcessingParameterFeatureSource(
self.INPUT, self.tr('Input layer'),
[QgsProcessing.TypeVectorPoint]))
self.addParameter(
QgsProcessingParameterNumber(
self.BUFFER, self.tr('Buffer region (% of extent)'),
minValue=0.0, defaultValue=0.0))
self.addParameter(
QgsProcessingParameterFeatureSink(
self.OUTPUT, self.tr('Voronoi polygons'),
type=QgsProcessing.TypeVectorPolygon))
def name(self):
return 'voronoipolygons'
def displayName(self):
return self.tr('Voronoi polygons')
def processAlgorithm(self, parameters, context, feedback):
source = self.parameterAsSource(parameters, self.INPUT, context)
if source is None:
raise QgsProcessingException(
self.invalidSourceError(parameters, self.INPUT))
buf = self.parameterAsDouble(parameters, self.BUFFER, context)
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT,
context, source.fields(),
QgsWkbTypes.Polygon,
source.sourceCrs())
if sink is None:
raise QgsProcessingException(
self.invalidSinkError(parameters, self.OUTPUT))
outFeat = QgsFeature()
extent = source.sourceExtent()
extraX = extent.width() * (buf / 100.0)
# Adjust the extent
extent.setXMinimum(extent.xMinimum() - extraX)
extent.setXMaximum(extent.xMaximum() + extraX)
extraY = extent.height() * (buf / 100.0)
extent.setYMinimum(extent.yMinimum() - extraY)
extent.setYMaximum(extent.yMaximum() + extraY)
height = extent.height()
width = extent.width()
c = voronoi.Context()
pts = []
ptDict = {}
ptNdx = -1
# Find the minimum and maximum x and y for the input points
xmin = width
xmax = 0
ymin = height
ymax = 0
features = source.getFeatures()
total = 100.0 / source.featureCount() if source.featureCount() else 0
for current, inFeat in enumerate(features):
if feedback.isCanceled():
break
geom = inFeat.geometry()
point = geom.asPoint()
x = point.x() - extent.xMinimum()
y = point.y() - extent.yMinimum()
pts.append((x, y))
ptNdx += 1
ptDict[ptNdx] = inFeat.id()
if x < xmin:
xmin = x
if y < ymin:
ymin = y
if x > xmax:
xmax = x
if y > ymax:
ymax = y
feedback.setProgress(int(current * total))
if xmin == xmax or ymin == ymax:
raise QgsProcessingException('The extent of the input points is '
'not a polygon (all the points are '
'on a vertical or horizontal line) '
'- cannot make a Voronoi diagram!')
xyminmax = [xmin, ymin, xmax, ymax]
if len(pts) < 3:
raise QgsProcessingException(
self.tr('Input file should contain at least 3 points. Choose '
'another file and try again.'))
# Eliminate duplicate points
uniqueSet = set(item for item in pts)
ids = [pts.index(item) for item in uniqueSet]
sl = voronoi.SiteList([voronoi.Site(i[0], i[1], sitenum=j)
for (j, i) in enumerate(uniqueSet)])
voronoi.voronoi(sl, c)
if len(c.polygons) == 0:
raise QgsProcessingException(
self.tr('There were no polygons created.'))
inFeat = QgsFeature()
current = 0
total = 100.0 / len(c.polygons)
# Clip each of the generated "polygons"
for (site, edges) in list(c.polygons.items()):
if feedback.isCanceled():
break
request = QgsFeatureRequest().setFilterFid(ptDict[ids[site]])
inFeat = next(source.getFeatures(request))
boundarypoints = self.clip_voronoi(edges, c, width,
height, extent,
inFeat.geometry().asPoint(),
xyminmax)
ptgeom = QgsGeometry.fromMultiPointXY(boundarypoints)
geom = QgsGeometry(ptgeom.convexHull())
outFeat.setGeometry(geom)
outFeat.setAttributes(inFeat.attributes())
sink.addFeature(outFeat, QgsFeatureSink.FastInsert)
current += 1
feedback.setProgress(int(current * total))
return {self.OUTPUT: dest_id}
def clip_voronoi(self, edges, c, width, height, extent, point, xyminmax):
"""Clip voronoi function based on code written for Inkscape.
Copyright (C) 2010 Alvin Penner, penner@vaxxine.com
Clips one Thiessen polygon (convex polygon) to extent
"""
pt_x = point.x() - extent.xMinimum()
pt_y = point.y() - extent.yMinimum()
(xmin, ymin, xmax, ymax) = xyminmax
def isclose(a, b, rel_tol=1e-9, abs_tol=0.0):
return abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)
def clip_line(x1, y1, x2, y2, w, h):
if x1 < 0 and x2 < 0:
# Completely to the left
return [0, 0, 0, 0]
if x1 > w and x2 > w:
# Completely to the right
return [0, 0, 0, 0]
if y1 < 0 and y2 < 0:
# Completely below
return [0, 0, 0, 0]
if y1 > h and y2 > h:
# Completely above
return [0, 0, 0, 0]
# Clip on the left envelope boundary
if x1 < 0:
# First point to the left
y1 = (y1 * x2 - y2 * x1) / (x2 - x1)
x1 = 0
if x2 < 0:
# Last point to the left
y2 = (y1 * x2 - y2 * x1) / (x2 - x1)
x2 = 0
# Clip on the right envelope boundary
if x1 > w:
# First point to the right
y1 = y1 + (w - x1) * (y2 - y1) / (x2 - x1)
x1 = w
if x2 > w:
# Last point to the right
y2 = y1 + (w - x1) * (y2 - y1) / (x2 - x1)
x2 = w
if isclose(x1, x2) and isclose(y1, y2):
return [0, 0, 0, 0]
# Clip on the bottom envelope boundary
if y1 < 0:
# First point below
x1 = (x1 * y2 - x2 * y1) / (y2 - y1)
y1 = 0
if y2 < 0:
# Second point below
x2 = (x1 * y2 - x2 * y1) / (y2 - y1)
y2 = 0
# Clip on the top envelope boundary
if y1 > h:
# First point above
x1 = x1 + (h - y1) * (x2 - x1) / (y2 - y1)
y1 = h
if y2 > h:
# Second point above
x2 = x1 + (h - y1) * (x2 - x1) / (y2 - y1)
y2 = h
if isclose(x1, x2) and isclose(y1, y2):
return [0, 0, 0, 0]
return [x1, y1, x2, y2]
bndpoints = []
hasXMin = False
hasYMin = False
hasXMax = False
hasYMax = False
XMinNumber = 0
XMaxNumber = 0
YMinNumber = 0
YMaxNumber = 0
# The same line may appear twice for collinear input points,
# so have to remember which lines have contributed
XMinLine = -1
XMaxLine = -1
YMinLine = -1
YMaxLine = -1
for edge in edges:
if edge[1] >= 0 and edge[2] >= 0:
# Two vertices
[x1, y1, x2, y2] = clip_line(
c.vertices[edge[1]][0],
c.vertices[edge[1]][1],
c.vertices[edge[2]][0],
c.vertices[edge[2]][1],
width,
height
)
elif edge[1] >= 0:
# Only one (left) vertex
if c.lines[edge[0]][1] == 0:
# Vertical line
#xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
xtemp = c.vertices[edge[1]][0]
ytemp = 0 - 1
#if c.vertices[edge[1]][1] > height / 2:
# ytemp = height
#else:
# ytemp = 0
else:
# Create an end of the line at the right edge - OK
xtemp = width
ytemp = (c.lines[edge[0]][2] - width *
c.lines[edge[0]][0]) / c.lines[edge[0]][1]
[x1, y1, x2, y2] = clip_line(
c.vertices[edge[1]][0],
c.vertices[edge[1]][1],
xtemp,
ytemp,
width,
height
)
elif edge[2] >= 0:
# Only one (right) vertex
if c.lines[edge[0]][1] == 0:
# Vertical line
#xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
xtemp = c.vertices[edge[2]][0]
ytemp = height + 1
#if c.vertices[edge[2]][1] > height / 2:
# ytemp = height
#else:
# ytemp = 0.0
else:
# End the line at the left edge - OK
xtemp = 0.0
ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1]
[x1, y1, x2, y2] = clip_line(
xtemp,
ytemp,
c.vertices[edge[2]][0],
c.vertices[edge[2]][1],
width,
height,
)
else:
# No vertex, only a line
if c.lines[edge[0]][1] == 0:
# Vertical line - should not happen
xtemp = c.lines[edge[0]][2] / c.lines[edge[0]][0]
ytemp = 0.0
xend = xtemp
yend = height
else:
# End the line at both edges - ???
xtemp = 0.0
ytemp = c.lines[edge[0]][2] / c.lines[edge[0]][1]
xend = width
yend = (c.lines[edge[0]][2] - width *
c.lines[edge[0]][0]) / c.lines[edge[0]][1]
[x1, y1, x2, y2] = clip_line(
xtemp,
ytemp,
xend,
yend,
width,
height,
)
if x1 or x2 or y1 or y2:
bndpoints.append(QgsPointXY(x1 + extent.xMinimum(),
y1 + extent.yMinimum()))
bndpoints.append(QgsPointXY(x2 + extent.xMinimum(),
y2 + extent.yMinimum()))
if 0 in (x1, x2):
hasXMin = True
if XMinLine != edge[0]:
XMinNumber = XMinNumber + 1
XMinLine = edge[0]
if 0 in (y1, y2):
hasYMin = True
if YMinLine != edge[0]:
YMinNumber = YMinNumber + 1
YMinLine = edge[0]
if height in (y1, y2):
hasYMax = True
if YMaxLine != edge[0]:
YMaxNumber = YMaxNumber + 1
YMaxLine = edge[0]
if width in (x1, x2):
hasXMax = True
if XMaxLine != edge[0]:
XMaxNumber = XMaxNumber + 1
XMaxLine = edge[0]
# Add auxiliary points for corner cases, if necessary (duplicate
# points is not a problem - will be ignored later).
# a) Extreme input points (lowest, leftmost, rightmost, highest)
# A point can be extreme on both axis
if pt_x == xmin: # leftmost point
if XMinNumber == 0:
bndpoints.append(QgsPointXY(extent.xMinimum(),
extent.yMinimum()))
bndpoints.append(QgsPointXY(extent.xMinimum(),
height + extent.yMinimum()))
elif XMinNumber == 1:
if hasYMin:
bndpoints.append(QgsPointXY(extent.xMinimum(),
extent.yMinimum()))
elif hasYMax:
bndpoints.append(QgsPointXY(extent.xMinimum(),
height + extent.yMinimum()))
elif pt_x == xmax: # rightmost point
if XMaxNumber == 0:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
extent.yMinimum()))
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
height + extent.yMinimum()))
elif XMaxNumber == 1:
if hasYMin:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
extent.yMinimum()))
elif hasYMax:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
height + extent.yMinimum()))
if pt_y == ymin: # lowest point
if YMinNumber == 0:
bndpoints.append(QgsPointXY(extent.xMinimum(),
extent.yMinimum()))
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
extent.yMinimum()))
elif YMinNumber == 1:
if hasXMin:
bndpoints.append(QgsPointXY(extent.xMinimum(),
extent.yMinimum()))
elif hasXMax:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
extent.yMinimum()))
elif pt_y == ymax: # highest point
if YMaxNumber == 0:
bndpoints.append(QgsPointXY(extent.xMinimum(),
height + extent.yMinimum()))
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
height + extent.yMinimum()))
elif YMaxNumber == 1:
if hasXMin:
bndpoints.append(QgsPointXY(extent.xMinimum(),
height + extent.yMinimum()))
elif hasXMax:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
height + extent.yMinimum()))
# b) Polygon that covers the x or the y extent:
if hasYMin and hasYMax:
if YMaxNumber > 1:
if hasXMin:
bndpoints.append(QgsPointXY(extent.xMinimum(),
extent.yMinimum()))
elif hasXMax:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
extent.yMinimum()))
elif YMinNumber > 1:
if hasXMin:
bndpoints.append(QgsPointXY(extent.xMinimum(),
height + extent.yMinimum()))
elif hasXMax:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
height + extent.yMinimum()))
elif hasXMin and hasXMax:
if XMaxNumber > 1:
if hasYMin:
bndpoints.append(QgsPointXY(extent.xMinimum(),
extent.yMinimum()))
elif hasYMax:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
extent.yMinimum()))
elif XMinNumber > 1:
if hasYMin:
bndpoints.append(QgsPointXY(extent.xMinimum(),
height + extent.yMinimum()))
elif hasYMax:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
height + extent.yMinimum()))
# c) Simple corners:
if XMinNumber == 1 and YMinNumber == 1 and not hasXMax and not hasYMax:
bndpoints.append(QgsPointXY(extent.xMinimum(),
extent.yMinimum()))
if XMinNumber == 1 and YMaxNumber == 1 and not hasXMax and not hasYMin:
bndpoints.append(QgsPointXY(extent.xMinimum(),
height + extent.yMinimum()))
if XMaxNumber == 1 and YMinNumber == 1 and not hasXMin and not hasYMax:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
extent.yMinimum()))
if XMaxNumber == 1 and YMaxNumber == 1 and not hasXMin and not hasYMin:
bndpoints.append(QgsPointXY(width + extent.xMinimum(),
height + extent.yMinimum()))
return bndpoints