QGIS API Documentation  2.15.0-Master (972fc9f)
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 "qgsgeometry.h"
17 #include "qgsspatialindex.h"
18 #include "qgsvectorfilewriter.h"
19 #include "qgsvectorlayer.h"
20 #include <QFile>
21 #include "mersenne-twister.h"
22 
23 
24 QgsPointSample::QgsPointSample( QgsVectorLayer* inputLayer, const QString& outputLayer, const QString& nPointsAttribute, const QString& minDistAttribute ): mInputLayer( inputLayer ),
25  mOutputLayer( outputLayer ), mNumberOfPointsAttribute( nPointsAttribute ), mMinDistanceAttribute( minDistAttribute ), mNCreatedPoints( 0 )
26 {
27 }
28 
30  : mInputLayer( nullptr )
31  , mNCreatedPoints( 0 )
32 {
33 }
34 
36 {
37  Q_UNUSED( pd );
38 
39  //create input layer from id (test if polygon, valid)
40  if ( !mInputLayer )
41  {
42  return 1;
43  }
44 
45  if ( mInputLayer->geometryType() != QGis::Polygon )
46  {
47  return 2;
48  }
49 
50  //delete output file if it already exists
51  if ( QFile::exists( mOutputLayer ) )
52  {
54  }
55 
56  //create vector file writer
57  QgsFields outputFields;
58  outputFields.append( QgsField( "id", QVariant::Int ) );
59  outputFields.append( QgsField( "station_id", QVariant::Int ) );
60  outputFields.append( QgsField( "stratum_id", QVariant::Int ) );
61  QgsVectorFileWriter writer( mOutputLayer, "UTF-8",
62  outputFields,
64  &( mInputLayer->crs() ) );
65 
66  //check if creation of output layer successfull
67  if ( writer.hasError() != QgsVectorFileWriter::NoError )
68  {
69  return 3;
70  }
71 
72  //init random number generator
73  mt_srand( QTime::currentTime().msec() );
74  QgsFeature fet;
75  int nPoints = 0;
76  double minDistance = 0;
77  mNCreatedPoints = 0;
78 
79  QgsFeatureIterator fIt = mInputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes(
80  QStringList() << mNumberOfPointsAttribute << mMinDistanceAttribute, mInputLayer->fields() ) );
81  while ( fIt.nextFeature( fet ) )
82  {
83  nPoints = fet.attribute( mNumberOfPointsAttribute ).toInt();
84  if ( !mMinDistanceAttribute.isEmpty() )
85  {
86  minDistance = fet.attribute( mMinDistanceAttribute ).toDouble();
87  }
88  addSamplePoints( fet, writer, nPoints, minDistance );
89  }
90 
91  return 0;
92 }
93 
94 void QgsPointSample::addSamplePoints( QgsFeature& inputFeature, QgsVectorFileWriter& writer, int nPoints, double minDistance )
95 {
96  if ( !inputFeature.constGeometry() )
97  return;
98 
99  const QgsGeometry* geom = inputFeature.constGeometry();
100  QgsRectangle geomRect = geom->boundingBox();
101  if ( geomRect.isEmpty() )
102  {
103  return;
104  }
105 
106  QgsSpatialIndex sIndex; //to check minimum distance
107  QMap< QgsFeatureId, QgsPoint > pointMapForFeature;
108 
109  int nIterations = 0;
110  int maxIterations = nPoints * 200;
111  int points = 0;
112 
113  double randX = 0;
114  double randY = 0;
115 
116  while ( nIterations < maxIterations && points < nPoints )
117  {
118  randX = (( double )mt_rand() / MD_RAND_MAX ) * geomRect.width() + geomRect.xMinimum();
119  randY = (( double )mt_rand() / MD_RAND_MAX ) * geomRect.height() + geomRect.yMinimum();
120  QgsPoint randPoint( randX, randY );
121  QgsGeometry* ptGeom = QgsGeometry::fromPoint( randPoint );
122  if ( ptGeom->within( geom ) && checkMinDistance( randPoint, sIndex, minDistance, pointMapForFeature ) )
123  {
124  //add feature to writer
125  QgsFeature f( mNCreatedPoints );
126  f.setAttribute( "id", mNCreatedPoints + 1 );
127  f.setAttribute( "station_id", points + 1 );
128  f.setAttribute( "stratum_id", inputFeature.id() );
129  f.setGeometry( ptGeom );
130  writer.addFeature( f );
131  sIndex.insertFeature( f );
132  pointMapForFeature.insert( mNCreatedPoints, randPoint );
133  ++points;
134  ++mNCreatedPoints;
135  }
136  else
137  {
138  delete ptGeom;
139  }
140  ++nIterations;
141  }
142 }
143 
144 bool QgsPointSample::checkMinDistance( QgsPoint& pt, QgsSpatialIndex& index, double minDistance, QMap< QgsFeatureId, QgsPoint >& pointMap )
145 {
146  if ( minDistance <= 0 )
147  {
148  return true;
149  }
150 
151  QList<QgsFeatureId> neighborList = index.nearestNeighbor( pt, 1 );
152  if ( neighborList.isEmpty() )
153  {
154  return true;
155  }
156 
157  QMap< QgsFeatureId, QgsPoint >::const_iterator it = pointMap.find( neighborList[0] );
158  if ( it == pointMap.constEnd() ) //should not happen
159  {
160  return true;
161  }
162 
163  QgsPoint neighborPt = it.value();
164  if ( neighborPt.sqrDist( pt ) < ( minDistance * minDistance ) )
165  {
166  return false;
167  }
168  return true;
169 }
170 
171 
172 
173 
QgsFeatureId id() const
Get the feature ID for this feature.
Definition: qgsfeature.cpp:65
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
bool isEmpty() const
test if rectangle is empty.
int createRandomPoints(QProgressDialog *pd)
Starts calculation of random points.
static bool deleteShapeFile(const QString &theFileName)
Delete a shapefile (and its accompanying shx / dbf / prf)
QgsFields fields() const
Returns the list of fields of this layer.
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest())
Query the provider for features specified in request.
QgsRectangle boundingBox() const
Returns the bounding box of this feature.
Container of fields for a vector layer.
Definition: qgsfield.h:193
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:76
bool setAttribute(int field, const QVariant &attr)
Set an attribute&#39;s value by field index.
Definition: qgsfeature.cpp:222
bool addFeature(QgsFeature &feature, QgsFeatureRendererV2 *renderer=nullptr, QGis::UnitType outputUnit=QGis::Meters)
Add feature to the currently opened data source.
bool exists() const
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:187
double sqrDist(double x, double y) const
Returns the squared distance between this point a specified x, y coordinate.
Definition: qgspoint.cpp:353
#define MD_RAND_MAX
int mt_rand()
void setGeometry(const QgsGeometry &geom)
Set this feature&#39;s geometry from another QgsGeometry object.
Definition: qgsfeature.cpp:124
int toInt(bool *ok) const
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:202
bool isEmpty() const
bool isEmpty() const
const_iterator constEnd() const
This class wraps a request for features to a vector layer (or directly its vector data provider)...
void mt_srand(unsigned value)
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: qgsfield.cpp:309
const T & value() const
QGis::GeometryType geometryType() const
Returns point, line or polygon.
Encapsulate a field in an attribute table or data source.
Definition: qgsfield.h:44
A class to represent a point.
Definition: qgspoint.h:117
QList< QgsFeatureId > nearestNeighbor(const QgsPoint &point, int neighbors) const
Returns nearest neighbors (their count is specified by second parameter)
static QgsGeometry * fromPoint(const QgsPoint &point)
Creates a new geometry from a QgsPoint object.
QTime currentTime()
bool insertFeature(const QgsFeature &f)
Add feature to index.
WriterError hasError()
Checks whether there were any errors in constructor.
QVariant attribute(const QString &name) const
Lookup attribute value from attribute name.
Definition: qgsfeature.cpp:271
bool within(const QgsGeometry *geometry) const
Test for if geometry is within another (uses GEOS)
const QgsGeometry * constGeometry() const
Gets a const pointer to the geometry object associated with this feature.
Definition: qgsfeature.cpp:82
const QgsCoordinateReferenceSystem & crs() const
Returns layer&#39;s spatial reference system.
double toDouble(bool *ok) const
iterator insert(const Key &key, const T &value)
bool nextFeature(QgsFeature &f)
QgsPointSample(QgsVectorLayer *inputLayer, const QString &outputLayer, const QString &nPointsAttribute, const QString &minDistAttribute=QString())
double width() const
Width of the rectangle.
Definition: qgsrectangle.h:207
Represents a vector layer which manages a vector based data sets.
double xMinimum() const
Get the x minimum value (left side of rectangle)
Definition: qgsrectangle.h:192
iterator find(const Key &key)
double height() const
Height of the rectangle.
Definition: qgsrectangle.h:212