Déplacer des points sur des lignes (~ quartier)

14

J'ai deux couches vectorielles, dont l'une est une couche ponctuelle basée sur des "événements" par télédétection et la seconde est une couche linéaire issue de la recherche locale.

Dans mon cas, ce sont des tremblements de terre et des failles tectoniques, mais je suppose que l'on pourrait simplement choisir "les accidents de voiture et les routes" comme exemple général.

Donc, ce que je voudrais faire, c'est déplacer / copier les points sur le point le plus proche des lignes, tant que c'est dans une distance de tolérance (disons 1-2 km ou 0,0xx °), avec la nouvelle couche de points (+ attr déplacé o / n).

Des idées ?

Linux, QGIS 1.8

Chris Pallasch
la source
Cherchez-vous une fonction totalement automatisée pour ce faire, ou une sorte d'outil de capture pour le faire à la main serait-elle OK?
Simbamangu
J'ai posé une question similaire, j'essayais d'aligner la ligne aux points mais je n'ai jamais trouvé de solution facile. gis.stackexchange.com/questions/52232/…
GreyHippo
Qu'en est-il de la triangulation et de l'appariement des distances?
huckfinn
J'ai trouvé cette question sur une méthode qui fonctionne dans ArcGIS à l'aide de Near. Je suis allé chercher QGIS Near équivalent et j'ai trouvé ce message dans le forum où quelqu'un a suggéré GRASS v.distance. Cela m'a amené à ce tutoriel qui peut identifier une méthode. Peut-être que quelque part là-dedans, quelqu'un a déjà écrit un plugin?
Chris W

Réponses:

13

Publié un extrait de code (testé dans la console python) qui fait ce qui suit

  1. Utilisez QgsSpatialIndex pour trouver l'entité linéaire la plus proche d'un point
  2. Trouvez le point le plus proche sur cette ligne du point. J'ai utilisé un paquet galbé comme raccourci pour cela. J'ai trouvé les méthodes QGis pour cela insuffisantes (ou très probablement je ne les comprends pas correctement)
  3. Élastiques ajoutés aux emplacements d'accrochage
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)

Edit: Je viens de découvrir que la méthode @radouxju utilisant le plus procheSegmentWithContext donne les mêmes résultats en moins de lignes de code. Je me demande pourquoi ils ont trouvé ce nom de méthode étrange? aurait dû être quelque chose comme le plus prochePointOnGeometry.

Nous pouvons donc éviter les galbes et faire comme,

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
vinayan
la source
1
courir dans le cauchemar en essayant de formater ce code python..argh !!
vinayan
5

voici un pseudo-code pour commencer. J'espère que cela aide et que quelqu'un aura le temps de fournir le code complet (je n'en ai pas pour le moment)

La première chose à faire est de boucler sur le point et de sélectionner les lignes situées à l'intérieur de la distance seuil de chaque point. Cela peut être fait avec QgsSpatialIndex

Dans la première boucle, la deuxième chose à faire est de boucler sur les lignes sélectionnées et de trouver le point le plus proche sur la ligne. Cela peut être fait directement sur la base de QgsGeometry :: narrowSegmentWithContext

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

Recherche le segment de géométrie le plus proche du point donné.

Paramètres point Spécifie le point de recherche

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: retourne si le point se trouve à gauche du côté droit du segment (<0 signifie gauche,> 0 signifie droite) epsilon epsilon pour l'accrochage de segment (ajouté en 1.8)

la troisième étape (au sein de la première boucle) consisterait à mettre à jour la géométrie du point avec la géométrie du minDistPoint avec la plus petite distance

mise à jour avec du code (sur 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)
radouxju
la source