¿Proximidad entre la capa de puntos y la capa de líneas?

Question

He hecho esta pregunta varias veces en stackoverflow e irc entre #qgis y #postgis y también he intentado codificarlo o implementarlo yo mismo en postgis sin respuesta real.

Mediante programación (preferiblemente python), me gustaría dibujar una línea desde una capa de puntos, hasta su proyección en la línea más cercana de una capa de líneas o polígonos.

Por ahora, la mayoría de mis datos están en los formatos shape y postgis de ESRI; sin embargo, prefiero mantenerme alejado de una solución postgis, ya que soy un usuario predominantemente de shp + qgis.

Una solución ideal sería implementar GDAL/OGR con bibliotecas de python o similares

• Utilizando las bibliotecas GDAL/OGR, ¿por dónde debería empezar? ¿Sería posible dar un plan de solución?
• ¿Puedo utilizar NetworkX para realizar el análisis del vecino más cercano?
• ¿Es esto realmente posible?

Si es más fácil, los puntos podrían conectarse al punto final del segmento en lugar de un punto proyectado

Answer 1

Véase el comentario de abajo sobre cómo mi respuesta no debe considerarse una solución fiable... Voy a dejar este post original aquí sólo para que otros puedan examinar el problema.

Si entiendo la pregunta, este procedimiento general debería funcionar.

Encontrar el camino más corto entre un punto (definido por x,y o x,y,z) y un polinomio (definido por un conjunto de conexión de x,y o x,y,z) dentro del espacio euclidiano:

1) A partir de un punto definido por el usuario (lo llamaré pt0), encuentra el vértice más cercano de la polilínea (pt1). OGRinfo puede sondear los vértices de una polilínea, y entonces los cálculos de distancia pueden hacerse mediante métodos estándar. Por ejemplo, iterar sobre un cálculo de distancia como: distancia_en_radios=2*math.asin(math.sqrt(math.pow((math.sin((pt0_radians-ptx_radians)/2)),2) + math.cos(pt0_radianes)*math.cos(ptx_radians)*math.pow((math.sin((pt0_radians)/2)),2))

2) Almacenar el valor de la distancia mínima asociada (d1) y (pt1)

3) mira los dos segmentos que se alejan de pt1 (en el ogrinfo linestring, estos serán los vértices anteriores y posteriores). Registra estos vértices (n2 y n3).

4) crear la fórmula y = mx + b para cada segmento

5) Relaciona tu punto (pt0) con la perpendicular para cada una de esas dos fórmulas

6) Calcula la distancia y las intersecciones (d2 y d3; pt2, pt3)

7) Compara las tres distancias (d1, d2, d3) para la más corta. Su pt0 al nodo asociado (pt1, pt2 o pt3) es el enlace más corto.

Esa es una respuesta de flujo de conciencia - espero que mi imagen mental del problema y la solución encaje.

Answer 2

Aquí hay un script en python para QGIS >2.0 hecho a partir de las pistas y soluciones dadas anteriormente. Funciona bien para una cantidad razonable de puntos y líneas. Pero no lo he probado con una gran cantidad de objetos.

Por supuesto había que copiarlo en idle o cualquier otra solución menos "pitónica" y guardarlo como "closest.point.py".

En la caja de herramientas de QGIS vaya a script, herramientas, añada un script, y elíjalo.

##Vector=group
##CLosest_Point_V2=name
##Couche_de_Points=vector
##Couche_de_Lignes=vector

"""
This script intent to provide a count as for the SQL Funciton CLosestPoint
Ce script vise a recréer dans QGIS la Focntion SQL : CLosest Point
It rely on the solutions provided in "Nearest neighbor between a point layer and a line layer"
  http://gis.stackexchange.com/questions/396/nearest-pojected-point-from-a-point-                               layer-on-a-line-or-polygon-outer-ring-layer
V2 du  8 aout 2016
jean-christophe.baudin@onema.fr
"""
from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import os
import sys
import unicodedata 
from osgeo import ogr
from math import sqrt
from sys import maxint

from processing import *

def magnitude(p1, p2):
    if p1==p2: return 1
    else:
        vect_x = p2.x() - p1.x()
        vect_y = p2.y() - p1.y()
        return sqrt(vect_x**2 + vect_y**2)

def intersect_point_to_line(point, line_start, line_end):
    line_magnitude =  magnitude(line_end, line_start)
    u = ((point.x()-line_start.x())*(line_end.x()-line_start.x())+(point.y()-line_start.y())*(line_end.y()-line_start.y()))/(line_magnitude**2)
    # closest point does not fall within the line segment, 
    # take the shorter distance to an endpoint
    if u < 0.0001 or u > 1:
        ix = magnitude(point, line_start)
        iy = magnitude(point, line_end)
        if ix > iy:
            return line_end
        else:
            return line_start
    else:
        ix = line_start.x() + u * (line_end.x() - line_start.x())
        iy = line_start.y() + u * (line_end.y() - line_start.y())
        return QgsPoint(ix, iy)

layerP = processing.getObject(Couche_de_Points)
providerP = layerP.dataProvider()
fieldsP = providerP.fields()
inFeatP = QgsFeature()

layerL = processing.getObject(Couche_de_Lignes)
providerL = layerL.dataProvider()
fieldsL = providerL.fields()
inFeatL = QgsFeature()

counterP = counterL= nElement=0

for featP in layerP.selectedFeatures():
    counterP+=1
if counterP==0:
    QMessageBox.information(None,"information:","Choose at least one point from point layer_"+ str(layerP.name())) 

indexLine=QgsSpatialIndex()
for featL in layerL.selectedFeatures():
    indexLine.insertFeature(featL)
    counterL+=1
if counterL==0:
    QMessageBox.information(None,"information:","Choose at least one line from point layer_"+ str(layerL.name()))
    #QMessageBox.information(None,"DEBUGindex:",str(indexBerge))     
ClosestP=QgsVectorLayer("Point", "Projected_ Points_From_"+ str(layerP.name()), "memory")
QgsMapLayerRegistry.instance().addMapLayer(ClosestP)
prClosestP = ClosestP.dataProvider()

for f in fieldsP:
    znameField= f.name()
    Type= str(f.typeName())
    if Type == 'Integer': prClosestP.addAttributes([ QgsField( znameField, QVariant.Int)])
    if Type == 'Real': prClosestP.addAttributes([ QgsField( znameField, QVariant.Double)])
    if Type == 'String': prClosestP.addAttributes([ QgsField( znameField, QVariant.String)])
    else : prClosestP.addAttributes([ QgsField( znameField, QVariant.String)])
prClosestP.addAttributes([QgsField("DistanceP", QVariant.Double),
                                        QgsField("XDep", QVariant.Double),
                                        QgsField("YDep", QVariant.Double),
                                        QgsField("XProj", QVariant.Double),
                                        QgsField("YProj", QVariant.Double),
                                        QgsField("Xmed", QVariant.Double),
                                        QgsField("Ymed", QVariant.Double)])
featsP = processing.features(layerP)
nFeat = len(featsP)
"""
for inFeatP in featsP:
    progress.setPercentage(int(100 * nElement / nFeatL))
    nElement += 1
    # pour avoir l'attribut d'un objet/feat .... 
    attributs = inFeatP.attributes()
"""

for inFeatP in layerP.selectedFeatures():
    progress.setPercentage(int(100 * nElement / counterL))
    nElement += 1
    attributs=inFeatP.attributes()
    geomP=inFeatP.geometry()
    nearest_point = None
    minVal=0.0
    counterSelec=1
    first= True
    nearestsfids=indexLine.nearestNeighbor(geomP.asPoint(),counterSelec)
    #http://blog.vitu.ch/10212013-1331/advanced-feature-requests-qgis
    #layer.getFeatures( QgsFeatureRequest().setFilterFid( fid ) )
    request = QgsFeatureRequest().setFilterFids( nearestsfids )
    #list = [ feat for feat in CoucheL.getFeatures( request ) ]
    # QMessageBox.information(None,"DEBUGnearestIndex:",str(list))
    NBNodes=0
    Dist=DistT=minValT=Distance=0.0
    for featL in  layerL.getFeatures(request):
        geomL=featL.geometry()
        firstM=True
        geomL2=geomL.asPolyline()
        NBNodes=len(geomL2)
        for i in range(1,NBNodes):
            lineStart,lineEnd=geomL2[i-1],geomL2[i]
            ProjPoint=intersect_point_to_line(geomP.asPoint(),QgsPoint(lineStart),QgsPoint(lineEnd))
            Distance=magnitude(geomP.asPoint(),ProjPoint)
            toto=''
            toto=toto+ 'lineStart :'+ str(lineStart)+ '  lineEnd : '+ str(lineEnd)+ '\n'+ '\n'
            toto=toto+ 'ProjPoint '+ str(ProjPoint)+ '\n'+ '\n'
            toto=toto+ 'Distance '+ str(Distance)
            #QMessageBox.information(None,"DEBUG", toto)
            if firstM:
                minValT,nearest_pointT,firstM = Distance,ProjPoint,False
            else:
                if Distance < minValT:
                    minValT=Distance
                    nearest_pointT=ProjPoint
            #at the end of the loop save the nearest point for a line object
            #min_dist=magnitude(ObjetPoint,PProjMin)
            #QMessageBox.information(None,"DEBUG", " Dist min: "+ str(minValT))
        if first:
            minVal,nearest_point,first = minValT,nearest_pointT,False
        else:
            if minValT < minVal:
                minVal=minValT
                nearest_point=nearest_pointT
                #at loop end give the nearest Projected points on Line nearest Line
    PProjMin=nearest_point
    Geom= QgsGeometry().fromPoint(PProjMin)
    min_dist=minVal
    PX=geomP.asPoint().x()
    PY=geomP.asPoint().y()
    Xmed=(PX+PProjMin.x())/2
    Ymed=(PY+PProjMin.y())/2
    newfeat = QgsFeature()
    newfeat.setGeometry(Geom)
    Values=[]
    #Values.extend(attributs)
    fields=layerP.pendingFields()
    Values=[attributs[i] for i in range(len(fields))]
    Values.append(min_dist)
    Values.append(PX)
    Values.append(PY)
    Values.append(PProjMin.x())
    Values.append(PProjMin.y())
    Values.append(Xmed)
    Values.append(Ymed)
    newfeat.setAttributes(Values)
    ClosestP.startEditing()  
    prClosestP.addFeatures([ newfeat ])
    #prClosestP.updateExtents()
ClosestP.commitChanges()
iface.mapCanvas().refresh()

¡¡¡!!! ¡¡¡ADVERTENCIA !!! Preste atención a que algunos puntos proyectados "raros" / erróneos pueden producirse debido a este comando de línea:

nearestsfids=indexLine.nearestNeighbor(geomP.asPoint(),counterSelec)

Le site counterSelec en él se establece cuántos nearestNeighbor se devuelven. De hecho, cada punto debería proyectarse a la menor distancia posible de cada objeto de línea; y la distancia mínima encontrada daría la línea correcta y el punto proyectado como los vecinos más cercanos que buscamos. Para reducir el tiempo de bucle, se utiliza el comando Vecino más cercano. La elección de un counterSelec reducir a 1 devolverá el "primer" objeto encontrado (su caja delimitadora más exactamente) y puede no ser el correcto. Los objetos de diferente tamaño de línea pueden obligar a elegir pueden ser 3 o 5, o incluso más objetos más cercanos con el fin de determinar la distancia más corta. Cuanto más alto sea el valor, más tiempo se necesita. Con cientos de puntos y líneas empieza a ser muy lento con 3 o 5 vecinos más cercanos, con miles puede fallar con esos valores.