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