Encontrar los puntos más cercanos en otro polígono

Question

Tengo tres capas:

1. polígonos (con un identificador único)
2. se detiene (con el identificador del polígono sobre el que se encuentra) (punto)
3. posibilidades (con el identificador del polígono sobre el que se encuentra) (punto)

Ahora quiero encontrar el punto de posibilidad más cercano de cada punto de parada, que se encuentra en otro polígono.

Podría encontrar fácilmente el punto más cercano en general, pero ¿cómo podría añadir el requisito de que este punto tenga un polígono-id diferente al que tiene el propio punto de parada?

He aquí un ejemplo de lo que tengo:

Y lo que me gustaría conseguir:

También podría utilizar los datos de la calle de OSM (tendido entre los edificios) si eso ayudara.

Preferiblemente me gustaría poder integrar este proceso en un modelo gráfico en QGIS 3.4. Pero si eso no es posible, una solución diferente a las herramientas de procesamiento (como plugins o capas virtuales) lo haría también.

Answer 1

La forma más sencilla de resolver tu problema es utilizar una expresión que cree la línea más corta entre los dos puntos.

Aquí un ejemplo de la expresión utilizada en el generador de Geometría que funciona perfectamente utilizando el algoritmo Processing Tools>Geometría vectorial > Geometría por expresión en el Modelador Gráfico.

Utilizando la siguiente expresión

shortest_line(
$geometry, -- stops

aggregate(
'possibilities',
'collect',
$geometry,
"polygonid" <> attribute(, 'polygonid')
))

puede crear una línea (la línea más corta) que conecte el punto de parada con el punto más cercano del siguiente polígono, filtrando la investigación a todos los polígonos excepto al polígono que interseca el punto de parada.

La imagen mostrada resultado y ajuste en el generador de Geometría.

Basta con añadir la función end_point a la expresión tendrás la representación del punto seleccionado.

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

Answer 2

Este script Python funciona. Pero sólo se puede utilizar en QGIS Python Editor. No es un script de Processing.

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

Resultado:

Answer 3

Puede añadir las siguientes secuencias de comandos como Script de procesamiento a Processing Toolbox. A continuación, se pueden utilizar en un modelo gráfico. Los scripts no incluyen otra herramienta de Processing. Por lo tanto, puedes usarlos en diferentes versiones de QGIS 3.x sin problemas. Lo he probado en 3.4, 3.8, 3.10, 3.11, 3.12. (Por favor, lea también los comentarios bajo processAlgorithm método)

Solución 1:

Creando una nueva capa, copiando el punto(NP) más cercano en la "capa de posibilidades" a la nueva capa y añadiendo nuevos campos. Creo que eso es lo que necesita.

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

La roja es la nueva capa.

Aquí está la Solución 1 un poco más 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.')

Solución 2

Puede obtener información sobre NP creando una nueva "capa de paradas" y añadiendo possibilityid y polygonid de NP a la nueva capa como nuevos campos.

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

Un modelo de ejemplo:

Answer 4

Por diversión, he probado un método con puras Expresiones QGIS.

Descargo de responsabilidad : No lo he probado con conjuntos de datos muy grandes.

---

Para empezar, puede ver la solución como un punto "generador de geometría" (puede utilizar "geometría por expresión" de la caja de herramientas para convertirla en geometría real) en la capa de símbolos stops capa de puntos. Aquí el código:

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

---

Pero es 100% visual y no vemos en que stops característica que se adjunta. Por lo tanto, podemos crear otro stops capa de símbolo, tipo de geometría de línea esta vez con:

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

Por lo tanto, incluso se puede estilizar como una flecha: el inicio es el stops punto y el final es el más cercano possibilities ¡punto!

---

Pero aquí, la verdadera respuesta es ocuparse de los atributos.

Ir a la stops cree un nuevo campo entero con la siguiente expresión:

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(, 'polygonid')
          ),
          $geometry
        )
      )
    )
  )
  , 'possibilityid'
)

Hago un array con el possibilities y filtrarlo con la igualdad entre su geometría WKT y la geometría WKT del possibilities punto más cercano (se encuentra el código de capa de símbolo de punto utilizado anteriormente).

La matriz filtrada devuelve, normalmente, una possibilities (la más cercana a la actual stops función). Tomo esta ( array_first ) y obtener su possibilityid atributo.

Answer 5

Execute SQL debe estar fácilmente disponible para el constructor del modelo gráfico, y ejecutar

SELECT  a.id AS stop_id,
        b.id AS poss_id,
        MIN(Distance(a.geometry, b.geometry))
FROM    <stops> AS a, <possibilites> AS b
WHERE   a.<polygonId> <> b.<polygonId>
GROUP BY
        1
ORDER BY
        3
;

produce una capa con stop_id el más cercano poss_id con cualquier otro <polygonId> y la distancia ( en unidades CRS ver SpatiaLite función de referencia para otras opciones) entre ellos.

En la práctica (pero con sólo 5 minutos de tiempo invertido), no he sido capaz de pas bloqueo de QGIS (3.12) durante el uso autónomo Execute SQL ...

---

Nota:

• espacial (K) Vecino más próximo las búsquedas no se pueden optimizar de ninguna manera en la versión de SQLite con la que se compila QGIS; la sintaxis anterior funcionará significativamente peor en grandes conjuntos de datos que, por ejemplo, una implementación basada en índices (por ejemplo, PostGIS o incluso algunas de las herramientas integradas).