Finding closest points on another polygon

Question

I have three layers:

polygons (having a unique id)
stops (having the polygon id it is laying on) (point)
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 — Ian Turton, Commented May 22, 2020 at 16:53

Val P · Accepted Answer · 2020-05-29 21:55:28Z

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

Kadir Şahbaz · Accepted Answer · 2020-06-01 22:59:05Z

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

# add new fields
if S.fields().indexFromName('near_pl_id') == -1:
    # field for polygonid of the nearest point
    S.addAttribute(QgsField('near_pl_id', QVariant.Int))
if S.fields().indexFromName('near_ps_id') == -1:
    # field for possibilityid of the nearest point
    S.addAttribute(QgsField('near_ps_id', QVariant.Int))

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. — pnz1337, 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. — Kadir Şahbaz, Commented Jun 1, 2020 at 22:55

Kadir Şahbaz · Accepted Answer · 2020-06-08 20:06:13Z

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):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.POSSIBILITY_LYR, self.tr('Possibility layer'), [QgsProcessing.TypeVectorPoint]))
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.STOP_LYR, self.tr('Stop layer'), [QgsProcessing.TypeVectorPoint]))
        self.addParameter(
            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
            sink.addFeature(new_feat, QgsFeatureSink.FastInsert)

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

        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.STOP_LYR, self.tr('Source (Find nearest Points for this Layer)'), [QgsProcessing.TypeVectorPoint]))
        self.addParameter(
            QgsProcessingParameterField(
                self.STOP_IDFIELD, self.tr('Unique ID Field of Source Layer (Any Datatype)'),'ANY','STOP_LYR'))        
        self.addParameter(
            QgsProcessingParameterField(
                self.STOP_POLYGONFIELD, self.tr('Matching ID Field of Source Layer (Numerical)'),'ANY','STOP_LYR',0))  
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.POSSIBILITY_LYR, self.tr('Possibilities (Find Points on this Layer)'), [QgsProcessing.TypeVectorPoint]))
        self.addParameter(
            QgsProcessingParameterField(
                self.POSSIBILITY_IDFIELD, self.tr('Unique Possibility ID Field (Any Datatype)'),'ANY','POSSIBILITY_LYR'))
        self.addParameter(
            QgsProcessingParameterField(
                self.POSSIBILITY_POLYGONFIELD, self.tr('Matching ID Field of Possibilities Layer (Numerical)'),'ANY','POSSIBILITY_LYR',0))
        self.addParameter(
            QgsProcessingParameterEnum(
                self.OPERATION, self.tr('Matching ID Operation (Currently only != and = do work)'), ['!=','=','<','>','<=','>='])) #Only != and = will work here due to expression below
        self.addParameter(
            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
            sink.addFeature(new_feat, QgsFeatureSink.FastInsert)

        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):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.POSSIBILITY_LYR, self.tr('Possibility layer'), [QgsProcessing.TypeVectorPoint]))
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.STOP_LYR, self.tr('Stop layer'), [QgsProcessing.TypeVectorPoint]))
        self.addParameter(
            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():
            sink.addFeature(feature, QgsFeatureSink.FastInsert)

        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
            new_stop_layer.addAttribute(QgsField("near_ps_id", QVariant.Int, len=10))
        if new_stop_layer.fields().indexFromName("near_pl_id") == -1:
            # field for polygonid of the nearest point
            new_stop_layer.addAttribute(QgsField("near_pl_id", QVariant.Int, len=10))

        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. — Kadir Şahbaz, 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? — MrXsquared, Commented Jun 2, 2020 at 19:31
Thanks a lot for your efforts! Solution 1 is perfect for my needs. Now tested with 10k stops and 50k possibilities; took around 50mins. — MrXsquared, 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. :) — Kadir Şahbaz, Commented Jun 3, 2020 at 12:16

MrXsquared · Accepted Answer · 2020-05-30 14:57:32Z

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 — MrXsquared, 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. — J. Monticolo, Commented Jun 1, 2020 at 8:38

Stack Exchange Network

Finding closest points on another polygon

4 Answers 4

Your Answer

Not the answer you're looking for? Browse other questions tagged
qgis
qgis-3
qgis-processing
spatial-analysis
or ask your own question.

Linked

Hot Network Questions

Finding closest points on another polygon

4 Answers 4

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged qgisqgis-3qgis-processingspatial-analysis or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
qgis
qgis-3
qgis-processing
spatial-analysis
or ask your own question.