Skip to content
Permalink
Browse files

Port Distance Matrix algorithm to new API

Enhancements:
- support source/target layers in different CRS
- output layers with geometry (i.e. keep input point geometry - avoids
need to rejoin result back to original table to get geometry)
- keep original data types for id fields
- don't fire off many single feature requests - instead request
multiple features at once to improve speed
  • Loading branch information
nyalldawson committed Jul 15, 2017
1 parent 7f58af1 commit b7f888bf5bf68412a4227996db3d231464265474
@@ -32,29 +32,38 @@
import math

from qgis.PyQt.QtGui import QIcon

from qgis.core import QgsFeatureRequest, QgsProject, QgsDistanceArea, QgsFeatureSink, QgsProcessingUtils
from qgis.PyQt.QtCore import QVariant

from qgis.core import (QgsFeatureRequest,
QgsField,
QgsFields,
QgsProject,
QgsFeature,
QgsGeometry,
QgsDistanceArea,
QgsFeatureSink,
QgsProcessingParameterFeatureSource,
QgsProcessing,
QgsProcessingParameterEnum,
QgsProcessingParameterField,
QgsProcessingParameterNumber,
QgsProcessingParameterFeatureSink,
QgsSpatialIndex,
QgsWkbTypes)

from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
from processing.core.parameters import ParameterNumber
from processing.core.parameters import ParameterVector
from processing.core.parameters import ParameterSelection
from processing.core.parameters import ParameterTableField
from processing.core.outputs import OutputTable
from processing.tools import dataobjects

pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]


class PointDistance(QgisAlgorithm):

INPUT_LAYER = 'INPUT_LAYER'
INPUT = 'INPUT'
INPUT_FIELD = 'INPUT_FIELD'
TARGET_LAYER = 'TARGET_LAYER'
TARGET = 'TARGET'
TARGET_FIELD = 'TARGET_FIELD'
MATRIX_TYPE = 'MATRIX_TYPE'
NEAREST_POINTS = 'NEAREST_POINTS'
DISTANCE_MATRIX = 'DISTANCE_MATRIX'
OUTPUT = 'OUTPUT'

def icon(self):
return QIcon(os.path.join(pluginPath, 'images', 'ftools', 'matrix.png'))
@@ -70,22 +79,26 @@ def initAlgorithm(self, config=None):
self.tr('Standard (N x T) distance matrix'),
self.tr('Summary distance matrix (mean, std. dev., min, max)')]

self.addParameter(ParameterVector(self.INPUT_LAYER,
self.tr('Input point layer'), [dataobjects.TYPE_VECTOR_POINT]))
self.addParameter(ParameterTableField(self.INPUT_FIELD,
self.tr('Input unique ID field'), self.INPUT_LAYER,
ParameterTableField.DATA_TYPE_ANY))
self.addParameter(ParameterVector(self.TARGET_LAYER,
self.tr('Target point layer'), dataobjects.TYPE_VECTOR_POINT))
self.addParameter(ParameterTableField(self.TARGET_FIELD,
self.tr('Target unique ID field'), self.TARGET_LAYER,
ParameterTableField.DATA_TYPE_ANY))
self.addParameter(ParameterSelection(self.MATRIX_TYPE,
self.tr('Output matrix type'), self.mat_types, 0))
self.addParameter(ParameterNumber(self.NEAREST_POINTS,
self.tr('Use only the nearest (k) target points'), 0, 9999, 0))

self.addOutput(OutputTable(self.DISTANCE_MATRIX, self.tr('Distance matrix')))
self.addParameter(QgsProcessingParameterFeatureSource(self.INPUT,
self.tr('Input point layer'),
[QgsProcessing.TypeVectorPoint]))
self.addParameter(QgsProcessingParameterField(self.INPUT_FIELD,
self.tr('Input unique ID field'),
parentLayerParameterName=self.INPUT,
type=QgsProcessingParameterField.Any))
self.addParameter(QgsProcessingParameterFeatureSource(self.TARGET,
self.tr('Target point layer'),
[QgsProcessing.TypeVectorPoint]))
self.addParameter(QgsProcessingParameterField(self.TARGET_FIELD,
self.tr('Target unique ID field'),
parentLayerParameterName=self.TARGET,
type=QgsProcessingParameterField.Any))
self.addParameter(QgsProcessingParameterEnum(self.MATRIX_TYPE,
self.tr('Output matrix type'), options=self.mat_types, defaultValue=0))
self.addParameter(QgsProcessingParameterNumber(self.NEAREST_POINTS,
self.tr('Use only the nearest (k) target points'), type=QgsProcessingParameterNumber.Integer, minValue=0, maxValue=9999, defaultValue=0))

self.addParameter(QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr('Distance matrix'), QgsProcessing.TypeVectorPoint))

def name(self):
return 'distancematrix'
@@ -94,65 +107,86 @@ def displayName(self):
return self.tr('Distance matrix')

def processAlgorithm(self, parameters, context, feedback):
inLayer = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.INPUT_LAYER), context)
inField = self.getParameterValue(self.INPUT_FIELD)
targetLayer = QgsProcessingUtils.mapLayerFromString(self.getParameterValue(self.TARGET_LAYER), context)
targetField = self.getParameterValue(self.TARGET_FIELD)
matType = self.getParameterValue(self.MATRIX_TYPE)
nPoints = self.getParameterValue(self.NEAREST_POINTS)

outputFile = self.getOutputFromName(self.DISTANCE_MATRIX)
source = self.parameterAsSource(parameters, self.INPUT, context)
source_field = self.parameterAsString(parameters, self.INPUT_FIELD, context)
target_source = self.parameterAsSource(parameters, self.TARGET, context)
target_field = self.parameterAsString(parameters, self.TARGET_FIELD, context)
matType = self.parameterAsEnum(parameters, self.MATRIX_TYPE, context)
nPoints = self.parameterAsInt(parameters, self.NEAREST_POINTS, context)

if nPoints < 1:
nPoints = QgsProcessingUtils.featureCount(targetLayer, context)

self.writer = outputFile.getTableWriter([])
nPoints = target_source.featureCount()

if matType == 0:
# Linear distance matrix
self.linearMatrix(context, inLayer, inField, targetLayer, targetField,
matType, nPoints, feedback)
return self.linearMatrix(parameters, context, source, source_field, target_source, target_field,
matType, nPoints, feedback)
elif matType == 1:
# Standard distance matrix
self.regularMatrix(context, inLayer, inField, targetLayer, targetField,
nPoints, feedback)
return self.regularMatrix(parameters, context, source, source_field, target_source, target_field,
nPoints, feedback)
elif matType == 2:
# Summary distance matrix
self.linearMatrix(context, inLayer, inField, targetLayer, targetField,
matType, nPoints, feedback)
return self.linearMatrix(parameters, context, source, source_field, target_source, target_field,
matType, nPoints, feedback)

def linearMatrix(self, context, inLayer, inField, targetLayer, targetField,
def linearMatrix(self, parameters, context, source, inField, target_source, targetField,
matType, nPoints, feedback):
inIdx = source.fields().lookupField(inField)
outIdx = target_source.fields().lookupField(targetField)

fields = QgsFields()
input_id_field = source.fields()[inIdx]
input_id_field.setName('InputID')
fields.append(input_id_field)
if matType == 0:
self.writer.addRecord(['InputID', 'TargetID', 'Distance'])
target_id_field = target_source.fields()[outIdx]
target_id_field.setName('TargetID')
fields.append(target_id_field)
fields.append(QgsField('Distance', QVariant.Double))
else:
self.writer.addRecord(['InputID', 'MEAN', 'STDDEV', 'MIN', 'MAX'])
fields.append(QgsField('MEAN', QVariant.Double))
fields.append(QgsField('STDDEV', QVariant.Double))
fields.append(QgsField('MIN', QVariant.Double))
fields.append(QgsField('MAX', QVariant.Double))

index = QgsProcessingUtils.createSpatialIndex(targetLayer, context)
out_wkb = QgsWkbTypes.multiType(source.wkbType()) if matType == 0 else source.wkbType()
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
fields, out_wkb, source.sourceCrs())

inIdx = inLayer.fields().lookupField(inField)
outIdx = targetLayer.fields().lookupField(targetField)
index = QgsSpatialIndex(target_source.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(source.sourceCrs())), feedback)

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

features = QgsProcessingUtils.getFeatures(inLayer, context)
total = 100.0 / inLayer.featureCount() if inLayer.featureCount() else 0
features = source.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([inIdx]))
total = 100.0 / source.featureCount() if source.featureCount() else 0
for current, inFeat in enumerate(features):
if feedback.isCanceled():
break

inGeom = inFeat.geometry()
inID = str(inFeat.attributes()[inIdx])
featList = index.nearestNeighbor(inGeom.asPoint(), nPoints)
distList = []
vari = 0.0
request = QgsFeatureRequest().setFilterFids(featList).setSubsetOfAttributes([outIdx])
for outFeat in targetLayer.getFeatures(request):
request = QgsFeatureRequest().setFilterFids(featList).setSubsetOfAttributes([outIdx]).setDestinationCrs(source.sourceCrs())
for outFeat in target_source.getFeatures(request):
if feedback.isCanceled():
break

outID = outFeat.attributes()[outIdx]
outGeom = outFeat.geometry()
dist = distArea.measureLine(inGeom.asPoint(),
outGeom.asPoint())

if matType == 0:
self.writer.addRecord([inID, str(outID), str(dist)])
out_feature = QgsFeature()
out_geom = QgsGeometry.unaryUnion([inFeat.geometry(), outFeat.geometry()])
out_feature.setGeometry(out_geom)
out_feature.setAttributes([inID, outID, dist])
sink.addFeature(out_feature, QgsFeatureSink.FastInsert)
else:
distList.append(float(dist))

@@ -161,44 +195,61 @@ def linearMatrix(self, context, inLayer, inField, targetLayer, targetField,
for i in distList:
vari += (i - mean) * (i - mean)
vari = math.sqrt(vari / len(distList))
self.writer.addRecord([inID, str(mean),
str(vari), str(min(distList)),
str(max(distList))])

out_feature = QgsFeature()
out_feature.setGeometry(inFeat.geometry())
out_feature.setAttributes([inID, mean, vari, min(distList), max(distList)])
sink.addFeature(out_feature, QgsFeatureSink.FastInsert)

feedback.setProgress(int(current * total))

def regularMatrix(self, context, inLayer, inField, targetLayer, targetField,
return {self.OUTPUT: dest_id}

def regularMatrix(self, parameters, context, source, inField, target_source, targetField,
nPoints, feedback):
index = QgsProcessingUtils.createSpatialIndex(targetLayer, context)

inIdx = inLayer.fields().lookupField(inField)
index = QgsSpatialIndex(target_source.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(source.sourceCrs())), feedback)
inIdx = source.fields().lookupField(inField)

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

first = True
features = QgsProcessingUtils.getFeatures(inLayer, context)
total = 100.0 / inLayer.featureCount() if inLayer.featureCount() else 0
sink = None
dest_id = None
features = source.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([inIdx]))
total = 100.0 / source.featureCount() if source.featureCount() else 0
for current, inFeat in enumerate(features):
if feedback.isCanceled():
break

inGeom = inFeat.geometry()
inID = str(inFeat.attributes()[inIdx])
featList = index.nearestNeighbor(inGeom.asPoint(), nPoints)
if first:
first = False
data = ['ID']
fields = QgsFields()
input_id_field = source.fields()[inIdx]
input_id_field.setName('ID')
fields.append(input_id_field)
for i in range(len(featList)):
data.append('DIST_{0}'.format(i + 1))
self.writer.addRecord(data)
fields.append(QgsField('DIST_{0}'.format(i + 1), QVariant.Double))
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
fields, source.wkbType(), source.sourceCrs())

data = [inID]
for i in featList:
request = QgsFeatureRequest().setFilterFid(i)
outFeat = next(targetLayer.getFeatures(request))
outGeom = outFeat.geometry()
for target in target_source.getFeatures(QgsFeatureRequest().setSubsetOfAttributes([]).setFilterFids(featList).setDestinationCrs(source.sourceCrs())):
if feedback.isCanceled():
break
outGeom = target.geometry()
dist = distArea.measureLine(inGeom.asPoint(),
outGeom.asPoint())
data.append(str(float(dist)))
self.writer.addRecord(data)

data.append(float(dist))
out_feature = QgsFeature()
out_feature.setGeometry(inGeom)
out_feature.setAttributes(data)
sink.addFeature(out_feature, QgsFeatureSink.FastInsert)
feedback.setProgress(int(current * total))

return {self.OUTPUT: dest_id}
@@ -69,6 +69,7 @@
from .MeanCoords import MeanCoords
from .Merge import Merge
from .NearestNeighbourAnalysis import NearestNeighbourAnalysis
from .PointDistance import PointDistance
from .PointsInPolygon import PointsInPolygon
from .PointsLayerFromTable import PointsLayerFromTable
from .PolygonsToLines import PolygonsToLines
@@ -93,7 +94,6 @@
from .ZonalStatistics import ZonalStatistics

# from .ExtractByLocation import ExtractByLocation
# from .PointDistance import PointDistance
# from .UniqueValues import UniqueValues
# from .ExportGeometryInfo import ExportGeometryInfo
# from .SinglePartsToMultiparts import SinglePartsToMultiparts
@@ -184,7 +184,7 @@ def __init__(self):

def getAlgs(self):
# algs = [
# UniqueValues(), PointDistance(),
# UniqueValues(),
# ExportGeometryInfo(),
# SinglePartsToMultiparts(),
# ExtractNodes(),
@@ -263,6 +263,7 @@ def getAlgs(self):
MeanCoords(),
Merge(),
NearestNeighbourAnalysis(),
PointDistance(),
PointsInPolygon(),
PointsLayerFromTable(),
PolygonsToLines(),
@@ -0,0 +1,33 @@
<GMLFeatureClassList>
<GMLFeatureClass>
<Name>distance_matrix</Name>
<ElementPath>distance_matrix</ElementPath>
<!--MULTIPOINT-->
<GeometryType>4</GeometryType>
<SRSName>EPSG:4326</SRSName>
<DatasetSpecificInfo>
<FeatureCount>81</FeatureCount>
<ExtentXMin>0.00000</ExtentXMin>
<ExtentXMax>8.00000</ExtentXMax>
<ExtentYMin>-5.00000</ExtentYMin>
<ExtentYMax>3.00000</ExtentYMax>
</DatasetSpecificInfo>
<PropertyDefn>
<Name>InputID</Name>
<ElementPath>InputID</ElementPath>
<Type>String</Type>
<Width>8</Width>
</PropertyDefn>
<PropertyDefn>
<Name>TargetID</Name>
<ElementPath>TargetID</ElementPath>
<Type>String</Type>
<Width>8</Width>
</PropertyDefn>
<PropertyDefn>
<Name>Distance</Name>
<ElementPath>Distance</ElementPath>
<Type>Real</Type>
</PropertyDefn>
</GMLFeatureClass>
</GMLFeatureClassList>

0 comments on commit b7f888b

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