QGIS API Documentation  2.99.0-Master (69af2f5)
qgspointsample.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgspointsample.cpp
3  ---------------------
4  begin : July 2013
5  copyright : (C) 2013 by Marco Hugentobler
6  email : marco dot hugentobler at sourcepole dot ch
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 #include "qgspointsample.h"
16 #include "qgsfeatureiterator.h"
17 #include "qgsgeometry.h"
18 #include "qgsspatialindex.h"
19 #include "qgsvectorfilewriter.h"
20 #include "qgsvectorlayer.h"
21 #include <QFile>
22 #include "mersenne-twister.h"
23 
24 
25 QgsPointSample::QgsPointSample( QgsVectorLayer *inputLayer, const QString &outputLayer, const QString &nPointsAttribute, const QString &minDistAttribute ): mInputLayer( inputLayer ),
26  mOutputLayer( outputLayer ), mNumberOfPointsAttribute( nPointsAttribute ), mMinDistanceAttribute( minDistAttribute ), mNCreatedPoints( 0 )
27 {
28 }
29 
31  : mInputLayer( nullptr )
32  , mNCreatedPoints( 0 )
33 {
34 }
35 
36 int QgsPointSample::createRandomPoints( QProgressDialog *pd )
37 {
38  Q_UNUSED( pd );
39 
40  //create input layer from id (test if polygon, valid)
41  if ( !mInputLayer )
42  {
43  return 1;
44  }
45 
46  if ( mInputLayer->geometryType() != QgsWkbTypes::PolygonGeometry )
47  {
48  return 2;
49  }
50 
51  //delete output file if it already exists
52  if ( QFile::exists( mOutputLayer ) )
53  {
55  }
56 
57  //create vector file writer
58  QgsFields outputFields;
59  outputFields.append( QgsField( QStringLiteral( "id" ), QVariant::Int ) );
60  outputFields.append( QgsField( QStringLiteral( "station_id" ), QVariant::Int ) );
61  outputFields.append( QgsField( QStringLiteral( "stratum_id" ), QVariant::Int ) );
62  QgsVectorFileWriter writer( mOutputLayer, QStringLiteral( "UTF-8" ),
63  outputFields,
65  mInputLayer->crs() );
66 
67  //check if creation of output layer successful
68  if ( writer.hasError() != QgsVectorFileWriter::NoError )
69  {
70  return 3;
71  }
72 
73  //init random number generator
74  mt_srand( QTime::currentTime().msec() );
75  QgsFeature fet;
76  int nPoints = 0;
77  double minDistance = 0;
78  mNCreatedPoints = 0;
79 
80  QgsFeatureIterator fIt = mInputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes(
81  QStringList() << mNumberOfPointsAttribute << mMinDistanceAttribute, mInputLayer->fields() ) );
82  while ( fIt.nextFeature( fet ) )
83  {
84  nPoints = fet.attribute( mNumberOfPointsAttribute ).toInt();
85  if ( !mMinDistanceAttribute.isEmpty() )
86  {
87  minDistance = fet.attribute( mMinDistanceAttribute ).toDouble();
88  }
89  addSamplePoints( fet, writer, nPoints, minDistance );
90  }
91 
92  return 0;
93 }
94 
95 void QgsPointSample::addSamplePoints( QgsFeature &inputFeature, QgsVectorFileWriter &writer, int nPoints, double minDistance )
96 {
97  if ( !inputFeature.hasGeometry() )
98  return;
99 
100  QgsGeometry geom = inputFeature.geometry();
101  QgsRectangle geomRect = geom.boundingBox();
102  if ( geomRect.isEmpty() )
103  {
104  return;
105  }
106 
107  QgsSpatialIndex sIndex; //to check minimum distance
108  QMap< QgsFeatureId, QgsPointXY > pointMapForFeature;
109 
110  int nIterations = 0;
111  int maxIterations = nPoints * 200;
112  int points = 0;
113 
114  double randX = 0;
115  double randY = 0;
116 
117  while ( nIterations < maxIterations && points < nPoints )
118  {
119  randX = ( ( double )mt_rand() / MD_RAND_MAX ) * geomRect.width() + geomRect.xMinimum();
120  randY = ( ( double )mt_rand() / MD_RAND_MAX ) * geomRect.height() + geomRect.yMinimum();
121  QgsPointXY randPoint( randX, randY );
122  QgsGeometry ptGeom = QgsGeometry::fromPoint( randPoint );
123  if ( ptGeom.within( geom ) && checkMinDistance( randPoint, sIndex, minDistance, pointMapForFeature ) )
124  {
125  //add feature to writer
126  QgsFeature f( mNCreatedPoints );
127  f.setAttribute( QStringLiteral( "id" ), mNCreatedPoints + 1 );
128  f.setAttribute( QStringLiteral( "station_id" ), points + 1 );
129  f.setAttribute( QStringLiteral( "stratum_id" ), inputFeature.id() );
130  f.setGeometry( ptGeom );
131  writer.addFeature( f );
132  sIndex.insertFeature( f );
133  pointMapForFeature.insert( mNCreatedPoints, randPoint );
134  ++points;
135  ++mNCreatedPoints;
136  }
137  ++nIterations;
138  }
139 }
140 
141 bool QgsPointSample::checkMinDistance( QgsPointXY &pt, QgsSpatialIndex &index, double minDistance, QMap< QgsFeatureId, QgsPointXY > &pointMap )
142 {
143  if ( minDistance <= 0 )
144  {
145  return true;
146  }
147 
148  QList<QgsFeatureId> neighborList = index.nearestNeighbor( pt, 1 );
149  if ( neighborList.isEmpty() )
150  {
151  return true;
152  }
153 
154  QMap< QgsFeatureId, QgsPointXY >::const_iterator it = pointMap.find( neighborList[0] );
155  if ( it == pointMap.constEnd() ) //should not happen
156  {
157  return true;
158  }
159 
160  QgsPointXY neighborPt = it.value();
161  if ( neighborPt.sqrDist( pt ) < ( minDistance * minDistance ) )
162  {
163  return false;
164  }
165  return true;
166 }
167 
168 
169 
170 
QgsFeatureId id
Definition: qgsfeature.h:70
Wrapper for iterator of features from vector data provider or vector layer.
A rectangle specified with double values.
Definition: qgsrectangle.h:38
int createRandomPoints(QProgressDialog *pd)
Starts calculation of random points.
QList< QgsFeatureId > nearestNeighbor(const QgsPointXY &point, int neighbors) const
Returns nearest neighbors (their count is specified by second parameter)
bool within(const QgsGeometry &geometry) const
Test for if geometry is within another (uses GEOS)
A class to represent a 2D point.
Definition: qgspointxy.h:42
QgsWkbTypes::GeometryType geometryType() const
Returns point, line or polygon.
Container of fields for a vector layer.
Definition: qgsfields.h:41
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:96
bool setAttribute(int field, const QVariant &attr)
Set an attribute&#39;s value by field index.
Definition: qgsfeature.cpp:204
A convenience class for writing vector files to disk.
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:61
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:190
#define MD_RAND_MAX
int mt_rand()
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspointxy.cpp:264
bool addFeature(QgsFeature &feature, QgsFeatureSink::Flags flags=0) override
Adds a single feature to the sink.
QgsFields fields() const override
Returns the list of fields of this layer.
bool isEmpty() const
Returns true if the rectangle is empty.
double width() const
Returns the width of the rectangle.
Definition: qgsrectangle.h:118
This class wraps a request for features to a vector layer (or directly its vector data provider)...
void mt_srand(unsigned value)
QgsCoordinateReferenceSystem crs() const
Returns the layer&#39;s spatial reference system.
bool append(const QgsField &field, FieldOrigin origin=OriginProvider, int originIndex=-1)
Append a field. The field must have unique name, otherwise it is rejected (returns false) ...
Definition: qgsfields.cpp:59
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:46
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:101
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const override
Query the layer for features specified in request.
double yMinimum() const
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:106
bool insertFeature(const QgsFeature &f)
Add feature to index.
QgsRectangle boundingBox() const
Returns the bounding box of the geometry.
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:137
static bool deleteShapeFile(const QString &fileName)
Delete a shapefile (and its accompanying shx / dbf / prf)
double xMinimum() const
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:96
static QgsGeometry fromPoint(const QgsPointXY &point)
Creates a new geometry from a QgsPointXY object.
bool nextFeature(QgsFeature &f)
QgsPointSample(QgsVectorLayer *inputLayer, const QString &outputLayer, const QString &nPointsAttribute, const QString &minDistAttribute=QString())
Represents a vector layer which manages a vector based data sets.
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:255
double height() const
Returns the height of the rectangle.
Definition: qgsrectangle.h:125