QGIS API Documentation 3.43.0-Master (e01d6d7c4c0)
qgsalgorithmrandompointsinpolygons.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsalgorithmrandompointsinpolygons.cpp
3 ---------------------
4 begin : March 2020
5 copyright : (C) 2020 by HÃ¥vard Tveite
6 email : havard dot tveite at nmbu dot no
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
18
20#include "qgsspatialindex.h"
21
22#include <random>
23
25
26QString QgsRandomPointsInPolygonsAlgorithm::name() const
27{
28 return QStringLiteral( "randompointsinpolygons" );
29}
30
31QString QgsRandomPointsInPolygonsAlgorithm::displayName() const
32{
33 return QObject::tr( "Random points in polygons" );
34}
35
36QStringList QgsRandomPointsInPolygonsAlgorithm::tags() const
37{
38 return QObject::tr( "seed,attributes,create" ).split( ',' );
39}
40
41QString QgsRandomPointsInPolygonsAlgorithm::group() const
42{
43 return QObject::tr( "Vector creation" );
44}
45
46QString QgsRandomPointsInPolygonsAlgorithm::groupId() const
47{
48 return QStringLiteral( "vectorcreation" );
49}
50
51void QgsRandomPointsInPolygonsAlgorithm::initAlgorithm( const QVariantMap & )
52{
53 addParameter( new QgsProcessingParameterFeatureSource( INPUT, QObject::tr( "Input polygon layer" ), QList<int>() << static_cast<int>( Qgis::ProcessingSourceType::VectorPolygon ) ) );
54 auto numberPointsParam = std::make_unique<QgsProcessingParameterNumber>( POINTS_NUMBER, QObject::tr( "Number of points for each feature" ), Qgis::ProcessingNumberParameterType::Integer, 1, false, 1 );
55 numberPointsParam->setIsDynamic( true );
56 numberPointsParam->setDynamicPropertyDefinition( QgsPropertyDefinition( POINTS_NUMBER, QObject::tr( "Number of points for each feature" ), QgsPropertyDefinition::IntegerPositive ) );
57 numberPointsParam->setDynamicLayerParameterName( QStringLiteral( "INPUT" ) );
58 addParameter( numberPointsParam.release() );
59
60 auto minDistParam = std::make_unique<QgsProcessingParameterDistance>( MIN_DISTANCE, QObject::tr( "Minimum distance between points" ), 0, INPUT, true, 0 );
61 minDistParam->setIsDynamic( true );
62 minDistParam->setDynamicPropertyDefinition( QgsPropertyDefinition( MIN_DISTANCE, QObject::tr( "Minimum distance between points" ), QgsPropertyDefinition::DoublePositive ) );
63 minDistParam->setDynamicLayerParameterName( QStringLiteral( "INPUT" ) );
64 addParameter( minDistParam.release() );
65
66 auto minDistGlobalParam = std::make_unique<QgsProcessingParameterDistance>( MIN_DISTANCE_GLOBAL, QObject::tr( "Global minimum distance between points" ), 0, INPUT, true, 0 );
67 minDistGlobalParam->setFlags( minDistGlobalParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
68 addParameter( minDistGlobalParam.release() );
69
70 auto maxAttemptsParam = std::make_unique<QgsProcessingParameterNumber>( MAX_TRIES_PER_POINT, QObject::tr( "Maximum number of search attempts (for Min. dist. > 0)" ), Qgis::ProcessingNumberParameterType::Integer, 10, true, 1 );
71 maxAttemptsParam->setFlags( maxAttemptsParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
72 maxAttemptsParam->setIsDynamic( true );
73 maxAttemptsParam->setDynamicPropertyDefinition( QgsPropertyDefinition( MAX_TRIES_PER_POINT, QObject::tr( "Maximum number of attempts per point (for Min. dist. > 0)" ), QgsPropertyDefinition::IntegerPositiveGreaterZero ) );
74 maxAttemptsParam->setDynamicLayerParameterName( QStringLiteral( "INPUT" ) );
75 addParameter( maxAttemptsParam.release() );
76
77 auto randomSeedParam = std::make_unique<QgsProcessingParameterNumber>( SEED, QObject::tr( "Random seed" ), Qgis::ProcessingNumberParameterType::Integer, QVariant(), true, 1 );
78 randomSeedParam->setFlags( randomSeedParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
79 addParameter( randomSeedParam.release() );
80
81 auto includePolygonAttrParam = std::make_unique<QgsProcessingParameterBoolean>( INCLUDE_POLYGON_ATTRIBUTES, QObject::tr( "Include polygon attributes" ), true );
82 includePolygonAttrParam->setFlags( includePolygonAttrParam->flags() | Qgis::ProcessingParameterFlag::Advanced );
83 addParameter( includePolygonAttrParam.release() );
84
85 addParameter( new QgsProcessingParameterFeatureSink( OUTPUT, QObject::tr( "Random points in polygons" ), Qgis::ProcessingSourceType::VectorPoint ) );
86
87 addOutput( new QgsProcessingOutputNumber( OUTPUT_POINTS, QObject::tr( "Total number of points generated" ) ) );
88 addOutput( new QgsProcessingOutputNumber( POINTS_MISSED, QObject::tr( "Number of missed points" ) ) );
89 addOutput( new QgsProcessingOutputNumber( POLYGONS_WITH_MISSED_POINTS, QObject::tr( "Number of polygons with missed points" ) ) );
90 addOutput( new QgsProcessingOutputNumber( FEATURES_WITH_EMPTY_OR_NO_GEOMETRY, QObject::tr( "Number of features with empty or no geometry" ) ) );
91}
92
93QString QgsRandomPointsInPolygonsAlgorithm::shortDescription() const
94{
95 return QObject::tr( "Creates a layer with a number of points placed randomly in each polygon of a given layer." );
96}
97
98QString QgsRandomPointsInPolygonsAlgorithm::shortHelpString() const
99{
100 return QObject::tr( "<p>This algorithm creates a point layer, with points placed randomly "
101 "in the polygons of the <i><b>Input polygon layer</b></i>.</p> "
102 "<ul><li>For each feature in the <i><b>Input polygon layer</b></i>, the algorithm attempts to add "
103 "the specified <i><b>Number of points for each feature</b></i> to the output layer.</li> "
104 "<li>A <i><b>Minimum distance between points</b></i> and a "
105 "<i><b>Global minimum distance between points</b></i> can be specified.<br> "
106 "A point will not be added if there is an already generated point within "
107 "this (Euclidean) distance from the generated location. "
108 "With <i>Minimum distance between points</i>, only points in the same "
109 "polygon feature are considered, while for <i>Global minimum distance "
110 "between points</i> all previously generated points are considered. "
111 "If the <i>Global minimum distance between points</i> is set equal to "
112 "or larger than the (local) <i>Minimum distance between points</i>, the "
113 "latter has no effect.<br> "
114 "If the <i>Minimum distance between points</i> is too large, "
115 "it may not be possible to generate the specified <i>Number of points "
116 "for each feature</i>, but all the generated points are returned.</li> "
117 "<li>The <i><b>Maximum number of attempts per point</b></i> can be specified.</li> "
118 "<li>The seed for the random generator can be provided (<b><i>Random seed</i></b> "
119 "- integer, greater than 0).</li> "
120 "<li>The user can choose not to <i><b>Include polygon feature attributes</b></i> in "
121 "the attributes of the generated point features.</li> "
122 "</ul> "
123 "The total number of points will be<br> <b>'number of input features'</b> * "
124 "<i><b>Number of points for each feature</b></i><br> if there are no misses. "
125 "The <i>Number of points for each feature</i>, <i>Minimum distance between points</i> "
126 "and <i>Maximum number of attempts per point</i> can be data defined. "
127 "<p>Output from the algorithm:</p> "
128 "<ul> "
129 "<li> The number of features with an empty or no geometry "
130 "(<code>FEATURES_WITH_EMPTY_OR_NO_GEOMETRY</code>).</li> "
131 "<li> A point layer containing the random points (<code>OUTPUT</code>).</li> "
132 "<li> The number of generated features (<code>OUTPUT_POINTS</code>).</li> "
133 "<li> The number of missed points (<code>POINTS_MISSED</code>).</li> "
134 "<li> The number of features with non-empty geometry and missing points "
135 "(<code>POLYGONS_WITH_MISSED_POINTS</code>).</li> "
136 "</ul>"
137 );
138}
139
140QgsRandomPointsInPolygonsAlgorithm *QgsRandomPointsInPolygonsAlgorithm::createInstance() const
141{
142 return new QgsRandomPointsInPolygonsAlgorithm();
143}
144
145bool QgsRandomPointsInPolygonsAlgorithm::prepareAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback * )
146{
147 mNumPoints = parameterAsInt( parameters, POINTS_NUMBER, context );
148 mDynamicNumPoints = QgsProcessingParameters::isDynamic( parameters, POINTS_NUMBER );
149 if ( mDynamicNumPoints )
150 mNumPointsProperty = parameters.value( POINTS_NUMBER ).value<QgsProperty>();
151
152 mMinDistance = parameterAsDouble( parameters, MIN_DISTANCE, context );
153 mDynamicMinDistance = QgsProcessingParameters::isDynamic( parameters, MIN_DISTANCE );
154 if ( mDynamicMinDistance )
155 mMinDistanceProperty = parameters.value( MIN_DISTANCE ).value<QgsProperty>();
156
157 mMaxAttempts = parameterAsInt( parameters, MAX_TRIES_PER_POINT, context );
158 mDynamicMaxAttempts = QgsProcessingParameters::isDynamic( parameters, MAX_TRIES_PER_POINT );
159 if ( mDynamicMaxAttempts )
160 mMaxAttemptsProperty = parameters.value( MAX_TRIES_PER_POINT ).value<QgsProperty>();
161
162 mMinDistanceGlobal = parameterAsDouble( parameters, MIN_DISTANCE_GLOBAL, context );
163
164 mUseRandomSeed = parameters.value( SEED ).isValid();
165 mRandSeed = parameterAsInt( parameters, SEED, context );
166 mIncludePolygonAttr = parameterAsBoolean( parameters, INCLUDE_POLYGON_ATTRIBUTES, context );
167 return true;
168}
169
170QVariantMap QgsRandomPointsInPolygonsAlgorithm::processAlgorithm( const QVariantMap &parameters, QgsProcessingContext &context, QgsProcessingFeedback *feedback )
171{
172 std::unique_ptr<QgsProcessingFeatureSource> polygonSource( parameterAsSource( parameters, INPUT, context ) );
173 if ( !polygonSource )
174 throw QgsProcessingException( invalidSourceError( parameters, INPUT ) );
175
176 QgsFields fields;
177 fields.append( QgsField( QStringLiteral( "rand_point_id" ), QMetaType::Type::LongLong ) );
178 if ( mIncludePolygonAttr )
179 fields.extend( polygonSource->fields() );
180
181 QString ldest;
182 std::unique_ptr<QgsFeatureSink> sink( parameterAsSink( parameters, OUTPUT, context, ldest, fields, Qgis::WkbType::Point, polygonSource->sourceCrs() ) );
183 if ( !sink )
184 throw QgsProcessingException( invalidSinkError( parameters, OUTPUT ) );
185
186 QgsExpressionContext expressionContext = createExpressionContext( parameters, context, polygonSource.get() );
187
188 // Initialize random engine -- note that we only use this if the user has specified a fixed seed
189 std::random_device rd;
190 std::mt19937 mt( !mUseRandomSeed ? rd() : mRandSeed );
191 const std::uniform_real_distribution<> uniformDist( 0, 1 );
192 std::uniform_int_distribution<> uniformIntDist( 1, 999999999 );
193
194 // Index for finding global close points (mMinDistance != 0)
195 QgsSpatialIndex globalIndex;
196 int indexPoints = 0;
197
198 int totNPoints = 0;
199 int missedPoints = 0;
200 int missedPolygons = 0;
201 int emptyOrNullGeom = 0;
202
203 long long attempts = 0; // used for unique feature IDs in the indexes
204 const long numberOfFeatures = polygonSource->featureCount();
205 long long desiredNumberOfPoints = 0;
206 const double featureProgressStep = 100.0 / ( numberOfFeatures > 0 ? numberOfFeatures : 1 );
207 double baseFeatureProgress = 0.0;
208 QgsFeature polyFeat;
209 QgsFeatureIterator fitL = mIncludePolygonAttr || mDynamicNumPoints || mDynamicMinDistance || mDynamicMaxAttempts ? polygonSource->getFeatures()
210 : polygonSource->getFeatures( QgsFeatureRequest().setNoAttributes() );
211 while ( fitL.nextFeature( polyFeat ) )
212 {
213 if ( feedback->isCanceled() )
214 {
215 feedback->setProgress( 0 );
216 break;
217 }
218 if ( !polyFeat.hasGeometry() )
219 {
220 // Increment invalid features count
221 emptyOrNullGeom++;
222 baseFeatureProgress += featureProgressStep;
223 feedback->setProgress( baseFeatureProgress );
224 continue;
225 }
226 const QgsGeometry polyGeom( polyFeat.geometry() );
227 if ( polyGeom.isEmpty() )
228 {
229 // Increment invalid features count
230 emptyOrNullGeom++;
231 baseFeatureProgress += featureProgressStep;
232 feedback->setProgress( baseFeatureProgress );
233 continue;
234 }
235 if ( mDynamicNumPoints || mDynamicMinDistance || mDynamicMaxAttempts )
236 {
237 expressionContext.setFeature( polyFeat );
238 }
239 // (Re)initialize the local (per polygon) index
240 QgsSpatialIndex localIndex;
241 int localIndexPoints = 0;
242 int pointsAddedForThisFeature = 0;
243 // Get data defined parameters
244 int numberPointsForThisFeature = mNumPoints;
245 if ( mDynamicNumPoints )
246 {
247 numberPointsForThisFeature = mNumPointsProperty.valueAsInt( expressionContext, 0 );
248 }
249 desiredNumberOfPoints += numberPointsForThisFeature;
250 int maxAttemptsForThisFeature = mMaxAttempts;
251 if ( mDynamicMaxAttempts )
252 maxAttemptsForThisFeature = mMaxAttemptsProperty.valueAsInt( expressionContext, maxAttemptsForThisFeature );
253 double minDistanceForThisFeature = mMinDistance;
254 if ( mDynamicMinDistance )
255 minDistanceForThisFeature = mMinDistanceProperty.valueAsDouble( expressionContext, minDistanceForThisFeature );
256 const double pointProgressIncrement = featureProgressStep / ( numberPointsForThisFeature * maxAttemptsForThisFeature );
257 double pointProgress = 0.0;
258 // Check if we can avoid using the acceptPoint function
259 if ( ( minDistanceForThisFeature == 0 ) && ( mMinDistanceGlobal == 0 ) )
260 {
261 QVector<QgsPointXY> newPoints = polyGeom.randomPointsInPolygon( numberPointsForThisFeature, mUseRandomSeed ? uniformIntDist( mt ) : 0 );
262 for ( int i = 0; i < newPoints.length(); i++ )
263 {
264 // add the point
265 const QgsPointXY pt = newPoints[i];
266 QgsFeature f = QgsFeature( totNPoints );
267 QgsAttributes pAttrs = QgsAttributes();
268 pAttrs.append( totNPoints );
269 if ( mIncludePolygonAttr )
270 {
271 pAttrs.append( polyFeat.attributes() );
272 }
273 f.setAttributes( pAttrs );
274 const QgsGeometry newGeom = QgsGeometry::fromPointXY( pt );
275 f.setGeometry( newGeom );
276 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
277 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
278 totNPoints++;
279 pointsAddedForThisFeature++;
280 pointProgress += pointProgressIncrement * ( maxAttemptsForThisFeature );
281 }
282 feedback->setProgress( baseFeatureProgress + pointProgress );
283 continue;
284 }
285 else
286 {
287 // Have to check for minimum distance, provide the acceptPoints function
288 QVector<QgsPointXY> newPoints = polyGeom.randomPointsInPolygon( numberPointsForThisFeature, [&]( const QgsPointXY &newPoint ) -> bool {
289 attempts++;
290 // May have to check minimum distance to existing points
291 // The first point can always be added
292 // Local first (if larger than global)
293 if ( minDistanceForThisFeature != 0 && mMinDistanceGlobal < minDistanceForThisFeature && localIndexPoints > 0 )
294 {
295 const QList<QgsFeatureId> neighbors = localIndex.nearestNeighbor( newPoint, 1, minDistanceForThisFeature );
296 //if ( totNPoints > 0 && !neighbors.empty() )
297 if ( !neighbors.empty() )
298 {
299 return false;
300 }
301 }
302 // The global
303 if ( mMinDistanceGlobal != 0.0 && indexPoints > 0 )
304 {
305 const QList<QgsFeatureId> neighbors = globalIndex.nearestNeighbor( newPoint, 1, mMinDistanceGlobal );
306 //if ( totNPoints > 0 && !neighbors.empty() )
307 if ( !neighbors.empty() )
308 {
309 return false;
310 }
311 }
312 // Point is accepted - add it to the indexes
313 QgsFeature f = QgsFeature( attempts );
314 QgsAttributes pAttrs = QgsAttributes();
315 pAttrs.append( attempts );
316 f.setAttributes( pAttrs );
317 const QgsGeometry newGeom = QgsGeometry::fromPointXY( newPoint );
318
319 f.setGeometry( newGeom );
320 //totNPoints++;
321
322 if ( minDistanceForThisFeature != 0 )
323 {
324 if ( !localIndex.addFeature( f ) )
325 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QString() ) );
326 localIndexPoints++;
327 }
328 if ( mMinDistanceGlobal != 0.0 )
329 {
330 if ( !globalIndex.addFeature( f ) )
331 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QString() ) );
332 indexPoints++;
333 }
334 return true; }, mUseRandomSeed ? uniformIntDist( mt ) : 0, feedback, maxAttemptsForThisFeature );
335
336 // create and output features for the generated points
337 for ( int i = 0; i < newPoints.length(); i++ )
338 {
339 const QgsPointXY pt = newPoints[i];
340 QgsFeature f = QgsFeature( totNPoints );
341 QgsAttributes pAttrs = QgsAttributes();
342 pAttrs.append( totNPoints );
343 if ( mIncludePolygonAttr )
344 {
345 pAttrs.append( polyFeat.attributes() );
346 }
347 f.setAttributes( pAttrs );
348 const QgsGeometry newGeom = QgsGeometry::fromPointXY( pt );
349 f.setGeometry( newGeom );
350 if ( !sink->addFeature( f, QgsFeatureSink::FastInsert ) )
351 throw QgsProcessingException( writeFeatureError( sink.get(), parameters, QStringLiteral( "OUTPUT" ) ) );
352 totNPoints++;
353 pointsAddedForThisFeature++;
354 pointProgress += pointProgressIncrement * ( maxAttemptsForThisFeature );
355 }
356 feedback->setProgress( baseFeatureProgress + pointProgress );
357 }
358
359 baseFeatureProgress += featureProgressStep;
360 if ( pointsAddedForThisFeature < numberPointsForThisFeature )
361 {
362 missedPolygons++;
363 }
364 feedback->setProgress( baseFeatureProgress );
365 } // while features
366 missedPoints = desiredNumberOfPoints - totNPoints;
367 feedback->pushInfo( QObject::tr( "Total number of points generated: "
368 "%1\nNumber of missed points: "
369 "%2\nPolygons with missing points: "
370 "%3\nFeatures with empty or missing "
371 "geometries: %4"
372 )
373 .arg( totNPoints )
374 .arg( missedPoints )
375 .arg( missedPolygons )
376 .arg( emptyOrNullGeom ) );
377
378 sink->finalize();
379
380 QVariantMap outputs;
381 outputs.insert( OUTPUT, ldest );
382 outputs.insert( OUTPUT_POINTS, totNPoints );
383 outputs.insert( POINTS_MISSED, missedPoints );
384 outputs.insert( POLYGONS_WITH_MISSED_POINTS, missedPolygons );
385 outputs.insert( FEATURES_WITH_EMPTY_OR_NO_GEOMETRY, emptyOrNullGeom );
386
387 return outputs;
388}
389
@ VectorPoint
Vector point layers.
@ VectorPolygon
Vector polygon layers.
@ Advanced
Parameter is an advanced parameter which should be hidden from users by default.
A vector of attributes.
Expression contexts are used to encapsulate the parameters around which a QgsExpression should be eva...
void setFeature(const QgsFeature &feature)
Convenience function for setting a feature for the context.
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.
Wraps a request for features to a vector layer (or directly its vector data provider).
QgsFeatureRequest & setNoAttributes()
Set that no attributes will be fetched.
@ 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
void setAttributes(const QgsAttributes &attrs)
Sets the feature's attributes.
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:70
void extend(const QgsFields &other)
Extends with fields from another QgsFields container.
A geometry is the spatial representation of a feature.
static QgsGeometry fromPointXY(const QgsPointXY &point)
Creates a new geometry from a QgsPointXY object.
Represents a 2D point.
Definition qgspointxy.h:60
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.
virtual void pushInfo(const QString &info)
Pushes a general informational message from the algorithm.
A numeric output for processing algorithms.
A feature sink output for processing algorithms.
An input feature source (such as vector layers) parameter for processing algorithms.
static bool isDynamic(const QVariantMap &parameters, const QString &name)
Returns true if the parameter with matching name is a dynamic parameter, and must be evaluated once f...
Definition for a property.
Definition qgsproperty.h:45
@ IntegerPositiveGreaterZero
Non-zero positive integer values.
Definition qgsproperty.h:54
@ IntegerPositive
Positive integer values (including 0)
Definition qgsproperty.h:53
@ DoublePositive
Positive double value (including 0)
Definition qgsproperty.h:56
A store for object properties.
double valueAsDouble(const QgsExpressionContext &context, double defaultValue=0.0, bool *ok=nullptr) const
Calculates the current value of the property and interprets it as a double.
int valueAsInt(const QgsExpressionContext &context, int defaultValue=0, bool *ok=nullptr) const
Calculates the current value of the property and interprets it as an integer.
A spatial index for QgsFeature objects.
QList< QgsFeatureId > nearestNeighbor(const QgsPointXY &point, int neighbors=1, double maxDistance=0) const
Returns nearest neighbors to a point.
bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=QgsFeatureSink::Flags()) override
Adds a feature to the index.