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

37

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

dassouki
źródło
czy linia może być ograniczona do prostopadłości do odcinka linii?
WolfOdrade
@wolfOdrade - Ogólnie rzecz biorąc, to nie ma znaczenia.
dassouki

Odpowiedzi:

33

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

scw
źródło
1
Zastanawiam się, czy istnieje sposób na indeksowanie punktów wielokąta, aby uniknąć jawnego wyliczenia ...
mlt
@mlt nie wiem dokładnie, o czym myślisz, ale istnieją pewne podejścia, które mogą pomóc w zależności od geometrii. Mógłby wykonać kilka podstawowych rzutów promieniami, aby określić odpowiednie najbliższe segmenty, gdyby wydajność była problemem. W tym duchu przeniesienie tego do C lub Pyrex poprawiłoby sytuację.
scw
Mam na myśli, pairsże z algorytmem jest O (n) czy coś. Rozwiązanie @eprand może być może zostać zmodyfikowane, aby korzystać z KNN, ale jak dotąd udało mi się żyć bez PostGIS ...
mlt
Nie mogę już edytować mojego poprzedniego komentarza :( Być może rozwiązanie Nicklasa Avéna z ST_Closestpoint i ST_Shortestline jest najszybsze, jeśli PostGIS jest opcją.
mlt
Racja, możesz użyć algorytmu KNN bezpośrednio w Pythonie . Nie wierzę, że ST_Shortestline używa KNN, po prostu iteruje również w oparciu o moje czytanie postgis.refractions.net/documentation/postgis-doxygen/d1/dbf/…
scw
8

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
eprand
źródło
4

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)
zupa alfabetyczna
źródło
4

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.

Nicklas Avén
źródło
3

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.

glennon
źródło
To w ogóle nie zadziała. Np. Punkt = (1,1), linia = ((0,2), (0,3), (3,0), (2,0)). Jeśli naszkicujesz to, zobaczysz, że „najbliższe” wierzchołki na linii nie sąsiadują z odcinkiem, który przechodzi najbliżej punktu ... Myślę, że jedynym sposobem na poradzenie sobie z tym jest sprawdzenie każdego segmentu (ewentualnie używając pól ograniczających, aby uniknąć zoptymalizuj to trochę). HTH.
Tom
3

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

!!! 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.

Jean-Christophe Baudin
źródło
3

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 .

podmrok
źródło
To jest na potrzeby podróży i prognozowania. Zazwyczaj dzielimy obszary na strefy analizy ruchu (wielokąty) i ustalamy centroid wielokąta jako „obojętny” inicjator całego ruchu w tej strefie. Następnie narysujemy x lub y „obojętne połączenie drogowe” z tego punktu do najbliższych dróg i rozdzielamy ruch z tej strefy równomiernie na te obojętne połączenia i na rzeczywistą warstwę drogową
dassouki
Ach, więc twoim celem jest zautomatyzowanie tworzenia tego „fałszywego łącza drogowego”?
podmroku
rzeczywiście :) lub fikcyjne linki
dassouki