Bagaimana cara mengambil jaringan jalan ke jaringan heksagonal di QGIS?

Saya mencoba menggunakan QGIS 2.14 untuk mengambil jaringan jalan ke jaringan heksagonal, tapi saya mendapatkan artefak yang aneh.

Saya telah membuat hex grid dengan MMQGIS , sel kira-kira 20 x 23 m. Saya telah buffer jaringan jalan dengan 1m dan memadatkannya sehingga ada node setiap beberapa meter. Anda dapat melihat apa yang saya coba capai di bawah ini. Seperti yang Anda lihat, saya bisa membuatnya berfungsi dalam beberapa kasus: -

• biru adalah jalan padat (garis buffered)
• merah adalah versi 'hexified' - ini yang ingin saya temukan
• abu-abu adalah kisi hex

Saya kemudian menggunakan fitur Snap geometri baru untuk mengambil node ke sudut segi enam terdekat. Hasilnya menjanjikan, tetapi tampaknya ada beberapa kasus tepi di mana garis mengembang untuk mengisi segi enam (atau bagian dari itu): -

Alasan buffer adalah bahwa geometri Snap tidak memungkinkan Anda untuk mengambil lapisan yang geometrinya berbeda. Sebagai contoh, Anda tidak dapat mengambil node pada layer LINE ke poin pada layer POINT). Tampaknya menjadi yang paling menyenangkan memotret POLYGON ke POLYGON.

Saya menduga jalan melebar ketika satu sisi dari garis jalan yang disangga melompat ke satu sisi sel hex, dan sisi lain melompat ke sisi lain dari sel hex. Dalam contoh saya, jalan-jalan yang melintasi barat-timur pada sudut yang akut tampaknya yang terburuk.

Hal yang sudah saya coba, tanpa hasil: -

• buffering jaringan jalan dengan jumlah kecil, sehingga tetap poligon tetapi sangat tipis.
• memadatkan sel hex (jadi ada node di sepanjang tepi, tidak hanya di sudut)
• memvariasikan jarak gertakan maksimum (ini memiliki efek terbesar, tetapi saya tidak dapat menemukan nilai yang ideal)
• menggunakan layer LINE, bukan POLYGON

Saya menemukan bahwa jika saya mengubah menggunakan hanya lapisan LINE, itu berfungsi untuk sementara waktu, kemudian macet. Tampaknya untuk menyimpan pekerjaannya - beberapa baris telah diproses sebagian.

Adakah yang tahu cara lain untuk mengambil titik pada garis ke titik terdekat pada garis / lapisan poligon lain, idealnya tanpa perlu menggunakan postgres / postgis (walaupun solusi dengan postgis akan diterima juga)?

EDIT

Bagi siapa saja yang ingin mencoba, saya telah meletakkan proyek QGIS pemula di sini di Dropbox . Ini termasuk Hex Grid dan layer garis Densified. (Jaringan jalan dari OSM, jadi dapat diunduh menggunakan QuickOSM misalnya jika Anda perlu mendapatkan dokumen asli untuk mengurangi jalan).

Perhatikan bahwa ini dalam OSGB (epsg: 27700) yang merupakan UTM lokal untuk Inggris, dengan satuan dalam meter.

Solusi saya melibatkan skrip PyQGIS yang lebih cepat dan lebih efektif daripada alur kerja yang melibatkan gertakan (saya mencobanya juga). Dengan menggunakan algoritma saya, saya telah memperoleh hasil ini:

Anda dapat menjalankan cuplikan kode berikut secara berurutan dari dalam QGIS (di konsol Python QGIS). Pada akhirnya Anda mendapatkan lapisan memori dengan rute yang diambil dimuat ke QGIS.

Satu-satunya prasyarat adalah membuat Shapefile jalan multipart (gunakan Processing->Singleparts to multipart, saya menggunakan bidang fictitiuossebagai Unique ID fieldparameter). Ini akan memberi kita roads_multipart.shpfile dengan fitur tunggal.

Berikut adalah algoritme yang dijelaskan:

1. Dapatkan sisi segi enam terdekat tempat rute melintas. Untuk setiap segi enam kami membuat 6 segitiga antara setiap pasangan simpul tetangga dan centroid yang sesuai. Jika ada jalan yang memotong segitiga, ruas yang dibagikan oleh segi enam dan segitiga ditambahkan ke rute terakhir yang terpotong. Ini adalah bagian yang lebih berat dari keseluruhan algoritma, dibutuhkan 35 detik untuk berjalan di mesin saya. Di dua baris pertama ada 2 jalur Shapefile, Anda harus menyesuaikan mereka agar sesuai dengan jalur file Anda sendiri.

   hexgrid = QgsVectorLayer("/docs/borrar/hex_grid_question/layers/normal-hexgrid.shp", "hexgrid", "ogr")
   roads = QgsVectorLayer("/docs/borrar/hex_grid_question/layers/roads_multipart.shp", "roads", "ogr")  # Must be multipart!
   
   roadFeat = roads.getFeatures().next() # We just have 1 geometry
   road = roadFeat.geometry() 
   indicesHexSides = ((0,1), (1,2), (2,3), (3,4), (4,5), (5,0))
   
   epsilon = 0.01
   # Function to compare whether 2 segments are equal (even if inverted)
   def isSegmentAlreadySaved(v1, v2):
       for segment in listSegments:        
           p1 = QgsPoint(segment[0][0], segment[0][1])
           p2 = QgsPoint(segment[1][0], segment[1][1])
           if v1.compare(p1, epsilon) and v2.compare(p2, epsilon) \
               or v1.compare(p2, epsilon) and v2.compare(p1, epsilon):
               return True
       return False
   
   # Let's find the nearest sides of hexagons where routes cross
   listSegments = []
   for hexFeat in hexgrid.getFeatures():
       hex = hexFeat.geometry()
       if hex.intersects( road ):
           for side in indicesHexSides:
               triangle = QgsGeometry.fromPolyline([hex.centroid().asPoint(), hex.vertexAt(side[0]), hex.vertexAt(side[1])])
               if triangle.intersects( road ):
                   # Only append new lines, we don't want duplicates!!!
                   if not isSegmentAlreadySaved(hex.vertexAt(side[0]), hex.vertexAt(side[1])): 
                       listSegments.append( [[hex.vertexAt(side[0]).x(), hex.vertexAt(side[0]).y()], [hex.vertexAt(side[1]).x(),hex.vertexAt(side[1]).y()]] )  

2. Singkirkan segmen yang terputus (atau 'terbuka') dengan menggunakan daftar Python, tuple, dan kamus . Pada titik ini, ada beberapa segmen yang terputus yang tersisa, yaitu, segmen yang memiliki satu simpul terputus tetapi yang lainnya terhubung ke setidaknya 2 segmen lainnya (lihat segmen merah pada gambar berikutnya). Kita harus menyingkirkan mereka.

   

   # Let's remove disconnected/open segments
   lstVertices = [tuple(point) for segment in listSegments for point in segment]
   dictConnectionsPerVertex = dict((tuple(x),lstVertices.count(x)-1) for x in set(lstVertices))
   
   # A vertex is not connected and the other one is connected to 2 segments
   def segmentIsOpen(segment):
       return dictConnectionsPerVertex[tuple(segment[0])] == 0 and dictConnectionsPerVertex[tuple(segment[1])] >= 2 \
           or dictConnectionsPerVertex[tuple(segment[1])] == 0 and dictConnectionsPerVertex[tuple(segment[0])] >= 2
   
   # Remove open segments
   segmentsToDelete = [segment for segment in listSegments if segmentIsOpen(segment)]        
   for toBeDeleted in segmentsToDelete:
       listSegments.remove( toBeDeleted )

3. Sekarang kita dapat membuat layer vektor dari daftar koordinat dan memuatnya ke peta QGIS :

   # Create a memory layer and load it to QGIS map canvas
   vl = QgsVectorLayer("LineString", "Snapped Routes", "memory")
   pr = vl.dataProvider()
   features = []
   for segment in listSegments:
       fet = QgsFeature()
       fet.setGeometry( QgsGeometry.fromPolyline( [QgsPoint(segment[0][0], segment[0][1]), QgsPoint(segment[1][0], segment[1][1])] ) )
       features.append(fet)
   
   pr.addFeatures( features )
   vl.updateExtents()
   QgsMapLayerRegistry.instance().addMapLayer(vl)

Bagian lain dari hasilnya:

Jika Anda memerlukan atribut dalam rute yang diambil, kami dapat menggunakan Indeks Spasial untuk mengevaluasi persimpangan dengan cepat (seperti di /gis//a/130440/4972 ), tapi itu cerita lain.

Semoga ini membantu!

Saya melakukannya di ArcGIS, pasti bisa diimplementasikan menggunakan QGIS atau hanya python dengan paket yang mampu membaca geometri. Pastikan jalan mewakili jaringan, yaitu saling berpotongan di ujungnya saja. Anda berhadapan dengan OSM, saya kira memang begitu.

• Konversikan poligon kedekatan menjadi garis dan planarise, sehingga menjadi jaringan geometris juga.
• Tempatkan poin di ujung mereka - Poin Voronoi:
• Tempatkan poin di jalan pada interval reguler 5 m, pastikan jalan jaringan memiliki nama unik yang bagus:

• Untuk setiap Road Point, temukan koordinat Titik Voronoi terdekat:
• Buat "Jalan" dengan menghubungkan titik-titik terdekat dalam urutan yang sama:

Jika Anda tidak ingin melihat ini:

Jangan mencoba menggunakan poin rantai pada Garis Voronoi. Saya khawatir itu hanya akan memperburuknya. Jadi, satu-satunya pilihan Anda adalah membuat jaringan dari jalur Voronoi dan menemukan rute antara titik akhir jalan, itu juga bukan masalah besar

Saya menyadari bahwa Anda meminta metode QGIS, tetapi berikan jawaban yang lebih jelas:

roads = 'clipped roads' # roads layer
hexgrid = 'normal-hexgrid' # hex grid layer
sr = arcpy.Describe('roads').spatialReference # spatial reference
outlines = [] # final output lines
points = [] # participating grid vertices
vert_dict = {} # vertex dictionary
hex_dict = {} # grid dictionary
with arcpy.da.SearchCursor(roads,["SHAPE@","OID@"], spatial_reference=sr) as r_cursor: # loop through roads
    for r_row in r_cursor:
        with arcpy.da.SearchCursor(hexgrid,["SHAPE@","OID@"], spatial_reference=sr) as h_cursor: # loop through hex grid
            for h_row in h_cursor:
                if not r_row[0].disjoint(h_row[0]): # check if the shapes overlap
                    hex_verts = []
                    for part in h_row[0]:
                        for pnt in part:
                            hex_verts.append(pnt) # add grid vertices to list
                    int_pts = r_row[0].intersect(h_row[0],1) # find all intersection points between road and grid
                    hex_bnd = h_row[0].boundary() # convert grid to line
                    hex_dict[h_row[1]] = hex_bnd # add grid geometry to dictionary
                    for int_pt in int_pts: # loop through intersection points
                        near_dist = 1000 # arbitrary large number
                        int_pt = arcpy.PointGeometry(int_pt,sr)
                        for hex_vert in hex_verts: # loop through hex vertices
                            if int_pt.distanceTo(hex_vert) < near_dist: # find shortest distance between intersection point and grid vertex
                                near_vert = hex_vert # remember geometry
                                near_dist = int_pt.distanceTo(hex_vert) # remember distance
                        vert_dict.setdefault(h_row[1],[]).append(arcpy.PointGeometry(near_vert,sr)) # store geometry in dictionary
                        points.append(arcpy.PointGeometry(near_vert,sr)) # add to points list
for k,v in vert_dict.iteritems(): # loop through participating vertices
    if len(v) < 2: # skip if there was only one vertex
        continue
    hex = hex_dict[k] # get hex grid geometry
    best_path = hex # longest line possible is hex grid boundary
    for part in hex:
        for int_vert in v: # loop through participating vertices
            for i,pnt in enumerate(part): # loop through hex grid vertices
                if pnt.equals(int_vert): # find vertex index on hex grid corresponding to current point
                    start_i = i
                    if start_i == 6:
                        start_i = 0
                    for dir in [[0,6,1],[5,-1,-1]]: # going to loop once clockwise, once counter-clockwise
                        past_pts = 0 # keep track of number of passed participating vertices
                        cur_line_arr = arcpy.Array() # polyline coordinate holder
                        cur_line_arr.add(part[start_i]) # add starting vertex to growing polyline
                        for j in range(dir[0],dir[1],dir[2]): # loop through hex grid vertices
                            if past_pts < len(v): # only make polyline until all participating vertices have been visited
                                if dir[2] == 1: # hex grid vertex index bookkeeping
                                    if start_i + j < 6:
                                        index = start_i + j
                                    else:
                                        index = (start_i - 6) + j
                                else:
                                    index = j - (5 - start_i)
                                    if index < 0:
                                        index += 6
                                cur_line_arr.add(part[index]) # add current vertex to growing polyline
                                for cur_pnt in v:
                                    if part[index].equals(cur_pnt): # check if the current vertex is a participating vertex
                                        past_pts += 1 # add to counter
                        if cur_line_arr.count > 1:
                            cur_line = arcpy.Polyline(cur_line_arr,sr)
                            if cur_line.length < best_path.length: # see if current polyline is shorter than any previous candidate
                                best_path = cur_line # if so, store polyline
    outlines.append(best_path) # add best polyline to list
arcpy.CopyFeatures_management(outlines, r'in_memory\outlines') # write list
arcpy.CopyFeatures_management(points, r'in_memory\mypoints') # write points, if you want

Catatan:

• Script ini berisi banyak loop dalam loop dan kursor bersarang. Pasti ada ruang untuk optimasi. Saya menelusuri kumpulan data Anda dalam beberapa menit, tetapi lebih banyak fitur akan menambah masalah.

Jika Anda membagi garis jalan menjadi segmen di mana setiap segmen benar-benar terkandung oleh segi enam, keputusan Anda di mana segmen garis segi enam yang akan digunakan adalah apakah jarak dari pusat massa segmen jalan split ke setiap titik tengah sisi segi enam kurang dari setengah diameter segi enam (atau kurang dari jari-jari lingkaran yang pas di dalam segi enam).

Jadi, jika Anda (satu segmen pada satu waktu) memilih segmen garis hexagon (di mana setiap segmen adalah sisi hexagon) yang berada dalam jarak jari-jari hexagon, Anda dapat menyalin geometri garis tersebut dan menggabungkannya pada pengidentifikasi unik apa pun yang Anda gunakan untuk dataset jalan Anda.

Jika Anda mengalami masalah saat menggabungkan pengidentifikasi unik, Anda dapat menerapkan buffer dan memilih berdasarkan lokasi hanya pada segmen tersebut untuk menerapkan atribut dataset jalan Anda; dengan begitu Anda tidak perlu khawatir membuat kecocokan salah dengan buffer yang terlalu besar.

Masalah dengan alat sekejap adalah bahwa alat itu terkunci secara membabi buta; sulit untuk menemukan toleransi yang sempurna untuk digunakan. Dengan metodologi ini, Anda akan mengidentifikasi dengan benar segmen garis heksagon mana yang akan digunakan, lalu mengganti geometri data jalan Anda (atau memasukkan geometri ke dalam dataset yang berbeda).

Juga, jika Anda masih memiliki masalah dengan segmen garis yang melompat dari satu sisi segi enam ke sisi yang lain, Anda dapat membagi garis menjadi segmen dengan simpul, menghitung panjang setiap garis, lalu menghapus segmen garis yang lebih besar dari panjang rata-rata satu sisi segi enam.

Kakap geometri di qgis 3.0 telah dikerjakan ulang dan sekarang memungkinkan pemotretan di antara berbagai jenis geometri. Ini juga memiliki banyak perbaikan. Anda dapat mencoba versi "snapshot harian" untuk mendapatkan akses ke kakap yang ditingkatkan sebelum 3.0 dirilis secara resmi.