sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-04-15 00:04:00 -04:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4861 from nyalldawson/nn

Browse Source 
Port nearest neighbour alg to new API
This commit is contained in:
Nyall Dawson 2017-07-15 17:40:32 +10:00 committed by GitHub
commit eaad18c6ad
17 changed files with 130 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
python/core/qgsspatialindex.sip
View File

@ -28,19 +28,28 @@ class QgsSpatialIndex
Constructor - creates R-tree
%End
explicit QgsSpatialIndex( const QgsFeatureIterator &fi );
explicit QgsSpatialIndex( const QgsFeatureIterator &fi, QgsFeedback *feedback = 0 );
%Docstring
Constructor - creates R-tree and bulk loads it with features from the iterator.
This is much faster approach than creating an empty index and then inserting features one by one.
The optional ``feedback`` object can be used to allow cancelation of bulk feature loading. Ownership
of ``feedback`` is not transferred, and callers must take care that the lifetime of feedback exceeds
that of the spatial index construction.
.. versionadded:: 2.8
%End
explicit QgsSpatialIndex( const QgsFeatureSource &source );
explicit QgsSpatialIndex( const QgsFeatureSource &source, QgsFeedback *feedback = 0 );
%Docstring
Constructor - creates R-tree and bulk loads it with features from the source.
This is much faster approach than creating an empty index and then inserting features one by one.
The optional ``feedback`` object can be used to allow cancelation of bulk feature loading. Ownership
of ``feedback`` is not transferred, and callers must take care that the lifetime of feedback exceeds
that of the spatial index construction.
.. versionadded:: 3.0
%End

2
python/plugins/processing/algs/qgis/Difference.py
View File

@ -83,7 +83,7 @@ class Difference(QgisAlgorithm):
featB = QgsFeature()
outFeat = QgsFeature()
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())))
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())), feedback)
total = 100.0 / (sourceA.featureCount() * sourceB.featureCount()) if sourceA.featureCount() and sourceB.featureCount() else 1
count = 0

2
python/plugins/processing/algs/qgis/Intersection.py
View File

@ -95,7 +95,7 @@ class Intersection(QgisAlgorithm):
fields, geomType, sourceA.sourceCrs())
outFeat = QgsFeature()
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())))
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())), feedback)
total = 100.0 / sourceA.featureCount() if sourceA.featureCount() else 1
count = 0

91
python/plugins/processing/algs/qgis/NearestNeighbourAnalysis.py
View File

@ -33,21 +33,25 @@ import codecs
from qgis.PyQt.QtGui import QIcon
from qgis.core import QgsFeatureRequest, QgsFeature, QgsDistanceArea, QgsProcessingUtils
from qgis.core import (QgsFeatureRequest,
QgsDistanceArea,
QgsProject,
QgsProcessing,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterFileDestination,
QgsProcessingOutputHtml,
QgsProcessingOutputNumber,
QgsSpatialIndex)
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
from processing.core.parameters import ParameterVector
from processing.core.outputs import OutputHTML
from processing.core.outputs import OutputNumber
from processing.tools import dataobjects
pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]
class NearestNeighbourAnalysis(QgisAlgorithm):
POINTS = 'POINTS'
OUTPUT = 'OUTPUT'
INPUT = 'INPUT'
OUTPUT_HTML_FILE = 'OUTPUT_HTML_FILE'
OBSERVED_MD = 'OBSERVED_MD'
EXPECTED_MD = 'EXPECTED_MD'
NN_INDEX = 'NN_INDEX'
@ -64,20 +68,21 @@ class NearestNeighbourAnalysis(QgisAlgorithm):
super().__init__()
def initAlgorithm(self, config=None):
self.addParameter(ParameterVector(self.POINTS,
self.tr('Points'), [dataobjects.TYPE_VECTOR_POINT]))
self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT,
self.tr('Input layer'), [QgsProcessing.TypeVectorPoint]))
self.addOutput(OutputHTML(self.OUTPUT, self.tr('Nearest neighbour')))
self.addParameter(QgsProcessingParameterFileDestination(self.OUTPUT_HTML_FILE, self.tr('Nearest neighbour'), self.tr('HTML files (*.html)'), None, True))
self.addOutput(QgsProcessingOutputHtml(self.OUTPUT_HTML_FILE, self.tr('Nearest neighbour')))
self.addOutput(OutputNumber(self.OBSERVED_MD,
self.tr('Observed mean distance')))
self.addOutput(OutputNumber(self.EXPECTED_MD,
self.tr('Expected mean distance')))
self.addOutput(OutputNumber(self.NN_INDEX,
self.tr('Nearest neighbour index')))
self.addOutput(OutputNumber(self.POINT_COUNT,
self.tr('Number of points')))
self.addOutput(OutputNumber(self.Z_SCORE, self.tr('Z-Score')))
self.addOutput(QgsProcessingOutputNumber(self.OBSERVED_MD,
self.tr('Observed mean distance')))
self.addOutput(QgsProcessingOutputNumber(self.EXPECTED_MD,
self.tr('Expected mean distance')))
self.addOutput(QgsProcessingOutputNumber(self.NN_INDEX,
self.tr('Nearest neighbour index')))
self.addOutput(QgsProcessingOutputNumber(self.POINT_COUNT,
self.tr('Number of points')))
self.addOutput(QgsProcessingOutputNumber(self.Z_SCORE, self.tr('Z-Score')))
def name(self):
return 'nearestneighbouranalysis'
@ -86,26 +91,30 @@ class NearestNeighbourAnalysis(QgisAlgorithm):
return self.tr('Nearest neighbour analysis')
def processAlgorithm(self, parameters, context, feedback):
layer = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.POINTS), context)
output = self.getOutputValue(self.OUTPUT)
source = self.parameterAsSource(parameters, self.INPUT, context)
output_file = self.parameterAsFileOutput(parameters, self.OUTPUT_HTML_FILE, context)
spatialIndex = QgsProcessingUtils.createSpatialIndex(layer, context)
spatialIndex = QgsSpatialIndex(source, feedback)
neighbour = QgsFeature()
distance = QgsDistanceArea()
distance.setSourceCrs(source.sourceCrs())
distance.setEllipsoid(QgsProject.instance().ellipsoid())
sumDist = 0.00
A = layer.extent()
A = source.sourceExtent()
A = float(A.width() * A.height())
features = QgsProcessingUtils.getFeatures(layer, context)
count = QgsProcessingUtils.featureCount(layer, context)
features = source.getFeatures()
count = source.featureCount()
total = 100.0 / count if count else 1
for current, feat in enumerate(features):
if feedback.isCanceled():
break
neighbourID = spatialIndex.nearestNeighbor(
feat.geometry().asPoint(), 2)[1]
request = QgsFeatureRequest().setFilterFid(neighbourID).setSubsetOfAttributes([])
neighbour = next(layer.getFeatures(request))
neighbour = next(source.getFeatures(request))
sumDist += distance.measureLine(neighbour.geometry().asPoint(),
feat.geometry().asPoint())
@ -117,20 +126,24 @@ class NearestNeighbourAnalysis(QgisAlgorithm):
SE = float(0.26136 / math.sqrt(count ** 2 / A))
zscore = float((do - de) / SE)
data = []
data.append('Observed mean distance: ' + str(do))
data.append('Expected mean distance: ' + str(de))
data.append('Nearest neighbour index: ' + str(d))
data.append('Number of points: ' + str(count))
data.append('Z-Score: ' + str(zscore))
results = {}
results[self.OBSERVED_MD] = do
results[self.EXPECTED_MD] = de
results[self.NN_INDEX] = d
results[self.POINT_COUNT] = count
results[self.Z_SCORE] = zscore
self.createHTML(output, data)
if output_file:
data = []
data.append('Observed mean distance: ' + str(do))
data.append('Expected mean distance: ' + str(de))
data.append('Nearest neighbour index: ' + str(d))
data.append('Number of points: ' + str(count))
data.append('Z-Score: ' + str(zscore))
self.createHTML(output_file, data)
results[self.OUTPUT_HTML_FILE] = output_file
self.setOutputValue(self.OBSERVED_MD, float(data[0].split(': ')[1]))
self.setOutputValue(self.EXPECTED_MD, float(data[1].split(': ')[1]))
self.setOutputValue(self.NN_INDEX, float(data[2].split(': ')[1]))
self.setOutputValue(self.POINT_COUNT, float(data[3].split(': ')[1]))
self.setOutputValue(self.Z_SCORE, float(data[4].split(': ')[1]))
return results
def createHTML(self, outputFile, algData):
with codecs.open(outputFile, 'w', encoding='utf-8') as f:

6
python/plugins/processing/algs/qgis/PointDistance.py
View File

@ -33,7 +33,7 @@ import math
from qgis.PyQt.QtGui import QIcon
from qgis.core import QgsFeatureRequest, QgsDistanceArea, QgsFeatureSink, QgsProcessingUtils
from qgis.core import QgsFeatureRequest, QgsProject, QgsDistanceArea, QgsFeatureSink, QgsProcessingUtils
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
from processing.core.parameters import ParameterNumber
@ -134,6 +134,8 @@ class PointDistance(QgisAlgorithm):
outIdx = targetLayer.fields().lookupField(targetField)
distArea = QgsDistanceArea()
distArea.setSourceCrs(inLayer.crs())
distArea.setEllipsoid(QgsProject.instance().ellipsoid())
features = QgsProcessingUtils.getFeatures(inLayer, context)
total = 100.0 / inLayer.featureCount() if inLayer.featureCount() else 0
@ -172,6 +174,8 @@ class PointDistance(QgisAlgorithm):
inIdx = inLayer.fields().lookupField(inField)
distArea = QgsDistanceArea()
distArea.setSourceCrs(inLayer.sourceCrs())
distArea.setEllipsoid(QgsProject.instance().ellipsoid())
first = True
features = QgsProcessingUtils.getFeatures(inLayer, context)

2
python/plugins/processing/algs/qgis/PointsInPolygon.py
View File

@ -111,7 +111,7 @@ class PointsInPolygon(QgisAlgorithm):
fields, poly_source.wkbType(), poly_source.sourceCrs())
spatialIndex = QgsSpatialIndex(point_source.getFeatures(
QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(poly_source.sourceCrs())))
QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(poly_source.sourceCrs())), feedback)
point_attribute_indices = []
if weight_field_index >= 0:

3
python/plugins/processing/algs/qgis/PointsToPaths.py
View File

@ -37,6 +37,7 @@ from qgis.core import (QgsApplication,
QgsField,
QgsGeometry,
QgsDistanceArea,
QgsProject,
QgsWkbTypes,
QgsProcessingUtils)
@ -120,6 +121,8 @@ class PointsToPaths(QgisAlgorithm):
feedback.setProgress(0)
da = QgsDistanceArea()
da.setSourceCrs(layer.sourceCrs())
da.setEllipsoid(QgsProject.instance().ellipsoid())
current = 0
total = 100.0 / len(points) if points else 1

5
python/plugins/processing/algs/qgis/QGISAlgorithmProvider.py
View File

@ -66,6 +66,7 @@ from .ImportIntoSpatialite import ImportIntoSpatialite
from .Intersection import Intersection
from .LinesToPolygons import LinesToPolygons
from .Merge import Merge
from .NearestNeighbourAnalysis import NearestNeighbourAnalysis
from .PointsInPolygon import PointsInPolygon
from .PointsLayerFromTable import PointsLayerFromTable
from .PolygonsToLines import PolygonsToLines
@ -90,7 +91,6 @@ from .VoronoiPolygons import VoronoiPolygons
from .ZonalStatistics import ZonalStatistics
# from .ExtractByLocation import ExtractByLocation
# from .NearestNeighbourAnalysis import NearestNeighbourAnalysis
# from .LinesIntersection import LinesIntersection
# from .MeanCoords import MeanCoords
# from .PointDistance import PointDistance
@ -183,7 +183,7 @@ class QGISAlgorithmProvider(QgsProcessingProvider):
self.externalAlgs = []
def getAlgs(self):
# algs = [NearestNeighbourAnalysis(), MeanCoords(),
# algs = [MeanCoords(),
# LinesIntersection(), UniqueValues(), PointDistance(),
# ExportGeometryInfo(),
# SinglePartsToMultiparts(),
@ -260,6 +260,7 @@ class QGISAlgorithmProvider(QgsProcessingProvider):
Intersection(),
LinesToPolygons(),
Merge(),
NearestNeighbourAnalysis(),
PointsInPolygon(),
PointsLayerFromTable(),
PolygonsToLines(),

4
python/plugins/processing/algs/qgis/RandomPointsAlongLines.py
View File

@ -41,6 +41,7 @@ from qgis.core import (QgsApplication,
QgsFeature,
QgsPointXY,
QgsMessageLog,
QgsProject,
QgsProcessingUtils)
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
@ -97,6 +98,9 @@ class RandomPointsAlongLines(QgisAlgorithm):
points = dict()
da = QgsDistanceArea()
da.setSourceCrs(layer.sourceCrs())
da.setEllipsoid(QgsProject.instance().ellipsoid())
request = QgsFeatureRequest()
random.seed()

5
python/plugins/processing/algs/qgis/RandomPointsPolygonsFixed.py
View File

@ -33,7 +33,8 @@ from qgis.PyQt.QtCore import QVariant
from qgis.core import (QgsFields, QgsFeatureSink, QgsField, QgsDistanceArea, QgsGeometry, QgsWkbTypes,
QgsSpatialIndex, QgsPointXY, QgsFeature,
QgsMessageLog,
QgsProcessingUtils)
QgsProcessingUtils,
QgsProject)
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
from processing.core.parameters import ParameterVector
@ -94,6 +95,8 @@ class RandomPointsPolygonsFixed(QgisAlgorithm):
writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(fields, QgsWkbTypes.Point, layer.crs(), context)
da = QgsDistanceArea()
da.setSourceCrs(layer.sourceCrs())
da.setEllipsoid(QgsProject.instance().ellipsoid())
features = QgsProcessingUtils.getFeatures(layer, context)
for current, f in enumerate(features):

3
python/plugins/processing/algs/qgis/RandomPointsPolygonsVariable.py
View File

@ -33,6 +33,7 @@ from qgis.PyQt.QtCore import QVariant
from qgis.core import (QgsFields, QgsFeatureSink, QgsField, QgsFeature, QgsPointXY, QgsWkbTypes,
QgsGeometry, QgsSpatialIndex, QgsDistanceArea,
QgsMessageLog,
QgsProject,
QgsProcessingUtils)
from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
@ -95,6 +96,8 @@ class RandomPointsPolygonsVariable(QgisAlgorithm):
writer = self.getOutputFromName(self.OUTPUT).getVectorWriter(fields, QgsWkbTypes.Point, layer.crs(), context)
da = QgsDistanceArea()
da.setSourceCrs(layer.sourceCrs())
da.setEllipsoid(QgsProject.instance().ellipsoid())
features = QgsProcessingUtils.getFeatures(layer, context)
for current, f in enumerate(features):

2
python/plugins/processing/algs/qgis/SelectByAttributeSum.py
View File

@ -83,7 +83,7 @@ class SelectByAttributeSum(QgisAlgorithm):
geom = ft.geometry()
attrSum = ft[fieldName]
idx = QgsSpatialIndex(layer.getFeatures(QgsFeatureRequest.setSubsetOfAttributes([])))
idx = QgsSpatialIndex(layer.getFeatures(QgsFeatureRequest.setSubsetOfAttributes([])), feedback)
req = QgsFeatureRequest()
completed = False
while not completed:

5
python/plugins/processing/algs/qgis/SumLines.py
View File

@ -35,6 +35,7 @@ from qgis.core import (QgsFeature,
QgsGeometry,
QgsFeatureRequest,
QgsDistanceArea,
QgsProject,
QgsProcessing,
QgsProcessingParameterString,
QgsProcessingParameterFeatureSource,
@ -100,9 +101,11 @@ class SumLines(QgisAlgorithm):
fields, poly_source.wkbType(), poly_source.sourceCrs())
spatialIndex = QgsSpatialIndex(line_source.getFeatures(
QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(poly_source.sourceCrs())))
QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(poly_source.sourceCrs())), feedback)
distArea = QgsDistanceArea()
distArea.setSourceCrs(poly_source.sourceCrs())
distArea.setEllipsoid(QgsProject.instance().ellipsoid())
features = poly_source.getFeatures()
total = 100.0 / poly_source.featureCount() if poly_source.featureCount() else 0

4
python/plugins/processing/algs/qgis/SymmetricalDifference.py
View File

@ -87,8 +87,8 @@ class SymmetricalDifference(QgisAlgorithm):
featB = QgsFeature()
outFeat = QgsFeature()
indexA = QgsSpatialIndex(sourceA)
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())))
indexA = QgsSpatialIndex(sourceA, feedback)
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())), feedback)
total = 100.0 / (sourceA.featureCount() * sourceB.featureCount()) if sourceA.featureCount() and sourceB.featureCount() else 1
count = 0

4
python/plugins/processing/algs/qgis/Union.py
View File

@ -97,8 +97,8 @@ class Union(QgisAlgorithm):
featB = QgsFeature()
outFeat = QgsFeature()
indexA = QgsSpatialIndex(sourceA)
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())))
indexA = QgsSpatialIndex(sourceA, feedback)
indexB = QgsSpatialIndex(sourceB.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(sourceA.sourceCrs())), feedback)
total = 100.0 / (sourceA.featureCount() * sourceB.featureCount()) if sourceA.featureCount() and sourceB.featureCount() else 1
count = 0

28
src/core/qgsspatialindex.cpp
View File

@ -21,6 +21,7 @@
#include "qgsrectangle.h"
#include "qgslogger.h"
#include "qgsfeaturesource.h"
#include "qgsfeedback.h"
#include "SpatialIndex.h"
@ -92,9 +93,10 @@ class QgsFeatureIteratorDataStream : public IDataStream
{
public:
//! constructor - needs to load all data to a vector for later access when bulk loading
explicit QgsFeatureIteratorDataStream( const QgsFeatureIterator &fi )
explicit QgsFeatureIteratorDataStream( const QgsFeatureIterator &fi, QgsFeedback *feedback = nullptr )
: mFi( fi )
, mNextData( nullptr )
, mFeedback( feedback )
{
readNextEntry();
}
@ -107,6 +109,9 @@ class QgsFeatureIteratorDataStream : public IDataStream
//! returns a pointer to the next entry in the stream or 0 at the end of the stream.
IData *getNext() override
{
if ( mFeedback && mFeedback->isCanceled() )
return nullptr;
RTree::Data *ret = mNextData;
mNextData = nullptr;
readNextEntry();
@ -141,6 +146,7 @@ class QgsFeatureIteratorDataStream : public IDataStream
private:
QgsFeatureIterator mFi;
RTree::Data *mNextData = nullptr;
QgsFeedback *mFeedback = nullptr;
};
@ -157,9 +163,17 @@ class QgsSpatialIndexData : public QSharedData
initTree();
}
explicit QgsSpatialIndexData( const QgsFeatureIterator &fi )
/**
* Constructor for QgsSpatialIndexData which bulk loads features from the specified feature iterator
* \a fi.
*
* The optional \a feedback object can be used to allow cancelation of bulk feature loading. Ownership
* of \a feedback is not transferred, and callers must take care that the lifetime of feedback exceeds
* that of the spatial index construction.
*/
explicit QgsSpatialIndexData( const QgsFeatureIterator &fi, QgsFeedback *feedback = nullptr )
{
QgsFeatureIteratorDataStream fids( fi );
QgsFeatureIteratorDataStream fids( fi, feedback );
initTree( &fids );
}
@ -224,14 +238,14 @@ QgsSpatialIndex::QgsSpatialIndex()
d = new QgsSpatialIndexData;
}
QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureIterator &fi )
QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureIterator &fi, QgsFeedback *feedback )
{
d = new QgsSpatialIndexData( fi );
d = new QgsSpatialIndexData( fi, feedback );
}
QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureSource &source )
QgsSpatialIndex::QgsSpatialIndex( const QgsFeatureSource &source, QgsFeedback *feedback )
{
d = new QgsSpatialIndexData( source.getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QgsAttributeList() ) ) );
d = new QgsSpatialIndexData( source.getFeatures( QgsFeatureRequest().setSubsetOfAttributes( QgsAttributeList() ) ), feedback );
}
QgsSpatialIndex::QgsSpatialIndex( const QgsSpatialIndex &other ) //NOLINT

14
src/core/qgsspatialindex.h
View File

@ -33,6 +33,7 @@ namespace SpatialIndex SIP_SKIP
}
}
class QgsFeedback;
class QgsFeature;
class QgsRectangle;
class QgsPointXY;
@ -64,17 +65,26 @@ class CORE_EXPORT QgsSpatialIndex
/** Constructor - creates R-tree and bulk loads it with features from the iterator.
* This is much faster approach than creating an empty index and then inserting features one by one.
*
* The optional \a feedback object can be used to allow cancelation of bulk feature loading. Ownership
* of \a feedback is not transferred, and callers must take care that the lifetime of feedback exceeds
* that of the spatial index construction.
*
* \since QGIS 2.8
*/
explicit QgsSpatialIndex( const QgsFeatureIterator &fi );
explicit QgsSpatialIndex( const QgsFeatureIterator &fi, QgsFeedback *feedback = nullptr );
/**
* Constructor - creates R-tree and bulk loads it with features from the source.
* This is much faster approach than creating an empty index and then inserting features one by one.
*
* The optional \a feedback object can be used to allow cancelation of bulk feature loading. Ownership
* of \a feedback is not transferred, and callers must take care that the lifetime of feedback exceeds
* that of the spatial index construction.
*
* \since QGIS 3.0
*/
explicit QgsSpatialIndex( const QgsFeatureSource &source );
explicit QgsSpatialIndex( const QgsFeatureSource &source, QgsFeedback *feedback = nullptr );
//! Copy constructor
QgsSpatialIndex( const QgsSpatialIndex &other );