Memotong garis untuk mendapatkan persimpangan menggunakan Python dengan QGIS?

10

Saya memiliki satu set garis yang mewakili jalur bus. Beberapa jalur tumpang tindih dan mengambil jalan yang sama.

masukkan deskripsi gambar di sini

Saya dapat mengekstrak node. masukkan deskripsi gambar di sini

Namun saya tertarik mengekstraksi penyeberangan seperti ini: masukkan deskripsi gambar di sini

Bagaimana saya bisa melakukan ini? Saya mencari cara dengan QGIS atau Python.

Saya mencoba metode persimpangan dari GDAL Python tetapi ini pada dasarnya mengembalikan saya hanya simpul.

The Jalur Persimpangan Metode dari QGIS kembali saya penyeberangan jika dua garis silang. Namun dalam kasus bahwa dua jalur bus pergi jauh dari rute mereka di jalan yang sama, itu tidak memberi saya titik di mana mereka bergabung.

ustroetz
sumber
Sudahkah Anda mencoba alat persimpangan garis di QGIS: Alat Analisis Vektor> Garis Interseksi ... Ini tidak akan memberi Anda titik akhir dan titik mulai dari sebuah garis, tetapi semua persimpangan.
Jakob
Ya saya menulis tentang ini dalam pertanyaan.
ustroetz
Saya tidak jelas tentang apa yang Anda tanyakan, sebagian karena semua garis dilambangkan dengan cara yang sama dalam gambar Anda - saya tidak dapat membedakan rute yang berbeda untuk memahami node apa yang Anda lihat atau mengapa ada begitu banyak di gambar kedua. Apakah rute itu bertepatan di jalan? Apakah mereka semua segmen garis dua titik, atau polyline kontinu? Saya perhatikan perilaku yang Anda gambarkan sama dengan alat Intersect ArcGIS - garis / garis dengan keluaran garis memberi Anda tumpang tindih, tetapi garis / garis dengan titik keluaran hanya memberikan penyilangan.
Chris W
Berdasarkan itu, untuk mendapatkan apa yang saya pikir Anda inginkan, Anda mungkin harus menggunakan kedua metode ini. Dapatkan penyeberangan (garis / garis = titik) dan kemudian tumpang tindih (garis / garis = garis) dan ekstrak mulai / akhir node untuk garis-garis yang tumpang tindih. Itu harus semua titik / node yang Anda cari.
Chris W

Jawaban:

20

Simpul:

Anda menginginkan dua hal, titik akhir dari polyline (tanpa titik tengah) dan titik persimpangan. Ada masalah tambahan, beberapa titik akhir polyline juga merupakan titik persimpangan:

masukkan deskripsi gambar di sini

Solusinya adalah dengan menggunakan Python dan modul Shapely dan Fiona

1) Baca shapefile:

from shapely.geometry import Point, shape
import fiona
lines = [shape(line['geometry']) for line in fiona.open("your_shapefile.shp")]

2) Temukan Poin akhir dari garis ( bagaimana cara mendapatkan titik akhir dari polyline? ):

endpts = [(Point(list(line.coords)[0]), Point(list(line.coords)[-1])) for line  in lines]
# flatten the resulting list to a simple list of points
endpts= [pt for sublist in endpts  for pt in sublist] 

masukkan deskripsi gambar di sini

3) Hitung persimpangan (iterasi melalui pasangan geometri di lapisan dengan modul itertools ). Hasil dari beberapa persimpangan adalah MultiPoints dan kami ingin daftar poin:

import itertools
inters = []
for line1,line2 in  itertools.combinations(lines, 2):
  if  line1.intersects(line2):
    inter = line1.intersection(line2)
    if "Point" == inter.type:
        inters.append(inter)
    elif "MultiPoint" == inter.type:
        inters.extend([pt for pt in inter])
    elif "MultiLineString" == inter.type:
        multiLine = [line for line in inter]
        first_coords = multiLine[0].coords[0]
        last_coords = multiLine[len(multiLine)-1].coords[1]
        inters.append(Point(first_coords[0], first_coords[1]))
        inters.append(Point(last_coords[0], last_coords[1]))
    elif "GeometryCollection" == inter.type:
        for geom in inter:
            if "Point" == geom.type:
                inters.append(geom)
            elif "MultiPoint" == geom.type:
                inters.extend([pt for pt in geom])
            elif "MultiLineString" == geom.type:
                multiLine = [line for line in geom]
                first_coords = multiLine[0].coords[0]
                last_coords = multiLine[len(multiLine)-1].coords[1]
                inters.append(Point(first_coords[0], first_coords[1]))
                inters.append(Point(last_coords[0], last_coords[1]))

masukkan deskripsi gambar di sini

4) Hilangkan duplikat antara titik akhir dan titik persimpangan (seperti yang Anda lihat pada gambar)

result = endpts.extend([pt for pt in inters  if pt not in endpts])

5) Simpan shapefile yang dihasilkan

from shapely.geometry import mapping
# schema of the shapefile
schema = {'geometry': 'Point','properties': {'test': 'int'}}
# creation of the shapefile
with fiona.open('result.shp','w','ESRI Shapefile', schema) as output:
    for i, pt in enumerate(result):
        output.write({'geometry':mapping(pt), 'properties':{'test':i}})

Hasil akhir:

masukkan deskripsi gambar di sini

Segmen garis

Jika Anda ingin juga segmen antara node, Anda perlu "planarize" ( Grafik Planar , tidak ada ujung yang saling bersilangan) shapefile Anda. Ini dapat dilakukan dengan fungsi unary_union dari Shapely

from shapely.ops import unary_union
graph = unary_union(lines)

masukkan deskripsi gambar di sini

gen
sumber
Terima kasih @ gen atas jawaban terperinci. Saya mengedit bagian di mana ia membahas berbagai jenis geometri. Dalam kasus saya, persimpangan juga mengembalikan garis, koleksi geometri, dll. Tapi ini tergantung pada data input. Saya tidak cukup jelas dalam pertanyaan saya.
ustroetz
Jawaban yang bagus Mungkin saya menambahkan bahwa tidak perlu melakukan hal berikut: result = endpts.extend([pt for pt in inters if pt not in endpts])karena tampaknya .extendmetode memodifikasi endpt. Dalam kasus saya result = Nonesetelah operasi itu. Ini endptsyang akhirnya berisi hasil sett
user32882