Movendo pontos para linhas (~ vizinhança)

14

Eu tenho duas camadas vetoriais, das quais uma é uma camada pontual baseada em "eventos" por sensoriamento remoto e a segunda é uma camada linear de pesquisa local.

No meu caso, são terremotos e falhas tectônicas, mas acho que alguém poderia simplesmente escolher "acidentes de carro e estradas" como um exemplo geral.

Então, o que eu gostaria de fazer é mover / copiar os pontos para o ponto mais próximo das linhas, desde que a uma distância de tolerância (digamos 1-2km ou 0,0xx °), com a nova camada de pontos (+ attr movido s / n).

Alguma ideia ?

Linux, QGIS 1.8

Chris Pallasch
fonte
Você está procurando uma função totalmente automatizada para fazer isso, ou algum tipo de ferramenta de encaixe para fazê-lo manualmente seria bom?
Simbamangu
Fiz uma pergunta semelhante, estava tentando ajustar a linha aos pontos, mas nunca encontrei uma solução fácil. gis.stackexchange.com/questions/52232/...
GreyHippo
E quanto à triangulação e à distância?
precisa saber é o seguinte
Eu encontrei esta pergunta sobre um método que funciona no ArcGIS usando o Near. Foi procurar o equivalente próximo ao QGIS e encontrou esta postagem no fórum em que alguém sugeriu o GRASS v.distance. Isso me levou a este tutorial que pode identificar um método. Talvez em algum lugar lá alguém tenha escrito um plugin até agora?
Chris W

Respostas:

13

Publicou um trecho de código (testado no console python) que faz o abaixo

  1. Use QgsSpatialIndex para encontrar o recurso de linha mais próximo de um ponto
  2. Encontre o ponto mais próximo nesta linha para o ponto. Eu usei o pacote bem torneado como um atalho para isso. Eu achei os métodos QGis para isso como insuficientes (ou provavelmente eu não os entendo corretamente)
  3. Bandas de borracha adicionadas aos locais de snap
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: Acabei de descobrir que o método @radouxju usando o closestSegmentWithContext fornece os mesmos resultados em menos linhas de código. Eu me pergunto por que eles inventaram esse nome estranho de método? deveria ter sido algo como o mais próximoPointOnGeometry.

Para que possamos evitar formas e fazer algo assim,

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
fonte
1
correndo em pesadelo tentando formatar esse código python..argh !!
Vinayan
5

Aqui está um pseudo-código para começar. Espero que isso ajude e que alguém tenha tempo para fornecer o código completo (não tenho no momento)

A primeira coisa a fazer é fazer um loop no ponto e selecionar as linhas que estão localizadas dentro da distância limite de cada ponto. Isso pode ser feito com QgsSpatialIndex

Dentro do primeiro loop, a segunda coisa a fazer é fazer um loop nas linhas selecionadas e encontrar o ponto mais próximo da linha. Isso pode ser feito diretamente com base em QgsGeometry :: closestSegmentWithContext

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

Procura o segmento mais próximo da geometria para o ponto especificado.

Parâmetros point Especifica o ponto para pesquisa

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: Retorna se o ponto estiver no lado esquerdo do lado direito do segmento (<0 significa esquerda,> 0 significa direita) epsilon epsilon para ajuste de segmento (adicionado em 1.8)

o terceiro passo (dentro do primeiro loop) consistiria em atualizar a geometria do ponto com a geometria do minDistPoint com a menor distância

atualizar com algum código (no 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
fonte