QGIS API Documentation 3.38.0-Grenoble (exported)
Loading...
Searching...
No Matches
qgsalgorithmvoronoipolygons.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmvoronoipolygons.cpp
3 ---------------------
4 begin : July 2023
5 copyright : (C) 2023 by Alexander Bruy
6 email : alexander dot bruy at gmail dot com
7 ***************************************************************************/
8
9/***************************************************************************
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 ***************************************************************************/
17
19#include "qgsspatialindex.h"
20#include "qgsmultipoint.h"
21#include "qgsgeometryengine.h"
22
24
25QString QgsVoronoiPolygonsAlgorithm::name() const
26{
27 return QStringLiteral( "voronoipolygons" );
28}
29
30QString QgsVoronoiPolygonsAlgorithm::displayName() const
31{
32 return QObject::tr( "Voronoi polygons" );
33}
34
35QStringList QgsVoronoiPolygonsAlgorithm::tags() const
36{
37 return QObject::tr( "voronoi,polygons,tessellation,diagram" ).split( ',' );
38}
39
40QString QgsVoronoiPolygonsAlgorithm::group() const
41{
42 return QObject::tr( "Vector geometry" );
43}
44
45QString QgsVoronoiPolygonsAlgorithm::groupId() const
46{
47 return QStringLiteral( "vectorgeometry" );
48}
49
50QString QgsVoronoiPolygonsAlgorithm::shortHelpString() const
51{
52 return QObject::tr( "Generates a polygon layer containing the Voronoi diagram corresponding to input points." );
53}
54
55QgsVoronoiPolygonsAlgorithm *QgsVoronoiPolygonsAlgorithm::createInstance() const
56{
57 return new QgsVoronoiPolygonsAlgorithm();
58}
59
60void QgsVoronoiPolygonsAlgorithm::initAlgorithm( const QVariantMap & )
61{
62 addParameter( new QgsProcessingParameterFeatureSource( QStringLiteral( "INPUT" ), QObject::tr( "Input layer" ), QList< int >() << static_cast< int >( Qgis::ProcessingSourceType::VectorPoint ) ) );
63 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "BUFFER" ), QObject::tr( "Buffer region (% of extent)" ), Qgis::ProcessingNumberParameterType::Double, 0, false, 0 ) );
64 addParameter( new QgsProcessingParameterNumber( QStringLiteral( "TOLERANCE" ), QObject::tr( "Tolerance" ), Qgis::ProcessingNumberParameterType::Double, 0, true, 0 ) );
65 addParameter( new QgsProcessingParameterBoolean( QStringLiteral( "COPY_ATTRIBUTES" ), QObject::tr( "Copy attributes from input features" ), true ) );
66 addParameter( new QgsProcessingParameterFeatureSink( QStringLiteral( "OUTPUT" ), QObject::tr( "Voronoi polygons" ), Qgis::ProcessingSourceType::VectorPolygon ) );
67}
68
69bool QgsVoronoiPolygonsAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
70{
71 Q_UNUSED( feedback );
72
73 mSource.reset( parameterAsSource( parameters, QStringLiteral( "INPUT" ), context ) );
74 if ( !mSource )
75 throw QgsProcessingException( invalidSourceError( parameters, QStringLiteral( "INPUT" ) ) );
76
77 if ( mSource->featureCount() < 3 )
78 throw QgsProcessingException( QObject::tr( "Input layer should contain at least 3 points." ) );
79
80 mBuffer = parameterAsDouble( parameters, QStringLiteral( "BUFFER" ), context );
81 mTolerance = parameterAsDouble( parameters, QStringLiteral( "TOLERANCE" ), context );
82 mCopyAttributes = parameterAsBool( parameters, QStringLiteral( "COPY_ATTRIBUTES" ), context );
83
84 return true;
85}
86
87QVariantMap QgsVoronoiPolygonsAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
88{
89 QString dest;
90 if ( mCopyAttributes )
91 {
92 dest = voronoiWithAttributes( parameters, context, feedback );
93 }
94 else
95 {
96 dest = voronoiWithoutAttributes( parameters, context, feedback );
97 }
98
99 QVariantMap outputs;
100 outputs.insert( QStringLiteral( "OUTPUT" ), dest );
101 return outputs;
102}
103
104QString QgsVoronoiPolygonsAlgorithm::voronoiWithAttributes( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
105{
106 const QgsFields fields = mSource->fields();
107
108 QString dest;
109 std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, Qgis::WkbType::Polygon, mSource->sourceCrs() ) );
110 if ( !sink )
111 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
112
113 QgsFeatureRequest request;
115
116 QgsGeometry allPoints;
117 QHash< QgsFeatureId, QgsAttributes > attributeCache;
118
119 long long i = 0;
120 const double step = mSource->featureCount() > 0 ? 50.0 / mSource->featureCount() : 1;
121
122 const QgsSpatialIndex index( it, [&]( const QgsFeature & f )->bool
123 {
124 i++;
125 if ( feedback->isCanceled() )
126 return false;
127
128 feedback->setProgress( i * step );
129
130 if ( !f.hasGeometry() )
131 return true;
132
133 const QgsAbstractGeometry *geom = f.geometry().constGet();
134 if ( QgsWkbTypes::isMultiType( geom->wkbType() ) )
135 {
136 const QgsMultiPoint mp( *qgsgeometry_cast< const QgsMultiPoint * >( geom ) );
137 for ( auto pit = mp.const_parts_begin(); pit != mp.const_parts_end(); ++pit )
138 {
139 allPoints.addPartV2( qgsgeometry_cast< QgsPoint * >( *pit )->clone(), Qgis::WkbType::Point );
140 }
141 }
142 else
143 {
144 allPoints.addPartV2( qgsgeometry_cast< QgsPoint * >( geom )->clone(), Qgis::WkbType::Point );
145 }
146
147 attributeCache.insert( f.id(), f.attributes() );
148
149 return true;
151
152 QgsRectangle extent = mSource->sourceExtent();
153 double delta = extent.width() * mBuffer / 100.0;
154 extent.setXMinimum( extent.xMinimum() - delta );
155 extent.setXMaximum( extent.xMaximum() + delta );
156 delta = extent.height() * mBuffer / 100.0;
157 extent.setYMinimum( extent.yMinimum() - delta );
158 extent.setYMaximum( extent.yMaximum() + delta );
159 const QgsGeometry clippingGeom = QgsGeometry::fromRect( extent );
160
161 const QgsGeometry voronoiDiagram = allPoints.voronoiDiagram( clippingGeom, mTolerance );
162
163 if ( !voronoiDiagram.isEmpty() )
164 {
165 std::unique_ptr< QgsGeometryEngine > engine;
166 std::unique_ptr< QgsGeometryEngine > extentEngine( QgsGeometry::createGeometryEngine( clippingGeom.constGet() ) );
167 const QVector< QgsGeometry > collection = voronoiDiagram.asGeometryCollection();
168 for ( int i = 0; i < collection.length(); i++ )
169 {
170 if ( feedback->isCanceled() )
171 {
172 break;
173 }
174 QgsFeature f;
175 f.setFields( fields );
176 f.setGeometry( QgsGeometry( extentEngine->intersection( collection[i].constGet() ) ) );
177 const QList< QgsFeatureId > intersected = index.intersects( collection[i].boundingBox() );
178 engine.reset( QgsGeometry::createGeometryEngine( collection[i].constGet() ) );
179 engine->prepareGeometry();
180 for ( const QgsFeatureId id : intersected )
181 {
182 if ( engine->intersects( index.geometry( id ).constGet() ) )
183 {
184 f.setAttributes( attributeCache.value( id ) );
185 break;
186 }
187 }
188 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
189 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
190 feedback->setProgress( 50 + i * step );
191 }
192 }
193
194 return dest;
195}
196
197QString QgsVoronoiPolygonsAlgorithm::voronoiWithoutAttributes( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
198{
199 QgsFields fields;
200 fields.append( QgsField( QStringLiteral( "id" ), QMetaType::Type::LongLong ) );
201
202 QString dest;
203 std::unique_ptr< QgsFeatureSink > sink( parameterAsSink( parameters, QStringLiteral( "OUTPUT" ), context, dest, fields, Qgis::WkbType::Polygon, mSource->sourceCrs() ) );
204 if ( !sink )
205 throw QgsProcessingException( invalidSinkError( parameters, QStringLiteral( "OUTPUT" ) ) );
206
207 std::unique_ptr< QgsMultiPoint > points = std::make_unique< QgsMultiPoint >();
208
209 long long i = 0;
210 const double step = mSource->featureCount() > 0 ? 50.0 / mSource->featureCount() : 1;
212 QgsFeature f;
213 while ( it.nextFeature( f ) )
214 {
215 i++;
216
217 if ( feedback->isCanceled() )
218 break;
219
220 feedback->setProgress( i * step );
221
222 if ( !f.hasGeometry() )
223 continue;
224
225 const QgsAbstractGeometry *geom = f.geometry().constGet();
226 if ( QgsWkbTypes::isMultiType( geom->wkbType() ) )
227 {
228 const QgsMultiPoint mp( *qgsgeometry_cast< const QgsMultiPoint * >( geom ) );
229 for ( auto pit = mp.const_parts_begin(); pit != mp.const_parts_end(); ++pit )
230 {
231 points->addGeometry( qgsgeometry_cast< QgsPoint * >( *pit )->clone() );
232 }
233 }
234 else
235 {
236 points->addGeometry( qgsgeometry_cast< QgsPoint * >( geom )->clone() );
237 }
238 }
239
240 QgsRectangle extent = mSource->sourceExtent();
241 double delta = extent.width() * mBuffer / 100.0;
242 extent.setXMinimum( extent.xMinimum() - delta );
243 extent.setXMaximum( extent.xMaximum() + delta );
244 delta = extent.height() * mBuffer / 100.0;
245 extent.setYMinimum( extent.yMinimum() - delta );
246 extent.setYMaximum( extent.yMaximum() + delta );
247 const QgsGeometry clippingGeom = QgsGeometry::fromRect( extent );
248
249 QgsGeometry allPoints = QgsGeometry( std::move( points ) );
250 const QgsGeometry voronoiDiagram = allPoints.voronoiDiagram( clippingGeom, mTolerance );
251
252 if ( !voronoiDiagram.isEmpty() )
253 {
254 std::unique_ptr< QgsGeometryEngine > engine;
255 std::unique_ptr< QgsGeometryEngine > extentEngine( QgsGeometry::createGeometryEngine( clippingGeom.constGet() ) );
256 const QVector< QgsGeometry > collection = voronoiDiagram.asGeometryCollection();
257 for ( int i = 0; i < collection.length(); i++ )
258 {
259 if ( feedback->isCanceled() )
260 {
261 break;
262 }
263 QgsFeature f;
264 f.setFields( fields );
265 f.setGeometry( QgsGeometry( extentEngine->intersection( collection[i].constGet() ) ) );
266 f.setAttributes( QgsAttributes() << i );
267 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
268 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
269 feedback->setProgress( i * step );
270 }
271 }
272
273 return dest;
274}
275
@ VectorPoint
Vector point layers.
@ VectorPolygon
Vector polygon layers.
@ SkipGeometryValidityChecks
Invalid geometry checks should always be skipped. This flag can be useful for algorithms which always...
@ Polygon
Polygon.
Abstract base class for all geometries.
Qgis::WkbType wkbType() const
Returns the WKB type of the geometry.
A vector of attributes.
Wrapper for iterator of features from vector data provider or vector layer.
bool nextFeature(QgsFeature &f)
Fetch next feature and stores in f, returns true on success.
This class wraps a request for features to a vector layer (or directly its vector data provider).
@ FastInsert
Use faster inserts, at the cost of updating the passed features to reflect changes made at the provid...
The feature class encapsulates a single feature including its unique ID, geometry and a list of field...
Definition qgsfeature.h:58
QgsAttributes attributes
Definition qgsfeature.h:67
QgsFeatureId id
Definition qgsfeature.h:66
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
void setFields(const QgsFields &fields, bool initAttributes=false)
Assigns a field map with the feature to allow attribute access by attribute name.
QgsGeometry geometry
Definition qgsfeature.h:69
bool hasGeometry() const
Returns true if the feature has an associated geometry.
void setGeometry(const QgsGeometry &geometry)
Set the feature's geometry.
bool isCanceled() const
Tells whether the operation has been canceled already.
Definition qgsfeedback.h:53
void setProgress(double progress)
Sets the current progress for the feedback object.
Definition qgsfeedback.h:61
Encapsulate a field in an attribute table or data source.
Definition qgsfield.h:53
Container of fields for a vector layer.
Definition qgsfields.h:46
bool append(const QgsField &field, Qgis::FieldOrigin origin=Qgis::FieldOrigin::Provider, int originIndex=-1)
Appends a field.
Definition qgsfields.cpp:60
A geometry is the spatial representation of a feature.
static QgsGeometry fromRect(const QgsRectangle &rect)
Creates a new geometry from a QgsRectangle.
QVector< QgsGeometry > asGeometryCollection() const
Returns contents of the geometry as a list of geometries.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
QgsGeometry voronoiDiagram(const QgsGeometry &extent=QgsGeometry(), double tolerance=0.0, bool edgesOnly=false) const
Creates a Voronoi diagram for the nodes contained within the geometry.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
Qgis::GeometryOperationResult addPartV2(const QVector< QgsPointXY > &points, Qgis::WkbType wkbType=Qgis::WkbType::Unknown)
Adds a new part to a the geometry.
static QgsGeometryEngine * createGeometryEngine(const QgsAbstractGeometry *geometry, double precision=0.0)
Creates and returns a new geometry engine representing the specified geometry using precision on a gr...
Multi point geometry collection.
Contains information about the context in which a processing algorithm is executed.
Custom exception class for processing related exceptions.
Base class for providing feedback from a processing algorithm.
A boolean parameter for processing algorithms.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
A numeric parameter for processing algorithms.
A rectangle specified with double values.
double xMinimum() const
Returns the x minimum value (left side of rectangle).
void setYMinimum(double y)
Set the minimum y value.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
void setXMinimum(double x)
Set the minimum x value.
double width() const
Returns the width of the rectangle.
double xMaximum() const
Returns the x maximum value (right side of rectangle).
double yMaximum() const
Returns the y maximum value (top side of rectangle).
void setYMaximum(double y)
Set the maximum y value.
void setXMaximum(double x)
Set the maximum x value.
double height() const
Returns the height of the rectangle.
A spatial index for QgsFeature objects.
@ FlagStoreFeatureGeometries
Indicates that the spatial index should also store feature geometries. This requires more memory,...
static bool isMultiType(Qgis::WkbType type)
Returns true if the WKB type is a multi type.
qint64 QgsFeatureId
64 bit feature ids negative numbers are used for uncommitted/newly added features