mirror of
https://github.com/qgis/QGIS.git
synced 2025-02-28 00:17:30 -05:00
2628c480c5
Prevent changes to files that weren't changed between releases. This eases review of the changes between releases significantly.
477 lines
20 KiB
Python
Executable File
477 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'
|
|
|
|
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
|