Memindahkan titik ke garis (~ lingkungan)

14

Saya memiliki dua lapisan vektor, yang satu adalah lapisan titik berdasarkan "peristiwa" dengan penginderaan jauh dan yang kedua adalah lapisan garis dari penelitian lokal.

Dalam kasus saya ini adalah gempa bumi dan kesalahan tektonik, tapi saya kira kita bisa memilih "kecelakaan mobil dan jalan" sebagai contoh umum.

Jadi yang ingin saya lakukan adalah memindahkan / menyalin titik ke titik terdekat dari garis, selama itu dalam jarak toleransi (katakanlah 1-2km atau 0,0xx °), dengan layer titik baru (+ attr dipindahkan y / n).

Ada ide?

Linux, QGIS 1.8

Chris Pallasch
sumber
3
Akan ada solusi PostGIS: PostGIS: titik terdekat pada linestring ke titik tertentu
underdark
Apakah Anda mencari fungsi yang sepenuhnya otomatis untuk melakukan ini, atau apakah alat gertakan untuk melakukannya dengan tangan boleh?
Simbamangu
Saya mengajukan pertanyaan yang sama, saya mencoba untuk mengubah garis ke poin tetapi tidak pernah menemukan solusi yang mudah. gis.stackexchange.com/questions/52232/...
GreyHippo
Bagaimana dengan triangulasi dan pencocokan jarak?
huckfinn
Saya menemukan pertanyaan ini tentang metode yang bekerja di ArcGIS menggunakan Near. Pergi mencari QGIS Dekat setara dan menemukan posting forum ini di mana seseorang menyarankan GRASS v.distance. Itu mengarahkan saya ke tutorial ini yang dapat mengidentifikasi metode. Mungkin di suatu tempat di sana seseorang telah menulis sebuah plugin sekarang?
Chris W

Jawaban:

13

Diposting potongan kode (diuji dalam konsol python) yang melakukan hal di bawah ini

  1. Gunakan QgsSpatialIndex untuk menemukan fitur garis terdekat ke suatu titik
  2. Temukan titik terdekat pada baris ini ke titik tersebut. Saya menggunakan paket rupawan sebagai jalan pintas untuk ini. Saya menemukan metode QGis untuk ini tidak memadai (atau kemungkinan besar saya tidak memahaminya dengan benar)
  3. Menambahkan rubberbands ke lokasi 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)

Sunting: Baru saja menemukan bahwa metode @radouxju menggunakan terdekatSegmentWithContext memberikan hasil yang sama dalam lebih sedikit baris kode. Saya bertanya-tanya mengapa mereka datang dengan nama metode yang aneh ini? seharusnya sesuatu yang terdekat denganPointOnGeometry.

Jadi kita bisa menghindari bentuk dan melakukan suka,

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
sumber
1
berlari ke dalam mimpi buruk mencoba memformat kode python ini..argh !!
vinayan
5

di sini adalah pseudo-code untuk memulai. Saya harap ini membantu dan seseorang akan punya waktu untuk memberikan kode lengkap (saya tidak punya saat ini)

Hal pertama yang harus dilakukan adalah mengulang pada titik dan memilih garis yang berada dalam jarak ambang batas untuk setiap titik. Ini dapat dilakukan dengan QgsSpatialIndex

Dalam loop pertama, hal kedua yang harus dilakukan adalah loop pada garis yang dipilih dan menemukan titik terdekat pada garis. Ini dapat dilakukan langsung berdasarkan QgsGeometry :: terdekatSegmentWithContext

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

Mencari segmen geometri terdekat ke titik tertentu.

Titik parameter Menentukan titik untuk pencarian

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 kiriOf Keluar: Mengembalikan jika titik terletak di sebelah kiri sisi kanan segmen (<0 berarti kiri,> 0 berarti kanan) epsilon epsilon untuk gertakan segmen (ditambahkan 1,8)

langkah ketiga (dalam loop pertama) akan terdiri dalam memperbarui geometri titik dengan geometri minDistPoint dengan jarak terkecil

perbarui dengan beberapa kode (pada 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
sumber