Mover puntos a líneas (~ vecindario)

14

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

En mi caso, estos son 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 un distancia de tolerancia (por ejemplo, 1-2 km o 0.0xx °), con la nueva capa de puntos (+ attr movido y / n).

¿Alguna idea?

Linux, QGIS 1.8

    
pregunta Chris Pallasch 03.01.2013 - 17:36

2 respuestas

12

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. Encuentre el punto más cercano en esta línea al punto. Utilicé el paquete bien formado como atajo para esto He encontrado los métodos QGis para esto como insuficiente (o muy probablemente no los entiendo bien)
  3. Se agregaron bandas de goma 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: Acabo de encontrar que el método @radouxju utilizando closestSegmentWithContext da los mismos resultados en menos líneas de código. Me pregunto por qué se les ocurrió este nombre de método extraño? debería haber sido algo así como closestPointOnGeometry.

Por lo tanto, podemos evitar bien y bien hacer,

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
    
respondido por el vinayan 29.04.2014 - 11:04
4

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

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

Dentro del primer bucle, lo segundo que hay que hacer es realizar un bucle en las líneas seleccionadas y encontrar el punto más cercano en la línea. Esto se puede hacer directamente en función de QgsGeometry :: closestSegmentWithContext

  

double QgsGeometry :: closestSegmentWithContext (const QgsPoint &     punto, QgsPoint & minDistPoint, int & afterVertex, doble *     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 para la 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 o al lado derecho del segmento (< 0 significa izquierda, > 0 significa derecha)       épsilon épsilon para segmentar segmentos (añadido 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 distancia más pequeña

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)
    
respondido por el radouxju 27.04.2014 - 22:00

Lea otras preguntas en las etiquetas