Przenoszenie punktów na linie (~ sąsiedztwo)

Mam dwie warstwy wektorowe, z których jedna jest warstwą punktową opartą na „zdarzeniach” za pomocą teledetekcji, a druga to warstwa liniowa z lokalnych badań.

W moim przypadku są to trzęsienia ziemi i wady tektoniczne, ale wydaje mi się, że można po prostu wybrać „wypadki samochodowe i drogi” jako ogólny przykład.

Chciałbym więc przenieść / skopiować punkty na najbliższy punkt linii, o ile mieści się w odległości tolerancji (powiedzmy 1-2 km lub 0.0xx °), z nową warstwą punktową (przesunięcie + attr t / n).

Jakieś pomysły ?

Linux, QGIS 1.8

Wysłałem fragment kodu (przetestowany w konsoli Pythona), który spełnia poniższe wymagania

1. Użyj QgsSpatialIndex, aby znaleźć element linii najbliżej punktu
2. Znajdź najbliższy punkt na tej linii do tego punktu. Użyłem do tego zgrabnego pakietu. Znalazłem w tym celu metody QGis jako niewystarczające (lub najprawdopodobniej nie rozumiem ich poprawnie)
3. Dodano gumki do lokalizacji zatrzasków

from shapely.wkt import *
from shapely.geometry import *
from qgis.gui import *
from PyQt4.QtCore import Qt
lineLayer = iface.mapCanvas().layer(0)
pointLayer =  iface.mapCanvas().layer(1)
canvas =  iface.mapCanvas()
spIndex = QgsSpatialIndex() #create spatial index object
lineIter =  lineLayer.getFeatures()
for lineFeature in lineIter:
    spIndex.insertFeature(lineFeature)        
pointIter =  pointLayer.getFeatures()
for feature in pointIter:
    ptGeom = feature.geometry()
    pt = feature.geometry().asPoint()
    nearestIds = spIndex.nearestNeighbor(pt,1) # we need only one neighbour
    featureId = nearestIds[0]
    nearestIterator = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId))
    nearFeature = QgsFeature()
    nearestIterator.nextFeature(nearFeature)
    shplyLineString = shapely.wkt.loads(nearFeature.geometry().exportToWkt())
    shplyPoint = shapely.wkt.loads(ptGeom.exportToWkt())
    #nearest distance from point to line
    dist = shplyLineString.distance(shplyPoint)
    print dist
    #the point on the road where the point should snap
    shplySnapPoint = shplyLineString.interpolate(shplyLineString.project(shplyPoint))
    #add rubber bands to the new points
    snapGeometry = QgsGeometry.fromWkt(shapely.wkt.dumps(shplySnapPoint))
    r = QgsRubberBand(canvas,QGis.Point)
    r.setColor(Qt.red)
    r.setToGeometry(snapGeometry,pointLayer)

Edycja: Właśnie odkryłem, że metoda @radouxju przy użyciu closestSegmentWithContext daje te same wyniki w mniejszej liczbie wierszy kodu. Zastanawiam się, dlaczego wymyślili tę dziwną nazwę metody? powinien być czymś w rodzaju najbliższejPointOnGeometry.

Abyśmy mogli uniknąć zgrabnych i robić to,

nearFeature = QgsFeature()
nearestIterator.nextFeature(nearFeature)   

closeSegResult = nearFeature.geometry().closestSegmentWithContext(ptGeom.asPoint())
closePoint = closeSegResult[1]
snapGeometry = QgsGeometry.fromPoint(QgsPoint(closePoint[0],closePoint[1])) 

p1 = ptGeom.asPoint()
p2 = snapGeometry.asPoint()

dist = math.hypot(p2.x() - p1.x(), p2.y() - p1.y())
print dist

tutaj jest pseudo-kod na początek. Mam nadzieję, że to pomoże i że ktoś będzie miał czas na dostarczenie pełnego kodu (obecnie nie mam)

pierwszą rzeczą do zrobienia jest zapętlenie punktu i wybranie linii, które znajdują się w odległości progowej do każdego punktu. Można to zrobić za pomocą QgsSpatialIndex

W ramach pierwszej pętli drugą rzeczą jest zapętlenie wybranych linii i znalezienie najbliższego punktu na linii. Można to zrobić bezpośrednio w oparciu o QgsGeometry :: closestSegmentWithContext

double QgsGeometry :: closestSegmentWithContext (const QgsPoint & point, QgsPoint & minDistPoint, int & afterVertex, double * leftOf = 0, double epsilon = DEFAULT_SEGMENT_EPSILON)

Wyszukuje najbliższy segment geometrii do danego punktu.

Parametry point Określa punkt wyszukiwania

minDistPoint  Receives the nearest point on the segment

afterVertex   Receives index of the vertex after the closest segment. The vertex before the closest segment is always afterVertex -

1 leftOf Out: Zwraca, jeśli punkt leży po lewej stronie prawej części segmentu (<0 oznacza lewą stronę,> 0 oznacza prawą) epsilon epsilon do przyciągania segmentu (dodany w 1.8)

trzeci krok (w ramach pierwszej pętli) polegałby na zaktualizowaniu geometrii punktu o geometrię minDistPoint o najmniejszej odległości

aktualizacja z pewnym kodem (na QGIS3)

pointlayer = QgsProject.instance().mapLayersByName('point')[0] #iface.mapCanvas().layer(0)
lineLayer = QgsProject.instance().mapLayersByName('lines')[0] # iface.mapCanvas().layer(1)

epsg = pointlayer.crs().postgisSrid()
uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=distance:double(20,2)&field=left:integer&index=yes"
snapped = QgsVectorLayer(uri,'snapped', 'memory')

prov = snapped.dataProvider()

testIndex = QgsSpatialIndex(lineLayer)
i=0

feats=[]

for p in pointlayer.getFeatures():
    i+=1
    mindist = 10000.
    near_ids = testIndex.nearestNeighbor(p.geometry().asPoint(),4) #nearest neighbor works with bounding boxes, so I need to take more than one closest results and further check all of them. 
    features = lineLayer.getFeatures(QgsFeatureRequest().setFilterFids(near_ids))
    for tline in features:
        closeSegResult = tline.geometry().closestSegmentWithContext(p.geometry().asPoint())
        if mindist > closeSegResult[0]:
            closePoint = closeSegResult[1]
            mindist = closeSegResult[0]
            side = closeSegResult[3]
    feat = QgsFeature()
    feat.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(closePoint[0],closePoint[1])))
    feat.setAttributes([i,mindist,side])
    feats.append(feat)

prov.addFeatures(feats)
QgsProject.instance().addMapLayer(snapped)