Nächster Nachbar zwischen Punktebene und Linienebene? [geschlossen]

Ich habe diese Frage mehrmals in Bezug auf Stackoverflow und IRC zwischen #qgis und #postgis gestellt und auch versucht, sie zu codieren oder selbst in postgis zu implementieren, ohne eine echte Antwort zu finden.

Mit der Programmierung (am besten mit Python) möchte ich eine Linie von einer Punktebene bis zu ihrer Projektion auf die nächste Linie einer Linie oder Polygonebene zeichnen.

Ab sofort liegen die meisten meiner Daten in den Formaten ESRI und Postgis vor. ich möchte mich jedoch lieber von einer postgis-lösung fernhalten, da ich überwiegend shp + qgis-benutzer bin.

Eine ideale Lösung wäre die Implementierung von GDAL / OGR mit Python oder ähnlichen Bibliotheken

• Wo soll ich mit der Verwendung der GDAL / OGR-Bibliotheken beginnen? Wäre es möglich, einen Lösungsplan zu erstellen?
• Kann ich mit NetworkX die Analyse des nächsten Nachbarn durchführen?
• Ist das überhaupt möglich?

Wenn es einfacher ist, können die Punkte anstelle eines projizierten Punkts mit dem Segmentendpunkt verbunden werden

Diese Frage stellte sich als etwas kniffliger heraus, als ich gedacht hatte. Es gibt viele Implementierungen für den kürzesten Abstand selbst, z. B. den Shapely Provided Distance (von GEOS). Nur wenige der Lösungen geben den Schnittpunkt selbst an, aber nur die Entfernung.

Bei meinem ersten Versuch wurde der Punkt durch den Abstand zwischen dem Punkt und dem Polygon gepuffert und nach Schnittpunkten gesucht, aber Rundungsfehler verhindern, dass dies eine genaue Antwort ergibt.

Hier ist eine vollständige Lösung mit Shapely, basierend auf diesen Gleichungen :

#!/usr/bin/env python
from shapely.geometry import Point, Polygon
from math import sqrt
from sys import maxint

# define our polygon of interest, and the point we'd like to test
# for the nearest location
polygon = Polygon(((0, 0), (0, 1), (1, 1), (1, 0), (0, 0)))
point = Point(0.5, 1.5)

# pairs iterator:
# http://stackoverflow.com/questions/1257413/1257446#1257446
def pairs(lst):
    i = iter(lst)
    first = prev = i.next()
    for item in i:
        yield prev, item
        prev = item
    yield item, first

# these methods rewritten from the C version of Paul Bourke's
# geometry computations:
# http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/
def magnitude(p1, p2):
    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.00001 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 Point([ix, iy])

nearest_point = None
min_dist = maxint

for seg_start, seg_end in pairs(list(polygon.exterior.coords)[:-1]):
    line_start = Point(seg_start)
    line_end = Point(seg_end)

    intersection_point = intersect_point_to_line(point, line_start, line_end)
    cur_dist =  magnitude(point, intersection_point)

    if cur_dist < min_dist:
        min_dist = cur_dist
        nearest_point = intersection_point

print "Closest point found at: %s, with a distance of %.2f units." % \
   (nearest_point, min_dist)

Für die Nachwelt sieht es so aus, als würde diese ArcView-Erweiterung dieses Problem ganz gut lösen. Schade, dass sie auf einer toten Plattform in einer toten Sprache geschrieben ist ...

Eine PostGIS-Antwort (für Multilinestring, wenn Linestring, entfernen Sie die Funktion st_geometryn)

select t2.gid as point_gid, t1.gid as line_gid, 
st_makeline(t2.geom,st_line_interpolate_point(st_geometryn(t1.geom,1),st_line_locate_point(st_geometryn(t1.geom,1),t2.geom))) as geom
from your_line_layer t1, your_point_layer t2, 
(
select gid as point_gid, 
(select gid 
from your_line_layer
order by st_distance(your_line_layer.geom, your_point_layer.geom)
limit 1 ) as line_gid
from your_point_layer
) as t3
where t1.gid = t3.line_gid
and t2.gid = t3.point_gid

Dies ist ein bisschen alt, aber ich habe heute nach Lösungen für dieses Problem gesucht (Punkt -> Linie). Die einfachste Lösung für dieses Problem ist:

>>> from shapely.geometry import Point, LineString
>>> line = LineString([(0, 0), (1, 1), (2, 2)])
>>> point = Point(0.3, 0.7)
>>> point
POINT (0.3000000000000000 0.7000000000000000)
>>> line.interpolate(line.project(point))
POINT (0.5000000000000000 0.5000000000000000)

Wenn ich Sie richtig verstehe, ist die Funktionalität, nach der Sie fragen, in PostGIS integriert.

Um einen Punkt auf eine Linie projiziert zu bekommen, können Sie ST_Closestpoint (auf PostGIS 1.5) verwenden.

Einige Hinweise zur Verwendung finden Sie hier: http://blog.jordogskog.no/2010/02/07/how-to-use-the-new-distance-related-functions-in-postgis-part1/

Sie können beispielsweise auch den nächsten Punkt auf einem Polygon zu einem anderen Polygon suchen.

Wenn Sie die Linie zwischen den beiden nächstgelegenen Punkten auf beiden Geometrien möchten, können Sie ST_Shortestline verwenden. ST_Closestpoint ist der erste Punkt in ST_Shortestline

Die Länge von ST_Shortestline zwischen zwei Geometrien entspricht der Länge von ST_Distance zwischen den Geometrien.

Im folgenden Kommentar wird erläutert, wie meine Antwort nicht als verlässliche Lösung angesehen werden kann. Ich werde diesen ursprünglichen Beitrag hier belassen, damit andere das Problem untersuchen können.

Wenn ich die Frage verstehe, sollte dieses allgemeine Verfahren funktionieren.

So finden Sie den kürzesten Weg zwischen einem Punkt (wie durch x, y oder x, y, z definiert) und einem Polyin (wie durch eine Verbindungsmenge von x, y oder x, y, z definiert) im euklidischen Raum:

1) Suchen Sie von einem bestimmten benutzerdefinierten Punkt aus (ich nenne es pt0) den nächsten Scheitelpunkt der Polylinie (pt1). OGRinfo kann die Eckpunkte einer Polylinie abfragen, und dann können Entfernungsberechnungen mit Standardmethoden durchgeführt werden. Beispiel: Iterieren Sie über eine Distanz. Berechnen Sie wie folgt: distance_in_radians = 2 * math.asin (math.sqrt (math.pow ((math.sin ((pt0_radians-ptx_radians) / 2)), 2) + math.cos (pt0_radians) * math.cos (ptx_radians) * math.pow ((math.sin ((pt0_radians-ptx_radians) / 2)), 2)))

2) Den zugehörigen Mindestabstandswert (d1) und (pt1) speichern

3) Sehen Sie sich die beiden Segmente an, die von pt1 abstammen (in der Ogrinfo-Linienfolge sind dies die vorherigen und nachfolgenden Scheitelpunkte). Notieren Sie diese Eckpunkte (n2 und n3).

4) Erstelle für jedes Segment eine y = mx + b Formel

5) Beziehen Sie Ihren Punkt (pt0) für jede dieser beiden Formeln auf die Senkrechte

6) Entfernungen und Schnittpunkte berechnen (d2 und d3; pt2, pt3)

7) Vergleichen Sie die drei Abstände (d1, d2, d3) für den kürzesten. Ihr pt0 zum zugehörigen Knoten (pt1, pt2 oder pt3) ist die kürzeste Verbindung.

Das ist ein Strom von Bewusstseinsantwort - hoffentlich passt mein mentales Bild des Problems und der Lösung.

Hier ist ein Python-Skript für QGIS> 2.0, das aus den oben angegebenen Hinweisen und Lösungen erstellt wurde. Es funktioniert gut für eine angemessene Anzahl von Punkten und Linien. Aber ich habe es nicht mit einer großen Anzahl von Objekten versucht.

Natürlich musste es im Leerlauf oder in einer anderen weniger "pythonischen Lösung" kopiert und als "next.point.py" gespeichert werden.

Wählen Sie in der QGIS-Toolbox Skript, Werkzeuge, fügen Sie ein Skript hinzu und wählen Sie es aus.

##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
[email protected]
"""
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()

!!! WARNUNG !!! Achten Sie darauf, dass einige "seltsame" / falsch projizierte Punkte aufgrund dieses Zeilenbefehls erzeugt werden können:

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

Der counterSelecWert darin legt fest, wie viele nextNeighbor zurückgegeben werden. Tatsächlich sollte jeder Punkt so weit wie möglich von jedem Linienobjekt entfernt projiziert werden. und der gefundene Mindestabstand würde die richtige Linie und den richtigen projizierten Punkt als die nächsten Nachbarn ergeben, die wir suchen. Um die Durchlaufzeit zu verkürzen, wird der nächste Nachbarbefehl verwendet. Wenn Sie den counterSelecWert auf 1 reduzieren, wird das "erste" Objekt zurückgegeben (genauer gesagt das Begrenzungsfeld), und es ist möglicherweise nicht das richtige. Objekte mit unterschiedlicher Liniengröße müssen möglicherweise 3 oder 5 auswählen, oder es müssen noch mehr Objekte in der Nähe ausgewählt werden, um die kürzeste Entfernung zu bestimmen. Je höher der Wert, desto länger dauert es. Mit Hunderten von Punkten und Linien wird es bei 3 oder 5 nächsten Nachbarn sehr langsam, bei Tausenden kann es Fehler bei solchen Werten geben.

Abhängig von Ihren Interessen und Ihrem Anwendungsfall kann es hilfreich sein, sich mit "Map Matching-Algorithmen" zu befassen. Es gibt beispielsweise ein RoadMatcher-Projekt im OSM-Wiki: http://wiki.openstreetmap.org/wiki/Roadmatcher .