QGIS API Documentation  2.5.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
qgspointsample.cpp
Go to the documentation of this file.
1 #include "qgspointsample.h"
2 #include "qgsgeometry.h"
3 #include "qgsspatialindex.h"
4 #include "qgsvectorfilewriter.h"
5 #include "qgsvectorlayer.h"
6 #include <QFile>
7 #include "mersenne-twister.h"
8 
9 
10 QgsPointSample::QgsPointSample( QgsVectorLayer* inputLayer, const QString& outputLayer, QString nPointsAttribute, QString minDistAttribute ): mInputLayer( inputLayer ),
11  mOutputLayer( outputLayer ), mNumberOfPointsAttribute( nPointsAttribute ), mMinDistanceAttribute( minDistAttribute ), mNCreatedPoints( 0 )
12 {
13 }
14 
16 {
17 }
18 
20 {
21 }
22 
23 int QgsPointSample::createRandomPoints( QProgressDialog* pd )
24 {
25  Q_UNUSED( pd );
26 
27  //create input layer from id (test if polygon, valid)
28  if ( !mInputLayer )
29  {
30  return 1;
31  }
32 
34  {
35  return 2;
36  }
37 
38  //delete output file if it already exists
39  if ( QFile::exists( mOutputLayer ) )
40  {
42  }
43 
44  //create vector file writer
45  QgsFields outputFields;
46  outputFields.append( QgsField( "id", QVariant::Int ) );
47  outputFields.append( QgsField( "station_id", QVariant::Int ) );
48  outputFields.append( QgsField( "stratum_id", QVariant::Int ) );
49  QgsVectorFileWriter writer( mOutputLayer, "UTF-8",
50  outputFields,
52  &( mInputLayer->crs() ) );
53 
54  //check if creation of output layer successfull
55  if ( writer.hasError() != QgsVectorFileWriter::NoError )
56  {
57  return 3;
58  }
59 
60  //init random number generator
61  mt_srand( QTime::currentTime().msec() );
62  QgsFeature fet;
63  int nPoints = 0;
64  double minDistance = 0;
65  mNCreatedPoints = 0;
66 
67  QgsFeatureIterator fIt = mInputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes(
69  while ( fIt.nextFeature( fet ) )
70  {
71  nPoints = fet.attribute( mNumberOfPointsAttribute ).toInt();
72  if ( !mMinDistanceAttribute.isEmpty() )
73  {
74  minDistance = fet.attribute( mMinDistanceAttribute ).toDouble();
75  }
76  addSamplePoints( fet, writer, nPoints, minDistance );
77  }
78 
79  return 0;
80 }
81 
82 void QgsPointSample::addSamplePoints( QgsFeature& inputFeature, QgsVectorFileWriter& writer, int nPoints, double minDistance )
83 {
84  QgsGeometry* geom = inputFeature.geometry();
85  if ( !geom )
86  {
87  return;
88  }
89 
90  QgsRectangle geomRect = geom->boundingBox();
91  if ( geomRect.isEmpty() )
92  {
93  return;
94  }
95 
96  QgsSpatialIndex sIndex; //to check minimum distance
97  QMap< QgsFeatureId, QgsPoint > pointMapForFeature;
98 
99  int nIterations = 0;
100  int maxIterations = nPoints * 200;
101  int points = 0;
102 
103  double randX = 0;
104  double randY = 0;
105 
106  while ( nIterations < maxIterations && points < nPoints )
107  {
108  randX = (( double )mt_rand() / MD_RAND_MAX ) * geomRect.width() + geomRect.xMinimum();
109  randY = (( double )mt_rand() / MD_RAND_MAX ) * geomRect.height() + geomRect.yMinimum();
110  QgsPoint randPoint( randX, randY );
111  QgsGeometry* ptGeom = QgsGeometry::fromPoint( randPoint );
112  if ( ptGeom->within( geom ) && checkMinDistance( randPoint, sIndex, minDistance, pointMapForFeature ) )
113  {
114  //add feature to writer
116  f.setAttribute( "id", mNCreatedPoints + 1 );
117  f.setAttribute( "station_id", points + 1 );
118  f.setAttribute( "stratum_id", inputFeature.id() );
119  f.setGeometry( ptGeom );
120  writer.addFeature( f );
121  sIndex.insertFeature( f );
122  pointMapForFeature.insert( mNCreatedPoints, randPoint );
123  ++points;
124  ++mNCreatedPoints;
125  }
126  else
127  {
128  delete ptGeom;
129  }
130  ++nIterations;
131  }
132 }
133 
134 bool QgsPointSample::checkMinDistance( QgsPoint& pt, QgsSpatialIndex& index, double minDistance, QMap< QgsFeatureId, QgsPoint >& pointMap )
135 {
136  if ( minDistance <= 0 )
137  {
138  return true;
139  }
140 
141  QList<QgsFeatureId> neighborList = index.nearestNeighbor( pt, 1 );
142  if ( neighborList.isEmpty() )
143  {
144  return true;
145  }
146 
147  QMap< QgsFeatureId, QgsPoint >::const_iterator it = pointMap.find( neighborList[0] );
148  if ( it == pointMap.constEnd() ) //should not happen
149  {
150  return true;
151  }
152 
153  QgsPoint neighborPt = it.value();
154  if ( neighborPt.sqrDist( pt ) < ( minDistance * minDistance ) )
155  {
156  return false;
157  }
158  return true;
159 }
160 
161 
162 
163