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