Skip to content
Permalink
Browse files

[FEATURE:] Add point sample class to analysis library

  • Loading branch information
mhugent committed Jan 23, 2014
1 parent b086195 commit f351c4f7f8eaab4cc4a93a7a89a62783a45f0fe3
@@ -10,6 +10,7 @@

%Include vector/qgsgeometryanalyzer.sip
%Include vector/qgsoverlayanalyzer.sip
%Include vector/qgspointsample.sip
%Include vector/qgstransectsample.sip
%Include vector/qgszonalstatistics.sip

@@ -0,0 +1,17 @@
/** \ingroup analysis
*/

class QgsPointSample
{
%TypeHeaderCode
#include <qgspointsample.h>
%End

public:
QgsPointSample( QgsVectorLayer* inputLayer, const QString& outputLayer, QString nPointsAttribute, QString minDistAttribute = QString() );
~QgsPointSample();

/**Starts calculation of random points
@return 0 in case of success*/
int createRandomPoints( QProgressDialog* pd );
};
@@ -36,6 +36,7 @@ SET(QGIS_ANALYSIS_SRCS
raster/qgsrastermatrix.cpp
vector/mersenne-twister.cpp
vector/qgsgeometryanalyzer.cpp
vector/qgspointsample.cpp
vector/qgstransectsample.cpp
vector/qgszonalstatistics.cpp
vector/qgsoverlayanalyzer.cpp
@@ -22,7 +22,6 @@

#include <stdio.h>
#include <stdint.h>
#include <limits>
#include "mersenne-twister.h"

/*
@@ -23,6 +23,7 @@
#define MERSENNE_TWISTER_H

#include <stdint.h>
#include <limits>

#ifdef __cplusplus
extern "C" {
@@ -31,7 +32,7 @@ extern "C" {
/*
* Maximum number you can get from rand().
*/
#define RAND_MAX INT32_MAX
#define RAND_MAX std::numeric_limits<int32_t>::max()

/*
* Initialize the number generator with given seed.
@@ -0,0 +1,163 @@
#include "qgspointsample.h"
#include "qgsgeometry.h"
#include "qgsspatialindex.h"
#include "qgsvectorfilewriter.h"
#include "qgsvectorlayer.h"
#include <QFile>
#include "mersenne-twister.h"


QgsPointSample::QgsPointSample( QgsVectorLayer* inputLayer, const QString& outputLayer, QString nPointsAttribute, QString minDistAttribute ): mInputLayer( inputLayer ),
mOutputLayer( outputLayer ), mNumberOfPointsAttribute( nPointsAttribute ), mMinDistanceAttribute( minDistAttribute ), mNCreatedPoints( 0 )
{
}

QgsPointSample::QgsPointSample()
{
}

QgsPointSample::~QgsPointSample()
{
}

int QgsPointSample::createRandomPoints( QProgressDialog* pd )
{
Q_UNUSED( pd );

//create input layer from id (test if polygon, valid)
if ( !mInputLayer )
{
return 1;
}

if ( mInputLayer->geometryType() != QGis::Polygon )
{
return 2;
}

//delete output file if it already exists
if ( QFile::exists( mOutputLayer ) )
{
QgsVectorFileWriter::deleteShapeFile( mOutputLayer );
}

//create vector file writer
QgsFields outputFields;
outputFields.append( QgsField( "id", QVariant::Int ) );
outputFields.append( QgsField( "station_id", QVariant::Int ) );
outputFields.append( QgsField( "stratum_id", QVariant::Int ) );
QgsVectorFileWriter writer( mOutputLayer, "UTF-8",
outputFields,
QGis::WKBPoint,
&( mInputLayer->crs() ) );

//check if creation of output layer successfull
if ( writer.hasError() != QgsVectorFileWriter::NoError )
{
return 3;
}

//init random number generator
mt_srand( QTime::currentTime().msec() );
QgsFeature fet;
int nPoints = 0;
double minDistance = 0;
mNCreatedPoints = 0;

QgsFeatureIterator fIt = mInputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes(
QStringList() << mNumberOfPointsAttribute << mMinDistanceAttribute, mInputLayer->pendingFields() ) );
while ( fIt.nextFeature( fet ) )
{
nPoints = fet.attribute( mNumberOfPointsAttribute ).toInt();
if ( !mMinDistanceAttribute.isEmpty() )
{
minDistance = fet.attribute( mMinDistanceAttribute ).toDouble();
}
addSamplePoints( fet, writer, nPoints, minDistance );
}

return 0;
}

void QgsPointSample::addSamplePoints( QgsFeature& inputFeature, QgsVectorFileWriter& writer, int nPoints, double minDistance )
{
QgsGeometry* geom = inputFeature.geometry();
if ( !geom )
{
return;
}

QgsRectangle geomRect = geom->boundingBox();
if ( geomRect.isEmpty() )
{
return;
}

QgsSpatialIndex sIndex; //to check minimum distance
QMap< QgsFeatureId, QgsPoint > pointMapForFeature;

int nIterations = 0;
int maxIterations = nPoints * 200;
int points = 0;

double randX = 0;
double randY = 0;

while ( nIterations < maxIterations && points < nPoints )
{
randX = (( double )mt_rand() / RAND_MAX ) * geomRect.width() + geomRect.xMinimum();
randY = (( double )mt_rand() / RAND_MAX ) * geomRect.height() + geomRect.yMinimum();
QgsPoint randPoint( randX, randY );
QgsGeometry* ptGeom = QgsGeometry::fromPoint( randPoint );
if ( ptGeom->within( geom ) && checkMinDistance( randPoint, sIndex, minDistance, pointMapForFeature ) )
{
//add feature to writer
QgsFeature f( mNCreatedPoints );
f.setAttribute( "id", mNCreatedPoints + 1 );
f.setAttribute( "station_id", points + 1 );
f.setAttribute( "stratum_id", inputFeature.id() );
f.setGeometry( ptGeom );
writer.addFeature( f );
sIndex.insertFeature( f );
pointMapForFeature.insert( mNCreatedPoints, randPoint );
++points;
++mNCreatedPoints;
}
else
{
delete ptGeom;
}
++nIterations;
}
}

bool QgsPointSample::checkMinDistance( QgsPoint& pt, QgsSpatialIndex& index, double minDistance, QMap< QgsFeatureId, QgsPoint >& pointMap )
{
if ( minDistance <= 0 )
{
return true;
}

QList<QgsFeatureId> neighborList = index.nearestNeighbor( pt, 1 );
if ( neighborList.isEmpty() )
{
return true;
}

QMap< QgsFeatureId, QgsPoint >::const_iterator it = pointMap.find( neighborList[0] );
if ( it == pointMap.constEnd() ) //should not happen
{
return true;
}

QgsPoint neighborPt = it.value();
if ( neighborPt.sqrDist( pt ) < ( minDistance * minDistance ) )
{
return false;
}
return true;
}




@@ -0,0 +1,42 @@
#ifndef QGSPOINTSAMPLE_H
#define QGSPOINTSAMPLE_H

#include "qgsfeature.h"
#include <QString>

class QgsFeature;
class QgsPoint;
class QgsSpatialIndex;
class QgsVectorFileWriter;
class QgsVectorLayer;
class QProgressDialog;

/**Creates random points in polygons / multipolygons*/
class ANALYSIS_EXPORT QgsPointSample
{
public:
QgsPointSample( QgsVectorLayer* inputLayer, const QString& outputLayer, QString nPointsAttribute, QString minDistAttribute = QString() );
~QgsPointSample();

/**Starts calculation of random points
@return 0 in case of success*/
int createRandomPoints( QProgressDialog* pd );

private:

QgsPointSample(); //default constructor is forbidden
void addSamplePoints( QgsFeature& inputFeature, QgsVectorFileWriter& writer, int nPoints, double minDistance );
bool checkMinDistance( QgsPoint& pt, QgsSpatialIndex& index, double minDistance, QMap< QgsFeatureId, QgsPoint >& pointMap );

/**Layer id of input polygon/multipolygon layer*/
QgsVectorLayer* mInputLayer;
/**Output path of result layer*/
QString mOutputLayer;
/**Attribute containing number of points per feature*/
QString mNumberOfPointsAttribute;
/**Attribute containing minimum distance between sample points (or -1 if no min. distance constraint)*/
QString mMinDistanceAttribute;
QgsFeatureId mNCreatedPoints; //helper to find free ids
};

#endif // QGSPOINTSAMPLE_H
@@ -188,7 +188,7 @@ int QgsTransectSample::createSample( QProgressDialog* pd )

while ( nCreatedTransects < nTransects && nIterations < nMaxIterations )
{
double randomPosition = (( double )mt_rand() / std::numeric_limits<int32_t>::max() ) * clippedBaseline->length();
double randomPosition = (( double )mt_rand() / RAND_MAX ) * clippedBaseline->length();
QgsGeometry* samplePoint = clippedBaseline->interpolate( randomPosition );
++nIterations;
if ( !samplePoint )

0 comments on commit f351c4f

Please sign in to comment.
You can’t perform that action at this time.