# 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 Commented 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. Commented 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. Commented 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. Commented 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? Commented Jun 2, 2020 at 19:31
• Yes, it is possible. I've edited the answer. Commented 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. Commented 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. :) Commented 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` Commented 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. Commented Jun 1, 2020 at 8:38