Mover puntos a líneas (~ vecindario)

14

Tengo dos capas de vectores, de las cuales una es una capa de puntos basada en "eventos" por teledetección y la segunda es una capa de línea de investigación local.

En mi caso, se trata de terremotos y fallas tectónicas, pero supongo que uno podría simplemente elegir "accidentes automovilísticos y carreteras" como ejemplo general.

Entonces, lo que me gustaría hacer es mover / copiar los puntos en el punto más cercano de las líneas, siempre que esté dentro de una distancia de tolerancia (digamos 1-2 km o 0.0xx °), con la nueva capa de puntos (+ attr movido y / n).

Algunas ideas ?

Linux, QGIS 1.8

Chris Pallasch
fuente
3
Habría una solución PostGIS: PostGIS: punto más cercano en una cadena
lineal
¿Está buscando una función totalmente automatizada para hacer esto, o algún tipo de herramienta de ajuste para hacerlo a mano estaría bien?
Simbamangu
Hice una pregunta similar, estaba tratando de ajustar la línea a los puntos pero nunca encontré una solución fácil. gis.stackexchange.com/questions/52232/…
GreyHippo
¿Qué pasa con la triangulación y la correspondencia de distancia?
huckfinn
Encontré esta pregunta sobre un método que funciona en ArcGIS usando Near. Fuimos a buscar QGIS Casi equivalente y encontré esta publicación en el foro donde alguien sugirió GRASS v.distance. Eso me llevó a este tutorial que puede identificar un método. ¿Quizás en algún lugar alguien ya ha escrito un complemento?
Chris W

Respuestas:

13

Publicó un fragmento de código (probado en la consola de Python) que hace lo siguiente

  1. Use QgsSpatialIndex para encontrar la entidad de línea más cercana a un punto
  2. Encuentra el punto más cercano en esta línea al punto. Utilicé el paquete bien proporcionado como atajo para esto. Encontré los métodos QGis para esto como insuficientes (o muy probablemente no los entiendo correctamente)
  3. Bandas de goma añadidas a las ubicaciones de ajuste
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)

Editar: justo ahora descubrí que el método @radouxju que utiliza el más cercanoSegmentWithContext da los mismos resultados en menos líneas de código. Me pregunto por qué se les ocurrió este extraño nombre de método. debería haber sido algo así comoarestPointOnGeometry.

Así podemos evitar bien y hacer me gusta,

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
fuente
1
toparse con una pesadilla tratando de formatear este código de Python ... ¡¡¡Argh !!
vinayan
5

Aquí hay un pseudocódigo para comenzar. Espero que esto ayude y que alguien tenga tiempo para proporcionar el código completo (no lo tengo en este momento)

Lo primero que debe hacer es recorrer el punto y seleccionar las líneas que se encuentran dentro de la distancia umbral a cada punto. Esto se puede hacer con QgsSpatialIndex

Dentro del primer bucle, lo segundo que debe hacer es recorrer las líneas seleccionadas y encontrar el punto más cercano en la línea. Esto se puede hacer directamente basado en QgsGeometry ::arestSegmentWithContext

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

Busca el segmento de geometría más cercano al punto dado.

Parámetros punto Especifica el punto de búsqueda

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: Devuelve si el punto se encuentra a la izquierda del lado derecho del segmento (<0 significa izquierda,> 0 significa derecha) epsilon epsilon para el ajuste de segmento (agregado en 1.8)

El tercer paso (dentro del primer bucle) consistiría en actualizar la geometría del punto con la geometría del minDistPoint con la menor distancia

actualizar con algún código (en 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
fuente