Skip to content
Permalink
Browse files

Merge pull request #4861 from nyalldawson/nn

Port nearest neighbour alg to new API
  • Loading branch information
nyalldawson committed Jul 15, 2017
2 parents 1cb2728 + 353d4fc commit eaad18c6ad40a0b4a65e5c83853b88ed7f61763c
@@ -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

@@ -83,7 +83,7 @@ def processAlgorithm(self, parameters, context, feedback):
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
@@ -95,7 +95,7 @@ def processAlgorithm(self, parameters, context, feedback):
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
@@ -33,21 +33,25 @@

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 @@ def __init__(self):
super().__init__()

def initAlgorithm(self, config=None):
self.addParameter(ParameterVector(self.POINTS,
self.tr('Points'), [dataobjects.TYPE_VECTOR_POINT]))

self.addOutput(OutputHTML(self.OUTPUT, 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.addParameter(QgsProcessingParameterFeatureSource(self.INPUT,
self.tr('Input layer'), [QgsProcessing.TypeVectorPoint]))

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(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 @@ def displayName(self):
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 @@ def processAlgorithm(self, parameters, context, feedback):
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))

self.createHTML(output, data)

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]))
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

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

return results

def createHTML(self, outputFile, algData):
with codecs.open(outputFile, 'w', encoding='utf-8') as f:
@@ -33,7 +33,7 @@

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 @@ def linearMatrix(self, context, inLayer, inField, targetLayer, targetField,
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 @@ def regularMatrix(self, context, inLayer, inField, targetLayer, targetField,
inIdx = inLayer.fields().lookupField(inField)

distArea = QgsDistanceArea()
distArea.setSourceCrs(inLayer.sourceCrs())
distArea.setEllipsoid(QgsProject.instance().ellipsoid())

first = True
features = QgsProcessingUtils.getFeatures(inLayer, context)
@@ -111,7 +111,7 @@ def processAlgorithm(self, parameters, context, feedback):
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:
@@ -37,6 +37,7 @@
QgsField,
QgsGeometry,
QgsDistanceArea,
QgsProject,
QgsWkbTypes,
QgsProcessingUtils)

@@ -120,6 +121,8 @@ def processAlgorithm(self, parameters, context, feedback):
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
@@ -66,6 +66,7 @@
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 .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 @@ def __init__(self):
self.externalAlgs = []

def getAlgs(self):
# algs = [NearestNeighbourAnalysis(), MeanCoords(),
# algs = [MeanCoords(),
# LinesIntersection(), UniqueValues(), PointDistance(),
# ExportGeometryInfo(),
# SinglePartsToMultiparts(),
@@ -260,6 +260,7 @@ def getAlgs(self):
Intersection(),
LinesToPolygons(),
Merge(),
NearestNeighbourAnalysis(),
PointsInPolygon(),
PointsLayerFromTable(),
PolygonsToLines(),
@@ -41,6 +41,7 @@
QgsFeature,
QgsPointXY,
QgsMessageLog,
QgsProject,
QgsProcessingUtils)

from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
@@ -97,6 +98,9 @@ def processAlgorithm(self, parameters, context, feedback):
points = dict()

da = QgsDistanceArea()
da.setSourceCrs(layer.sourceCrs())
da.setEllipsoid(QgsProject.instance().ellipsoid())

request = QgsFeatureRequest()

random.seed()
@@ -33,7 +33,8 @@
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 @@ def processAlgorithm(self, parameters, context, feedback):
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):
@@ -33,6 +33,7 @@
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 @@ def processAlgorithm(self, parameters, context, feedback):
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):
@@ -83,7 +83,7 @@ def processAlgorithm(self, parameters, context, feedback):
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:
@@ -35,6 +35,7 @@
QgsGeometry,
QgsFeatureRequest,
QgsDistanceArea,
QgsProject,
QgsProcessing,
QgsProcessingParameterString,
QgsProcessingParameterFeatureSource,
@@ -100,9 +101,11 @@ def processAlgorithm(self, parameters, context, feedback):
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
@@ -87,8 +87,8 @@ def processAlgorithm(self, parameters, context, feedback):
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
@@ -97,8 +97,8 @@ def processAlgorithm(self, parameters, context, feedback):
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

0 comments on commit eaad18c

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