Secara terprogram menemukan poligon yang> 90% tumpang tindih dengan layer vektor poligon lain menggunakan QGIS?

9

Contoh lapisan

Saya mencoba mencari cara bagaimana menggunakan python untuk mengekstraksi poligon dalam satu vektor yang tumpang tindih dengan> 90% oleh vektor lain. Saya kemudian ingin memiliki vektor / peta yang hanya akan menunjukkan poligon itu. Contoh gambar menunjukkan lapisan saya. Saya ingin semua poligon abu-abu yang> 90% merah.

Saya perlu melakukan ini semua melalui python (atau metode otomatis serupa). Saya punya ~ 1000 peta untuk diproses dengan cara yang sama.

Luar biasa
sumber
Anda ingin melakukan overlay 'union' (lihat infogeoblog.wordpress.com/2013/01/08/geo-processing-in-qgis untuk beberapa dasar) lalu untuk setiap poligon asli hitung statistik 'in' dan statistik 'out' gis.stackexchange.com/questions/43037/… untuk menentukan persen overlay ... petunjuk: Anda harus memiliki pengukuran area gis.stackexchange.com/questions/23355/…
Michael Stimson
Terima kasih atas tipsnya. Itu pendekatan yang sama yang saya coba. Saya dapat melakukan penyatuan melalui konsol python dengan cukup mudah. Sudah ditambahkan di nilai atribut area. Ini langkah selanjutnya yang saya tidak yakin. Bagaimana cara menggunakan python untuk menghitung statistik 'masuk' dan 'keluar' sehingga saya dapat mengidentifikasi / memilih / memotong dll. Poligon> 90%?
Luar biasa
Saya pikir itu mungkin tanpa python. Apakah Anda memerlukan python absolutly atau solusi dengan lapisan virtual yang baik untuk Anda?
Pierma
Area 'in' akan memiliki atribut dari kedua poligon, area 'out' hanya memiliki atribut dari satu set poligon. Dapatkan kedua set statistik area dan bergabung kembali ke poligon asli, tambahkan bidang untuk 'masuk', 'keluar' dan cakupan, menghitung nilai untuk 'masuk' dan 'keluar' dari jumlah area kemudian bagikan 'masuk' dengan area asli (atau 'dalam' + 'keluar') untuk menghitung persen.
Michael Stimson
1
Pierma - Saya hanya perlu metode otomatis untuk menemukan poligon.
Luar biasa

Jawaban:

3

Kode selanjutnya berfungsi di Konsol Python saya di QGIS. Ini menghasilkan lapisan memori dengan poligon yang> 90% tumpang tindih dengan area merah.

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for polygon_intersects
feats_lyr1 = [ feat for feat in layers[0].getFeatures() ]

#for xwRcl
feats_lyr2 = [ feat for feat in layers[1].getFeatures() ]

selected_feats = []

for i, feat1 in enumerate(feats_lyr1):
    area1 = 0
    area2 = 0
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area = feat1.geometry().intersection(feat2.geometry()).area()
            print i, j, area, feat2.attribute('class')
            if feat2.attribute('class') == 1:
                area1 += area
            else:
                area2 += area
    crit = area1/(area1 + area2)
    print crit
    if crit > 0.9:
        selected_feats.append(feat1)

epsg = layers[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "mem_layer",
                           "memory")

prov = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Saya mencoba kode dengan dua lapisan vektor ini:

masukkan deskripsi gambar di sini

Setelah menjalankan kode di Konsol Python QGIS, untuk hasil yang menguatkan, ada dicetak indeks i, j dari fitur yang terlibat, area persimpangan, atribut bidang dalam polygons_intersects (1 untuk area merah dan 2 untuk area abu-abu) dan kriteria yang tumpang tindih .

0 0 9454207.56892 1
0 1 17429206.7906 2
0 2 10326705.2376 2
0 4 40775341.6814 1
0 5 26342803.0964 2
0 7 11875753.3216 2
0.432253120382
1 6 1198411.02558 2
1 7 1545489.96614 2
1 10 27511427.9909 1
0.90930850584
2 7 750262.940888 2
2 8 12012343.5859 1
0.941213972294
3 6 23321277.5158 2
0.0

Lapisan memori yang dibuat (fitur hijau) dapat diamati pada gambar berikutnya. Seperti yang diharapkan.

masukkan deskripsi gambar di sini

xunilk
sumber
6

Berikut solusi yang tidak memerlukan python.

Tambahkan layer virtual baru dengan permintaan seperti:

WITH r AS (
SELECT 
    Basins800.rowid AS idGray, 
    area(Basins800.geometry) AS areaGray, 
    area(Intersection(Basins800.geometry, Severity.geometry)) AS aeraInter, 
    Basins800.geometry AS geomGray 
  FROM Basins800, Severity
)

SELECT *, areaInterSum/areaGray  AS overlap , geomGray 
    FROM (
        SELECT 
           idGray, 
           areaGray, 
           sum(areaInter) AS areaInterSum, 
           geomGray 
        FROM r 
        GROUP BY idGray) 
     WHERE areaInterSum/areaGray > 0.9

Dengan:

  • Basins800 sebagai layer yang ingin Anda filter dengan poligon abu-abu

  • Tingkat keparahan: lapisan merah Anda tumpang tindih.

Hasilnya akan menjadi layer baru dengan hanya semua poligon abu-abu> 90% tumpang tindih oleh poligon merah, dengan bidang baru yang berisi persen tumpang tindih.

masukkan deskripsi gambar di sini

Semoga ini berhasil. Saya dapat menambahkan rincian lebih lanjut tentang permintaan jika diperlukan.

Catatan: Data Anda mengandung poligon yang sangat kecil (berasal dari pemrosesan raster Anda dan sesuai dengan piksel raster (pada gambar, kita dapat melihat 4 poligon tetapi ada 25 poligon kecil lainnya). Ini membuat kueri sangat lambat dijalankan (fungsi titik-temu) menghasilkan satu fitur untuk setiap pasangan fitur dari dua lapisan).

Pierma
sumber
Saya mendapatkan kesalahan saat menjalankan kueri melalui tombol 'buat lapisan virtual'. "Kesalahan eksekusi permintaan pada CREATE TEMP VIEW _tview AS WITH r AS (" .... sisa kode di sini ... diikuti oleh: "1 - near" WITH ": syntax error" Saya cukup baru untuk QGIS. Dapatkah saya buat layer virtual ini secara programatis? Terima kasih atas bantuan Anda!
dnormous
Berikut tautan untuk mengunduh shapefile: tautan
dnormous
Maaf, salinan buruk antara abu-abu dan abu-abu (maaf untuk bahasa Inggris aproximative saya). Saya mengedit kueri. Seharusnya berfungsi sekarang. Mengapa Anda ingin membuat layer secara progamatic? Kelebihan dari lapisan virtual adalah tidak merusak, dan jika Anda mengedit data Anda (poligon abu-abu atau merah), lapisan virtual akan diperbarui secara otomatis.
Pierma
Ini hanya satu bagian kecil dari proses. Saya memiliki ~ 1000 peta ini untuk dikerjakan, jadi mengotomatiskan prosesnya akan sangat membantu.
Luar biasa
Saya masih mendapatkan kesalahan yang sama -> "1 - dekat" WITH ": kesalahan sintaks". Saya menghubungkan nama-nama lokal untuk setiap lapisan di untuk grayLayer dan redLayer. Apakah nama lokal yang harus saya gunakan? yaitu: lapisan abu-abu dilabeli sebagai "Basins_800", jadi saya punya kode seperti "Basins_800.geometry"
dnormous
2

Setelah melihat tautan ke shapefile Severity dan Basins800 , saya bisa memahami geoproses yang diperlukan. Saya mengubah kode di:

Secara terprogram menemukan poligon yang> 90% tumpang tindih dengan layer vektor poligon lain menggunakan QGIS?

untuk mendapatkan yang ini:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for Severity
feats_lyr1 = [ feat for feat in layers[0].getFeatures() ]

#for Basins800
feats_lyr2 = [ feat for feat in layers[1].getFeatures() ]

selected_feats = []

print "processing..."

for i, feat1 in enumerate(feats_lyr1):
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area1 = feat1.geometry().intersection(feat2.geometry()).area()
            area2 = feat1.geometry().area()
            print i, j, area1, area2
    crit = area1/area2
    print crit
    if crit > 0.9:
        selected_feats.append(feat1)

epsg = layers[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "mem_layer",
                           "memory")

prov = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Setelah menjalankan kode, dengan shapefile ini di Python Console dari QGIS, dalam beberapa menit saya mendapat hasil yang mirip dengan Pierma ; di mana lapisan memori memiliki 31 fitur (berbeda dari 29 poligon yang didapatnya).

masukkan deskripsi gambar di sini

Saya tidak akan men-debug hasil karena ada 1901 * 3528 = 6706728 interaksi untuk fitur. Namun, kodenya terlihat menjanjikan.

xunilk
sumber