Najbliższy sąsiad między warstwą punktową a warstwową? [Zamknięte]

Zadawałem to pytanie kilka razy na przepełnieniu stosu i IRC między #qgis a #postgis, a także próbowałem go zakodować lub zaimplementować samodzielnie w postgis bez prawdziwej odpowiedzi.

Korzystając z programowania (najkorzystniej w Pythonie), chciałbym narysować linię z warstwy punktowej do jej rzutowania na najbliższą linię warstwy lub warstwy wielokąta.

Obecnie większość moich danych jest w kształcie i formacie ESRI; wolałbym jednak trzymać się z dala od rozwiązania postgis, ponieważ jestem głównie użytkownikiem shp + qgis.

Idealnym rozwiązaniem byłoby wdrożenie GDAL / OGR za pomocą Pythona lub podobnych bibliotek

• Korzystanie z bibliotek GDAL / OGR od czego zacząć? czy można podać plan rozwiązania?
• Czy mogę użyć NetworkX do wykonania analizy najbliższego sąsiada?
• Czy to rzeczywiście możliwe?

Jeśli jest to łatwiejsze, punkty mogą łączyć się z punktem końcowym segmentu zamiast punktu rzutowanego

To pytanie okazało się nieco trudniejsze, niż się spodziewałem. Istnieje wiele implementacji samej najkrótszej odległości, takich jak Shapely zapewniona odległość (z GEOS). Niewiele rozwiązań zapewnia jednak sam punkt przecięcia, ale tylko odległość.

Moja pierwsza próba zbuforowała punkt o odległość między punktem a wielokątem i szukałam skrzyżowań, ale błędy zaokrąglania uniemożliwiają udzielenie dokładnej odpowiedzi.

Oto kompletne rozwiązanie wykorzystujące Shapely, oparte na tych równaniach :

#!/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)

Dla potomnych wygląda na to, że to rozszerzenie ArcView całkiem dobrze radzi sobie z tym problemem, szkoda, że ​​jest na martwej platformie napisanej martwym językiem ...

Odpowiedź PostGIS (dla multilinestring, jeśli linestring, usuń funkcję 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

To jest trochę stare, ale szukałem dziś rozwiązania tego problemu (punkt -> linia). Najprostszym rozwiązaniem, z jakim natknąłem się na ten powiązany problem, jest:

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

Jeśli dobrze cię rozumiem, funkcjonalność, o którą prosisz, jest wbudowana w PostGIS.

Aby uzyskać punkt rzutowany na linię, możesz użyć ST_Closestpoint (na PostGIS 1.5)

Kilka wskazówek na temat korzystania z niego można znaleźć tutaj: http://blog.jordogskog.no/2010/02/07/how-to-use-the-new-distance-related-functions-in-postgis-part1/

Na przykład przydatne jest znalezienie najbliższego punktu wielokąta do innego wielokąta.

Jeśli chcesz uzyskać linię między dwoma najbliższymi punktami na obu geometriach, możesz użyć ST_Shortestline. ST_Closestpoint to pierwszy punkt w ST_Shortestline

Długość ST_Shortestline między dwiema geometriami jest taka sama jak ST_Distance między geometriami.

Zobacz komentarz poniżej, w jaki sposób moja odpowiedź nie powinna być uważana za wiarygodne rozwiązanie ... Zostawię ten oryginalny post tutaj, aby inni mogli zbadać problem.

Jeśli rozumiem pytanie, ta ogólna procedura powinna zadziałać.

Aby znaleźć najkrótszą ścieżkę między punktem (zdefiniowanym przez x, y lub x, y, z) i poliną (zdefiniowanym przez zestaw łączący x, y lub x, y, z) w przestrzeni euklidesowej:

1) Z danego punktu zdefiniowanego przez użytkownika (nazywam go pt0), znajdź najbliższy wierzchołek polilinii (pt1). OGRinfo może sondować wierzchołki polilinii, a następnie można wykonać obliczenia odległości za pomocą standardowych metod. Na przykład, iteruj w odległościach obliczeniowych, takich jak: distance_in_radians = 2 * matath.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) Zapisz powiązaną minimalną odległość (d1) i (pt1)

3) spójrz na dwa segmenty odchodzące od pt1 (w linii ogrinfo będą to wcześniejsze i następne wierzchołki). Zapisz te wierzchołki (n2 i n3).

4) utwórz formułę y = mx + b dla każdego segmentu

5) Odnieś swój punkt (pt0) do prostopadłej dla każdej z tych dwóch formuł

6) Oblicz odległość i skrzyżowania (d2 i d3; pt2, pt3)

7) Porównaj najkrótsze trzy odległości (d1, d2, d3). Twój pt0 do powiązanego węzła (pt1, pt2 lub pt3) jest najkrótszym linkiem.

To jest strumień odpowiedzi świadomości - mam nadzieję, że mój mentalny obraz problemu i rozwiązania pasuje.

Oto skrypt Pythona dla QGIS> 2.0 wykonany na podstawie podanych wyżej wskazówek i rozwiązań. Działa dobrze dla rozsądnej liczby punktów i linii. Ale nie próbowałem tego z dużą ilością przedmiotów.

Oczywiście musiał on zostać skopiowany w trybie bezczynności lub w jakimkolwiek innym „pythonicznym rozwiązaniu” i zapisać go jako „closest.point.py”.

W przyborniku QGIS przejdź do skryptu, narzędzi, dodaj skrypt, wybierz go.

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

!!! OSTRZEŻENIE !!! Zwróć uwagę, że niektóre „dziwne” / nieprawidłowe rzutowane punkty mogą być generowane z powodu tego polecenia linii:

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

counterSelecWartość w nim ustawić ile nearestneighbor są zwracane. W rzeczywistości każdy punkt powinien być rzutowany w możliwie najkrótszej odległości do każdego obiektu linii; a znaleziona minimalna odległość dałaby prawidłową linię i rzutowany punkt jako najbliżsi sąsiedzi, których szukamy. Aby skrócić czas zapętlenia, używana jest najbliższa komenda Sąsiad. Wybranie counterSelecwartości zmniejszonej do 1 zwróci „pierwszy” spełniony obiekt (a konkretnie jego ramkę ograniczającą) i może nie być właściwy. Obiekty o różnych rozmiarach linii mogą zobowiązać do wyboru 3 lub 5, a nawet najbliższe obiekty w celu ustalenia najkrótszej odległości. Im wyższa wartość, tym dłużej trwa. Z setkami punktów i linii zaczyna się bardzo powoli z 3 lub 5 najbliższymi sąsiadami, z tysiącami może powodować błędy z takimi wartościami.

W zależności od zainteresowań i przypadku użycia przydatne może być zapoznanie się z „algorytmami dopasowywania map”. Na przykład istnieje projekt RoadMatcher na wiki OSM: http://wiki.openstreetmap.org/wiki/Roadmatcher .