QGIS API Documentation  2.99.0-Master (37c43df)
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 successfull
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, QgsPoint > 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  QgsPoint 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( QgsPoint& pt, QgsSpatialIndex& index, double minDistance, QMap< QgsFeatureId, QgsPoint >& 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, QgsPoint >::const_iterator it = pointMap.find( neighborList[0] );
155  if ( it == pointMap.constEnd() ) //should not happen
156  {
157  return true;
158  }
159 
160  QgsPoint 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:139
Wrapper for iterator of features from vector data provider or vector layer.
static unsigned index
A rectangle specified with double values.
Definition: qgsrectangle.h:35
int createRandomPoints(QProgressDialog *pd)
Starts calculation of random points.
static bool deleteShapeFile(const QString &theFileName)
Delete a shapefile (and its accompanying shx / dbf / prf)
bool within(const QgsGeometry &geometry) const
Test for if geometry is within another (uses GEOS)
bool addFeature(QgsFeature &feature, QgsFeatureRenderer *renderer=nullptr, QgsUnitTypes::DistanceUnit outputUnit=QgsUnitTypes::DistanceMeters)
Add feature to the currently opened data source.
QList< QgsFeatureId > nearestNeighbor(const QgsPoint &point, int neighbors) const
Returns nearest neighbors (their count is specified by second parameter)
QgsWkbTypes::GeometryType geometryType() const
Returns point, line or polygon.
Container of fields for a vector layer.
Definition: qgsfields.h:36
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:78
bool setAttribute(int field, const QVariant &attr)
Set an attribute&#39;s value by field index.
Definition: qgsfeature.cpp:228
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:135
bool hasGeometry() const
Returns true if the feature has an associated geometry.
Definition: qgsfeature.cpp:214
#define MD_RAND_MAX
int mt_rand()
QgsFields fields() const
Returns the list of fields of this layer.
bool isEmpty() const
test if rectangle is empty.
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:211
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest()) const
Query the layer for features specified in request.
static QgsGeometry fromPoint(const QgsPoint &point)
Creates a new geometry from a QgsPoint object.
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:61
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:47
QgsGeometry geometry() const
Returns the geometry associated with this feature.
Definition: qgsfeature.cpp:113
A class to represent a point.
Definition: qgspoint.h:111
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:206
bool insertFeature(const QgsFeature &f)
Add feature to index.
QgsRectangle boundingBox() const
Returns the bounding box of this feature.
void setGeometry(const QgsGeometry &geometry)
Set the feature&#39;s geometry.
Definition: qgsfeature.cpp:162
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:196
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspoint.cpp:348
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:277
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:216