Skip to content

Commit

Permalink
Rework Select by Location algorithm
Browse files Browse the repository at this point in the history 
Changes:

- handle different CRS transparently

- don't build a spatial index on the selection layer. Instead
only use feature requests to fetch features which are within
the desired bounds, and rely on the presence of an appropriate
spatial index at the provider's backend. Otherwise, we force
every user of this algorithm to have a full iteration of the
source table, regardless of how large the table is. That means
that trying to select a set of addresses which fall within
a specific locality from a table which contains the addresses
for a whole state will FORCE every address in the state to
be initially read before any calculation begins. With this
change only those features within the bounding box of the
selected localities will ever be fetched from the provider,
resulting in huge speed improvements for the algorithm.

- use prepared geometries for the spatial relation tests.
This dramatically speeds up the algorithm in the case
where the intersection layer features cover multiple
features from the 'selection' layer.

- Add a 'select within current selection' mode

- Optimise feature requests for efficiency (especially
with respect to the 'disjoint' selection mode)
  • Loading branch information
@nyalldawson
nyalldawson committed Sep 8, 2017
1 parent 20d8244 commit 24a4ab7
Show file tree
Hide file tree
Showing 2 changed files with 81 additions and 63 deletions.
4 changes: 2 additions & 2 deletions python/plugins/processing/algs/qgis/QGISAlgorithmProvider.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,6 +134,7 @@
from .SaveSelectedFeatures import SaveSelectedFeatures
from .SelectByAttribute import SelectByAttribute
from .SelectByExpression import SelectByExpression
from .SelectByLocation import SelectByLocation
from .ServiceAreaFromLayer import ServiceAreaFromLayer
from .ServiceAreaFromPoint import ServiceAreaFromPoint
from .SetMValue import SetMValue
Expand Down Expand Up @@ -167,7 +168,6 @@
from .ZonalStatistics import ZonalStatistics

# from .ExtractByLocation import ExtractByLocation
# from .SelectByLocation import SelectByLocation
# from .SpatialJoin import SpatialJoin

pluginPath = os.path.normpath(os.path.join(
Expand All @@ -183,7 +183,6 @@ def __init__(self):

def getAlgs(self):
# algs = [
# SelectByLocation(),
# ExtractByLocation(),
# SpatialJoin(),
# ]
Expand Down Expand Up @@ -281,6 +280,7 @@ def getAlgs(self):
SaveSelectedFeatures(),
SelectByAttribute(),
SelectByExpression(),
SelectByLocation(),
ServiceAreaFromLayer(),
ServiceAreaFromPoint(),
SetMValue(),
Expand Down
140 changes: 79 additions & 61 deletions python/plugins/processing/algs/qgis/SelectByLocation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,13 +29,19 @@

from qgis.PyQt.QtGui import QIcon

from qgis.core import QgsGeometry, QgsFeatureRequest, QgsProcessingUtils
from qgis.core import (QgsGeometry,
QgsFeatureRequest,
QgsProcessingUtils,
QgsProcessing,
QgsProcessingParameterVectorLayer,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterEnum,
QgsProcessingParameterNumber,
QgsProcessingOutputVectorLayer,
QgsVectorLayer)


from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
from processing.core.parameters import ParameterSelection
from processing.core.parameters import ParameterVector
from processing.core.parameters import ParameterNumber
from processing.core.outputs import OutputVector
from processing.tools import vector

pluginPath = os.path.split(os.path.split(os.path.dirname(__file__))[0])[0]
Expand Down Expand Up @@ -63,32 +69,43 @@ def initAlgorithm(self, config=None):
self.predicates = (
('intersects', self.tr('intersects')),
('contains', self.tr('contains')),
('disjoint', self.tr('disjoint')),
('equals', self.tr('equals')),
('disjoint', self.tr('is disjoint')),
('isEqual', self.tr('equals')),
('touches', self.tr('touches')),
('overlaps', self.tr('overlaps')),
('within', self.tr('within')),
('crosses', self.tr('crosses')))

self.reversed_predicates = {'intersects': 'intersects',
'contains': 'within',
'disjoint': 'disjoint',
'isEqual': 'isEqual',
'touches': 'touches',
'overlaps': 'overlaps',
'within': 'contains',
'crosses': 'crosses'}

self.methods = [self.tr('creating new selection'),
self.tr('adding to current selection'),
self.tr('select within current selection'),
self.tr('removing from current selection')]

self.addParameter(ParameterVector(self.INPUT,
self.tr('Layer to select from')))
self.addParameter(ParameterVector(self.INTERSECT,
self.tr('Additional layer (intersection layer)')))
self.addParameter(ParameterSelection(self.PREDICATE,
self.tr('Geometric predicate'),
self.predicates,
multiple=True))
self.addParameter(ParameterNumber(self.PRECISION,
self.tr('Precision'),
0.0, None, 0.0))
self.addParameter(ParameterSelection(self.METHOD,
self.tr('Modify current selection by'),
self.methods, 0))
self.addOutput(OutputVector(self.OUTPUT, self.tr('Selected (location)'), True))
self.addParameter(QgsProcessingParameterVectorLayer(self.INPUT,
self.tr('Select features from'), types=[QgsProcessing.TypeVectorAnyGeometry]))
self.addParameter(QgsProcessingParameterEnum(self.PREDICATE,
self.tr('Where the features are (geometric predicate)'),
options=[p[1] for p in self.predicates],
allowMultiple=True, defaultValue=[0]))
self.addParameter(QgsProcessingParameterFeatureSource(self.INTERSECT,
self.tr('By comparing to the features from'), types=[QgsProcessing.TypeVectorAnyGeometry]))
self.addParameter(QgsProcessingParameterNumber(self.PRECISION,
self.tr('Precision'), type=QgsProcessingParameterNumber.Double,
minValue=0.0, defaultValue=0.0))
self.addParameter(QgsProcessingParameterEnum(self.METHOD,
self.tr('Modify current selection by'),
options=self.methods, defaultValue=0))

self.addOutput(QgsProcessingOutputVectorLayer(self.OUTPUT, self.tr('Selected (by location)')))

def name(self):
return 'selectbylocation'
Expand All @@ -97,60 +114,61 @@ def displayName(self):
return self.tr('Select by location')

def processAlgorithm(self, parameters, context, feedback):
filename = self.getParameterValue(self.INPUT)
inputLayer = QgsProcessingUtils.mapLayerFromString(filename, context)
method = self.getParameterValue(self.METHOD)
filename2 = self.getParameterValue(self.INTERSECT)
selectLayer = QgsProcessingUtils.mapLayerFromString(filename2, context)
predicates = self.getParameterValue(self.PREDICATE)
precision = self.getParameterValue(self.PRECISION)

oldSelection = set(inputLayer.selectedFeatureIds())
inputLayer.removeSelection()
index = QgsProcessingUtils.createSpatialIndex(inputLayer, context)
select_layer = self.parameterAsVectorLayer(parameters, self.INPUT, context)
method = QgsVectorLayer.SelectBehavior(self.parameterAsEnum(parameters, self.METHOD, context))
intersect_source = self.parameterAsSource(parameters, self.INTERSECT, context)
# build a list of 'reversed' predicates, because in this function
# we actually test the reverse of what the user wants (allowing us
# to prepare geometries and optimise the algorithm)
predicates = [self.reversed_predicates[self.predicates[i][0]] for i in self.parameterAsEnums(parameters, self.PREDICATE, context)]
precision = self.parameterAsDouble(parameters, self.PRECISION, context)

if 'disjoint' in predicates:
disjoinSet = []
for feat in QgsProcessingUtils.getFeatures(inputLayer, context):
disjoinSet.append(feat.id())

geom = QgsGeometry()
selectedSet = []
features = QgsProcessingUtils.getFeatures(selectLayer, context)
total = 100.0 / selectLayer.featureCount() if selectLayer.featureCount() else 0
disjoint_set = select_layer.allFeatureIds()
else:
disjoint_set = None

selected_set = set()
request = QgsFeatureRequest().setSubsetOfAttributes([]).setDestinationCrs(select_layer.crs())
features = intersect_source.getFeatures(request)
total = 100.0 / intersect_source.featureCount() if intersect_source.featureCount() else 0
for current, f in enumerate(features):
geom = vector.snapToPrecision(f.geometry(), precision)
bbox = geom.boundingBox()
if feedback.isCanceled():
break

if not f.hasGeometry():
continue

engine = QgsGeometry.createGeometryEngine(f.geometry().geometry())
engine.prepareGeometry()
bbox = f.geometry().boundingBox()
bbox.grow(0.51 * precision)
intersects = index.intersects(bbox)

request = QgsFeatureRequest().setFilterFids(intersects).setSubsetOfAttributes([])
for feat in inputLayer.getFeatures(request):
tmpGeom = vector.snapToPrecision(feat.geometry(), precision)
request = QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry).setFilterRect(bbox).setSubsetOfAttributes([])
for test_feat in select_layer.getFeatures(request):
if feedback.isCanceled():
break

if test_feat in selected_set:
# already added this one, no need for further tests
continue

res = False
for predicate in predicates:
if predicate == 'disjoint':
if tmpGeom.intersects(geom):
if test_feat.geometry().intersects(f.geometry()):
try:
disjoinSet.remove(feat.id())
disjoint_set.remove(test_feat.id())
except:
pass # already removed
else:
res = getattr(tmpGeom, predicate)(geom)
if res:
selectedSet.append(feat.id())
if getattr(engine, predicate)(test_feat.geometry().geometry()):
selected_set.add(test_feat.id())
break

feedback.setProgress(int(current * total))

if 'disjoint' in predicates:
selectedSet = selectedSet + disjoinSet

if method == 1:
selectedSet = list(oldSelection.union(selectedSet))
elif method == 2:
selectedSet = list(oldSelection.difference(selectedSet))
selected_set = list(selected_set) + disjoint_set

inputLayer.selectByIds(selectedSet)
self.setOutputValue(self.OUTPUT, filename)
select_layer.selectByIds(list(selected_set), method)
return {self.OUTPUT: parameters[self.INPUT]}

0 comments on commit 24a4ab7

Please sign in to comment.