# Finding closest points on another polygon

I have three layers:

1. polygons (having a unique id)
2. stops (having the polygon id it is laying on) (point)
3. possibilities (having the polygon id it is laying on) (point)

I now want to find the closest possibility-point of each stop-point, which lays on another polygon.

I could easily find the closest point overall, but how could I add the requirement that this point has a different polygon-id than the stop-point itself has?

Here is an example of what I have:

And what I would like to achieve:

I could also use street data from OSM (laying between the buildings) if that would help.

Preferably I would like to be able to integrate this process in a graphical model in QGIS 3.4. But if that's not possible, a different solution than processing tools (like plugins or virtual layers) would do it as well.

• do you already know the polygon id it lies on? if so just exclude that one from the search with an AND otherwise look up the polygon id it lies on (CONTAINS) and then exclude it from the search May 22, 2020 at 16:53

The easiest way to solve your issue is to use an expression that creates the shortest line between the two points.

Here an example of the expression used in the Geometry generator that works perfectly using the algorithm Processing Tools>Vector geometry > Geometry by expression on the Graphical Modeler.

Using the following expression

shortest_line(
$geometry, -- stops aggregate( 'possibilities', 'collect',$geometry,
"polygonid" <> attribute(@parent, 'polygonid')
))

you can create a linestring (the shortest line) that connect the stop point with the closest point of the next polygon, filtering the research to all the polygons except the polygon that intersects the stop point.

The image shown result and setting in the Geometry generator.

Just adding the function end_point to the expression you will have the representation of the selected point.

end_point(shortest_line(
$geometry, -- stops aggregate( 'possibilities', 'collect',$geometry,
"polygonid" <> attribute(@parent, 'polygonid')
)))

This Python script works. But it can be used only in QGIS Python Editor. It is not a Processing script.

from PyQt5.QtCore import QVariant

# layers
P = QgsProject.instance().mapLayersByName("Posibility_Layer")[0]
S = QgsProject.instance().mapLayersByName("Stop_Layer")[0]

S.startEditing()

if S.fields().indexFromName('near_pl_id') == -1:
# field for polygonid of the nearest point
if S.fields().indexFromName('near_ps_id') == -1:
# field for possibilityid of the nearest point

for s in S.getFeatures():
# get points which have different polygon-id than the stop-point itself has
request = 'polygonid != ' + str(s["polygonid"])
distances = {p: s.geometry().distance(p.geometry()) for p in P.getFeatures(request)}

# get the feature which has the minimum distance value
nearest_point = min(distances, key=distances.get)

# you may need to change 'possibilit' as 'possibilityid', (character limit of shapefile)
s['near_ps_id'] = nearest_point["possibilit"]
s['near_pl_id'] = nearest_point["polygonid"]

S.updateFeature(s)

S.commitChanges()

Result:

• I have to note it can be easier with PyQGIS: basically just load the two layers with QgisVectorLayer() class, after just iterate on the two layers - a nested for loop (outer stops, inner possibilities), and a condition, in the inner loop: only calc distance if attributes aren't the same. It should less than 14 rows :) - Answer should be accepted btw. May 29, 2020 at 18:57
• @pnz this script can be used only in QGIS Python Editor. MrXsquared needs a Processing script, so that he can use it in Processing Modeler. I will add a new answer which can be used in a Processing Modeler. Jun 1, 2020 at 22:55

You can add the following scripts as Processing script to Processing Toolbox. Then they can be used in a graphical model. The scripts don't include another Processing tool. Therefore, you can use them in different QGIS 3.x versions without concern. I've tested it in 3.4, 3.8, 3.10, 3.11, 3.12. (Please, also read the comments under processAlgorithm method)

Solution 1:

Creating a new layer, copying the nearest point(NP) in "possibilities layer" to the new layer and adding new fields. I think that is what you need.

from PyQt5.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsField, QgsFeature, QgsProcessing, QgsExpression,
QgsFeatureSink, QgsFeatureRequest, QgsProcessingAlgorithm,
QgsProcessingParameterFeatureSink, QgsProcessingParameterFeatureSource)

class ClosestPointOnAnotherPolygon(QgsProcessingAlgorithm):
POSSIBILITY_LYR = 'POSSIBILITY_LYR'
STOP_LYR = 'STOP_LYR'
OUTPUT = 'OUTPUT'

def initAlgorithm(self, config=None):
QgsProcessingParameterFeatureSource(
self.POSSIBILITY_LYR, self.tr('Possibility layer'), [QgsProcessing.TypeVectorPoint]))
QgsProcessingParameterFeatureSource(
self.STOP_LYR, self.tr('Stop layer'), [QgsProcessing.TypeVectorPoint]))
QgsProcessingParameterFeatureSink(
self.OUTPUT, self.tr('Output Layer'), QgsProcessing.TypeVectorPoint))

def processAlgorithm(self, parameters, context, feedback):
# Get Parameters
possibility_layer = self.parameterAsSource(parameters, self.POSSIBILITY_LYR, context)
stop_layer = self.parameterAsSource(parameters, self.STOP_LYR, context)

fields = possibility_layer.fields()
fields.append(QgsField("stopid", QVariant.Int, len=10))
fields.append(QgsField("join_dist", QVariant.Double, len=20, prec=5))

(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
fields, possibility_layer.wkbType(),
possibility_layer.sourceCrs())

# iterate over stop features
for stop_feat in stop_layer.getFeatures():
# request string for points which have different polygonid
request = QgsFeatureRequest(QgsExpression('polygonid != ' + str(stop_feat["polygonid"])))
distances = {p: stop_feat.geometry().distance(p.geometry())
for p in possibility_layer.getFeatures(request)}

# get the feature which has the minimum distance value
nearest_point = min(distances, key=distances.get)

# create a new feature, set geometry and populate the fields
new_feat = QgsFeature(fields)
new_feat.setGeometry(nearest_point.geometry())
new_feat["possibilid"] = nearest_point["possibilid"]
new_feat["polygonid"] = nearest_point["polygonid"]
new_feat["stopid"] = stop_feat["stopid"]
new_feat["join_dist"] = distances[nearest_point]

# add nearest_point feature to the new layer

return {self.OUTPUT: dest_id}

def tr(self, string):
return QCoreApplication.translate('Processing', string)

def createInstance(self):
return ClosestPointOnAnotherPolygon()

def name(self):
return 'ClosestPointOnAnotherPolygon'

def displayName(self):
return self.tr('Closest Point on Another Polygon')

def group(self):
return self.tr('FROM GISSE')

def groupId(self):
return 'from_gisse'

def shortHelpString(self):
return self.tr('finds closest point on another polygon')

The red one is the new layer.

Here is Solution 1 a little more universal:

from PyQt5.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsField, QgsFeature, QgsProcessing, QgsExpression,
QgsFeatureSink, QgsFeatureRequest, QgsProcessingAlgorithm,
QgsProcessingParameterFeatureSink, QgsProcessingParameterField, QgsProcessingParameterFeatureSource, QgsProcessingParameterEnum)

class ClosestPointWithAttributeCondition(QgsProcessingAlgorithm):
POSSIBILITY_LYR = 'POSSIBILITY_LYR'
POSSIBILITY_IDFIELD = 'POSSIBILITY_IDFIELD'
POSSIBILITY_POLYGONFIELD = 'POSSIBILITY_POLYGONFIELD'
STOP_LYR = 'STOP_LYR'
STOP_IDFIELD = 'STOP_IDFIELD'
STOP_POLYGONFIELD = 'STOP_POLYGONFIELD'
OPERATION = 'OPERATION'
OUTPUT = 'OUTPUT'

def initAlgorithm(self, config=None):

QgsProcessingParameterFeatureSource(
self.STOP_LYR, self.tr('Source (Find nearest Points for this Layer)'), [QgsProcessing.TypeVectorPoint]))
QgsProcessingParameterField(
self.STOP_IDFIELD, self.tr('Unique ID Field of Source Layer (Any Datatype)'),'ANY','STOP_LYR'))
QgsProcessingParameterField(
self.STOP_POLYGONFIELD, self.tr('Matching ID Field of Source Layer (Numerical)'),'ANY','STOP_LYR',0))
QgsProcessingParameterFeatureSource(
self.POSSIBILITY_LYR, self.tr('Possibilities (Find Points on this Layer)'), [QgsProcessing.TypeVectorPoint]))
QgsProcessingParameterField(
self.POSSIBILITY_IDFIELD, self.tr('Unique Possibility ID Field (Any Datatype)'),'ANY','POSSIBILITY_LYR'))
QgsProcessingParameterField(
self.POSSIBILITY_POLYGONFIELD, self.tr('Matching ID Field of Possibilities Layer (Numerical)'),'ANY','POSSIBILITY_LYR',0))
QgsProcessingParameterEnum(
self.OPERATION, self.tr('Matching ID Operation (Currently only != and = do work)'), ['!=','=','<','>','<=','>='])) #Only != and = will work here due to expression below
QgsProcessingParameterFeatureSink(
self.OUTPUT, self.tr('Output Layer'), QgsProcessing.TypeVectorPoint))

def processAlgorithm(self, parameters, context, feedback):
# Get Parameters
possibility_layer = self.parameterAsSource(parameters, self.POSSIBILITY_LYR, context)
possibility_idfield = self.parameterAsFields(parameters, self.POSSIBILITY_IDFIELD, context)
possibility_polygonfield = self.parameterAsFields(parameters, self.POSSIBILITY_POLYGONFIELD, context)
stop_layer = self.parameterAsSource(parameters, self.STOP_LYR, context)
stop_polygonfield = self.parameterAsFields(parameters, self.STOP_POLYGONFIELD, context)
stop_idfield = self.parameterAsFields(parameters, self.STOP_IDFIELD, context)
operation = self.parameterAsString(parameters, self.OPERATION, context)
operationlist = [' != ',' = ',' < ',' > ',' <= ',' >= ']
expressionoperator = str(operationlist[int(operation[0])])

fields = possibility_layer.fields()
fields.append(QgsField(stop_idfield[0]))
fields.append(QgsField("join_dist", QVariant.Double, len=20, prec=5))

(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
fields, possibility_layer.wkbType(),
possibility_layer.sourceCrs())

# iterate over stop features
for stop_feat in stop_layer.getFeatures():
# request string for points which have different polygonid
request = QgsFeatureRequest(QgsExpression(possibility_polygonfield[0] + expressionoperator + str(stop_feat[stop_polygonfield[0]])))
distances = {p: stop_feat.geometry().distance(p.geometry())
for p in possibility_layer.getFeatures(request)}

# get the feature which has the minimum distance value
nearest_point = min(distances, key=distances.get)

# create a new feature, set geometry and populate the fields
new_feat = QgsFeature(fields)
new_feat.setGeometry(nearest_point.geometry())
new_feat[possibility_idfield[0]] = nearest_point[possibility_idfield[0]]
new_feat[possibility_polygonfield[0]] = nearest_point[possibility_polygonfield[0]]
new_feat[stop_idfield[0]] = stop_feat[stop_idfield[0]]
new_feat["join_dist"] = distances[nearest_point]

# add nearest_point feature to the new layer

return {self.OUTPUT: dest_id}

def tr(self, string):
return QCoreApplication.translate('Processing', string)

def createInstance(self):
return ClosestPointWithAttributeCondition()

def name(self):
return 'ClosestPointWithAttributeCondition'

def displayName(self):
return self.tr('Closest Point With Attribute Condition')

def group(self):
return self.tr('FROM GISSE')

def groupId(self):
return 'from_gisse'

def shortHelpString(self):
return self.tr('This Algorithm finds the Sourcelayer`s closest Possibility-Points according the Operation on the Matching ID. The result is an extraction of the Possibilitylayer having the Possibility ID, Matching ID, Source ID and Join Distance.')

Solution 2

You can get NP information by creating a new "stops layer" and adding possibilityid and polygonid of NP to the new layer as new fields.

from PyQt5.QtCore import QCoreApplication, QVariant
from qgis.core import (QgsField, QgsProcessing, QgsExpression,
QgsFeatureSink, QgsFeatureRequest, QgsProcessingUtils,
QgsProcessingAlgorithm, QgsProcessingParameterString,
QgsProcessingParameterFeatureSink, QgsProcessingParameterFeatureSource)

class ClosestPointOnAnotherPolygon2(QgsProcessingAlgorithm):
POSSIBILITY_LYR = 'POSSIBILITY_LYR'
STOP_LYR = 'STOP_LYR'
OUTPUT = 'OUTPUT'

def initAlgorithm(self, config=None):
QgsProcessingParameterFeatureSource(
self.POSSIBILITY_LYR, self.tr('Possibility layer'), [QgsProcessing.TypeVectorPoint]))
QgsProcessingParameterFeatureSource(
self.STOP_LYR, self.tr('Stop layer'), [QgsProcessing.TypeVectorPoint]))
QgsProcessingParameterFeatureSink(
self.OUTPUT, self.tr('Output Layer'), QgsProcessing.TypeVectorPoint))

def processAlgorithm(self, parameters, context, feedback):
# Get Parameters
possibility_layer = self.parameterAsSource(parameters, self.POSSIBILITY_LYR, context)
stop_layer = self.parameterAsSource(parameters, self.STOP_LYR, context)
(sink, dest_id) = self.parameterAsSink(parameters, self.OUTPUT, context,
stop_layer.fields(), stop_layer.wkbType(),
stop_layer.sourceCrs())

# Copy stop features to new stop layer
for feature in stop_layer.getFeatures():

new_stop_layer = QgsProcessingUtils.mapLayerFromString(dest_id, context)
new_stop_layer.startEditing()

# if not exist, add new fields
if new_stop_layer.fields().indexFromName("near_ps_id") == -1:
# field for possibilityid of the nearest point
if new_stop_layer.fields().indexFromName("near_pl_id") == -1:
# field for polygonid of the nearest point

for feature in new_stop_layer.getFeatures():
# get points which have different polygon-id than the stop-point itself has
request = QgsFeatureRequest(QgsExpression('polygonid != ' + str(feature["polygonid"])))
distances = {p: feature.geometry().distance(p.geometry())
for p in possibility_layer.getFeatures(request)}

# get the feature which has the minimum distance value
nearest_point = min(distances, key=distances.get)

# you may need to change 'possibilit' as 'possibilityid', (character limit of shapefile)
feature['near_ps_id'] = nearest_point["possibilit"]
feature['near_pl_id'] = nearest_point["polygonid"]

new_stop_layer.updateFeature(feature)

new_stop_layer.commitChanges()

return {self.OUTPUT: dest_id}

def tr(self, string):
return QCoreApplication.translate('Processing', string)

def createInstance(self):
return ClosestPointOnAnotherPolygon2()

def name(self):
return 'ClosestPointOnAnotherPolygon2'

def displayName(self):
return self.tr('Closest Point on Another Polygon 2')

def group(self):
return self.tr('FROM GISSE')

def groupId(self):
return 'from_gisse'

def shortHelpString(self):
return self.tr('finds closest point on another polygon')

A Sample Model:

• I added as the second answer. Because the script in my first answer is used only in QGIS Python Editor. So I think it must be a separate answer. Jun 1, 2020 at 23:45
• Very neat solution. Thanks a lot! Will do some more testing tomorrow. I don't want to ask for too much here, but do you think it is possible to include the related stopid and join distance in solution 1? Jun 2, 2020 at 19:31
• Yes, it is possible. I've edited the answer. Jun 2, 2020 at 22:02
• Thanks a lot for your efforts! Solution 1 is perfect for my needs. Now tested with 10k stops and 50k possibilities; took around 50mins. Jun 3, 2020 at 11:10
• You're welcome. A feedback can be used to update the progress bar. But I'm new in Processing script in QGIS 3. Its structure is so different than QGIS 2. Your question motivated me to learn Processing script. :) Jun 3, 2020 at 12:16

For fun, I've tested a method with pure QGIS Expressions.

Disclaimer: I didn't test it on very large datasets.

For start, you can see the solution as a point 'geometry generator' (you can use 'geometry by expression' from toolbox to make it actual geometry) symbol layer in the stops point layer. Here the code:

closest_point(
aggregate(
layer:='possibilities',
aggregate:='collect',
expression:=$geometry, filter:="polygonid" <> attribute(@parent, 'polygonid') ),$geometry
)

But it's 100 % visual and we don't see at which stops feature it's attached. So, we can create an other stops symbol layer, line geometry type this time with:

make_line(
$geometry, closest_point( aggregate( layer:='possibilities', aggregate:='collect', expression:=$geometry,
filter:="polygonid" <> attribute(@parent, 'polygonid')
),
$geometry ) ) So, you can even style it as an arrow: the start is the stops point and the end is the closest possibilities point! But here, the real answer is to deal with attributes. Go to the stops field calculator, create a new integer field with the expression below: attribute( array_first( array_filter( aggregate( layer:= 'possibilities', aggregate:= 'array_agg', expression:=$currentfeature
),
geom_to_wkt(geometry(@element)) =
geom_to_wkt(
closest_point(
aggregate(
layer:='possibilities',
aggregate:='collect',
expression:=$geometry, filter:="polygonid" <> attribute(@parent, 'polygonid') ),$geometry
)
)
)
)
, 'possibilityid'
)

I make an array with the possibilities features and filter it with the equality between its geometry WKT and the geometry WKT of the possibilities closest point (you find the point symbol layer code used earlier).

The array filtered return, normally, one possibilities feature (the closest to the current stops feature). I take this one (array_first) and get its possibilityid attribute.

• The first part of your answer works just fine. The second one returns NULL, seems like geometry(@element) does not return anything. Also tried \$geometry, @current_feature, @current_geometry instead, but overall result seems to be always NULL May 30, 2020 at 14:54
• It's strange, I just test it on my test dataset (buildings (polygonid) polygons (not really used here), (possibilities (possibilityid, polygonid) points, stops (stopid, polygonid) points). All the fields are in integer. The projection is EPSG: 3857. Jun 1, 2020 at 8:38