Meledak yang tumpang tindih dengan poligon baru yang tidak tumpang tindih?

10

Mengingat beberapa poligon yang tumpang tindih dengan berbagai cara, saya ingin mengekspor dari fitur ini semua poligon yang tidak tumpang tindih dengan yang lain, secara iteratif.

Produk akan menjadi sejumlah fitur tanpa tumpang tindih yang, ketika dirangkum bersama, menjadi yang asli.

Produk kemudian dapat digunakan sebagai input ke Statistik Zonal, dan ini akan jauh lebih cepat daripada iterasi Statistik Zonal atas setiap poligon.

Saya telah mencoba kode ini di ArcPy tanpa hasil.

Apakah kode untuk melakukan ini sudah ada?

ndimhypervol
sumber
Maksud Anda, Anda ingin 'meratakan' data ke dalam set topologi yang benar?
nagytech
@Geoist ZonalStats membutuhkan poligon yang tidak tumpang tindih. Ketika Anda memiliki koleksi yang tumpang tindih, solusi yang jelas tetapi tidak efisien adalah untuk mengulangi polys dan menghitung statistik zona satu per satu. Akan lebih efisien untuk memilih subset polys yang tidak tumpang tindih, menerapkan zonalstats padanya, dan beralih. Pertanyaannya adalah bagaimana cara membuat pilihan seperti itu secara efisien.
whuber
whuber - Saya pikir @Geoist menyarankan untuk membuat seperangkat poligon yang tidak tumpang tindih dari persimpangan poligon input. Lihatlah gambar ini - (tidak dapat memposting gambar dalam komentar?). Input ada di sebelah kiri. Seluruh wilayah ditutupi oleh tiga poligon, yang masing-masing berpotongan dengan yang lain. Satu-satunya himpunan bagian yang tidak tumpang tindih adalah lajang dan ini tidak memenuhi persyaratan Gotanuki bahwa serikat pekerja mengisi ruang tersebut. Saya pikir Geoist menyarankan untuk membuat set daerah non-berpotongan di sebelah kanan yang berlaku untuk zonalstats
Llaves
Saya pikir ada beberapa kebingungan mengenai apa produk akhirnya seharusnya. Bisakah Anda memberikan contoh? Interpretasi saya adalah Anda ingin output menjadi pilihan poligon yang tidak tumpang tindih - sambil membuang atau melarutkan poligon yang tersisa. Apakah Anda bekerja dengan satu atau beberapa kelas fitur?
Aaron
1
Bagi saya sepertinya @gotanuki ingin membuat jumlah minimum dari kelas fitur yang hanya berisi poligon yang tidak tumpang tindih dari kelas fitur poligon dengan poligon yang tumpang tindih
PolyGeo

Jawaban:

14

Ini adalah masalah pewarnaan grafik .

Ingatlah bahwa pewarnaan grafik adalah penugasan warna ke simpul grafik sedemikian rupa sehingga tidak ada dua simpul yang berbagi tepi juga akan memiliki warna yang sama. Secara khusus, simpul (abstrak) grafik adalah poligon. Dua simpul terhubung dengan tepi (tidak terarah) setiap kali mereka berpotongan (sebagai poligon). Jika kita mengambil solusi untuk masalah ini - yang merupakan urutan dari (katakanlah k ) koleksi poligon yang terpisah - dan berikan warna yang unik untuk setiap koleksi dalam urutan, maka kita akan mendapatkan k -warna grafik. . Sangat diinginkan untuk menemukan k kecil .

Masalah ini cukup sulit dan tetap tidak terpecahkan untuk grafik yang berubah-ubah. Pertimbangkan solusi perkiraan yang mudah dikodekan. Algoritma berurutan harus dilakukan. Algoritme Welsh-Powell adalah solusi serakah berdasarkan pada urutan menurun dari titik demi titik. Diterjemahkan ke bahasa poligon asli, pertama-tama urutkan poligon dalam urutan jumlah poligon lain yang tumpang tindih. Agar berhasil, beri poligon pertama warna awal. Pada setiap langkah berturut-turut, cobalah untuk mewarnai poligon berikutnya dengan warna yang ada: yaitu, pilih warna yang tidaksudah digunakan oleh semua tetangga poligon itu. (Ada banyak cara untuk memilih di antara warna yang tersedia; coba salah satu yang paling sedikit digunakan atau pilih salah satu secara acak.) Jika poligon berikutnya tidak dapat diwarnai dengan warna yang ada, buat warna baru dan warnai dengan itu.

Setelah Anda mencapai pewarnaan dengan sejumlah kecil warna, lakukan zonalstats warna dengan warna: dengan konstruksi, Anda dijamin tidak akan memiliki dua poligon warna yang tumpang tindih.


Berikut kode contoh di R. (Kode Python tidak akan jauh berbeda.) Pertama, kami menggambarkan tumpang tindih di antara tujuh poligon yang ditampilkan.

Peta tujuh poligon

edges <- matrix(c(1,2, 2,3, 3,4, 4,5, 5,1, 2,6, 4,6, 4,7, 5,7, 1,7), ncol=2, byrow=TRUE)

Yaitu, poligon 1 dan 2 tumpang tindih, dan begitu juga poligon 2 dan 3, 3 dan 4, ..., 1 dan 7.

Urutkan simpul berdasarkan derajat menurun:

vertices <- unique(as.vector(edges))
neighbors <- function(i) union(edges[edges[, 1]==i,2], edges[edges[, 2]==i,1])
nbrhoods <- sapply(vertices, neighbors)
degrees <- sapply(nbrhoods, length)
v <- vertices[rev(order(degrees))]

Algoritma pewarnaan berurutan (kasar) menggunakan warna paling awal yang tersedia yang belum digunakan oleh poligon yang tumpang tindih:

color <- function(i) {
  n <- neighbors(i)
  candidate <- min(setdiff(1:color.next, colors[n]))
  if (candidate==color.next) color.next <<- color.next+1
  colors[i] <<- candidate
}

Inisialisasi struktur data ( colorsdan color.next) dan terapkan algoritma:

colors <- rep(0, length(vertices))
color.next <- 1
temp <- sapply(v, color)

Bagi poligon menjadi kelompok-kelompok sesuai dengan warna:

split(vertices, colors)

Output dalam contoh ini menggunakan empat warna:

$`1`
[1] 2 4

$`2`
[1] 3 6 7

$`3`
[1] 5

$`4`
[1] 1

Empat warna poligon

Ini telah mempartisi poligon menjadi empat kelompok yang tidak tumpang tindih. Dalam hal ini solusinya tidak optimal ({{3,6,5}, {2,4}, {1,7}} adalah tiga warna untuk grafik ini). Secara umum solusi yang didapat seharusnya tidak terlalu buruk.

whuber
sumber
Saya tidak yakin apakah ini menjawab pertanyaan, atau apa pertanyaannya, tetapi itu adalah jawaban yang baik.
nagytech
@Geoist Apakah ada cara agar ilustrasi bisa lebih jelas atau menjelaskan masalahnya dengan lebih baik?
whuber
6

Metodologi yang direkomendasikan oleh #whuber mengilhami saya untuk mengambil arah baru, dan di sini adalah solusi sederhana saya, dalam dua fungsi. Yang pertama, disebut countOverlaps, membuat dua bidang, "tumpang tindih" dan "ovlpCount" untuk merekam setiap polis yang tumpang tindih dengannya, dan berapa banyak tumpang tindih yang terjadi. Fungsi kedua, explodeOverlaps, menciptakan bidang ketiga, "expl", yang memberikan bilangan bulat unik untuk setiap kelompok polis yang tidak tumpang tindih. Pengguna kemudian dapat mengekspor fc baru berdasarkan bidang ini. Proses ini dipecah menjadi dua fungsi karena saya pikir alat countOverlaps dapat terbukti berguna dengan sendirinya. Maafkan kecerobohan kode (dan konvensi penamaan yang ceroboh), karena ini cukup awal, tetapi berhasil. Juga pastikan bahwa "idName" bidang mewakili bidang dengan ID unik (hanya diuji dengan ID integer). Terima kasih Whuber untuk menyediakan saya dengan kerangka kerja yang diperlukan untuk mendekati masalah ini!

def countOverlaps(fc,idName):
    intersect = arcpy.Intersect_analysis(fc,'intersect')
    findID = arcpy.FindIdentical_management(intersect,"explFindID","Shape")
    arcpy.MakeFeatureLayer_management(intersect,"intlyr")
    arcpy.AddJoin_management("intlyr",arcpy.Describe("intlyr").OIDfieldName,findID,"IN_FID","KEEP_ALL")
    segIDs = {}
    featseqName = "explFindID.FEAT_SEQ"
    idNewName = "intersect."+idName

    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        featseqVal = row.getValue(featseqName)
        segIDs[featseqVal] = []
    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        featseqVal = row.getValue(featseqName)
        segIDs[featseqVal].append(idVal)

    segIDs2 = {}
    for row in arcpy.SearchCursor("intlyr"):
        idVal = row.getValue(idNewName)
        segIDs2[idVal] = []

    for x,y in segIDs.iteritems():
        for segID in y:
            segIDs2[segID].extend([k for k in y if k != segID])

    for x,y in segIDs2.iteritems():
        segIDs2[x] = list(set(y))

    arcpy.RemoveJoin_management("intlyr",arcpy.Describe(findID).name)

    if 'overlaps' not in [k.name for k in arcpy.ListFields(fc)]:
        arcpy.AddField_management(fc,'overlaps',"TEXT")
    if 'ovlpCount' not in [k.name for k in arcpy.ListFields(fc)]:
        arcpy.AddField_management(fc,'ovlpCount',"SHORT")

    urows = arcpy.UpdateCursor(fc)
    for urow in urows:
        idVal = urow.getValue(idName)
        if segIDs2.get(idVal):
            urow.overlaps = str(segIDs2[idVal]).strip('[]')
            urow.ovlpCount = len(segIDs2[idVal])
        urows.updateRow(urow)

def explodeOverlaps(fc,idName):

    countOverlaps(fc,idName)

    arcpy.AddField_management(fc,'expl',"SHORT")

    urows = arcpy.UpdateCursor(fc,'"overlaps" IS NULL')
    for urow in urows:
        urow.expl = 1
        urows.updateRow(urow)

    i=1
    lyr = arcpy.MakeFeatureLayer_management(fc)
    while int(arcpy.GetCount_management(arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",'"expl" IS NULL')).getOutput(0)) > 0:
        ovList=[]
        urows = arcpy.UpdateCursor(fc,'"expl" IS NULL','','','ovlpCount D')
        for urow in urows:
            ovVal = urow.overlaps
            idVal = urow.getValue(idName)
            intList = ovVal.replace(' ','').split(',')
            for x in intList:
                intList[intList.index(x)] = int(x)
            if idVal not in ovList:
                urow.expl = i
            urows.updateRow(urow)
            ovList.extend(intList)
        i+=1
ndimhypervol
sumber
2
Untuk menghubungkan ini ke solusi saya: Anda countOverlapsberkorespondensi dengan dua baris nbrhoods <- sapply(vertices, neighbors); degrees <- sapply(nbrhoods, length)dalam kode saya: degreesadalah jumlah tumpang tindih. Tentu saja kode Anda lebih panjang karena mencerminkan sebagian besar analisis GIS yang diterima begitu saja dalam solusi saya: yaitu, bahwa Anda pertama kali mengidentifikasi poligon yang tumpang tindih, dan pada akhirnya Anda menggunakan solusi untuk menghasilkan kumpulan data poligon. Ini akan menjadi ide yang baik untuk merangkum perhitungan grafik-teori, jadi jika Anda pernah menemukan algoritma pewarnaan yang lebih baik, akan lebih mudah untuk dihubungkan.
whuber
1

Sudah lama, tapi saya menggunakan kode ini untuk aplikasi saya sendiri dan sudah bekerja dengan baik - terima kasih. Saya menulis ulang sebagian untuk memperbaruinya, menerapkannya pada baris (dengan toleransi) dan mempercepatnya secara signifikan (di bawah - saya menjalankannya pada 50 juta buffer yang bersinggungan dan hanya membutuhkan beberapa jam).

def ExplodeOverlappingLines(fc, tolerance, keep=True):
        print('Buffering lines...')
        idName = "ORIG_FID"
        fcbuf = arcpy.Buffer_analysis(fc, fc+'buf', tolerance, line_side='FULL', line_end_type='FLAT')
        print('Intersecting buffers...')
        intersect = arcpy.Intersect_analysis(fcbuf,'intersect')

        print('Creating dictionary of overlaps...')
        #Find identical shapes and put them together in a dictionary, unique shapes will only have one value
        segIDs = defaultdict(list)
        with arcpy.da.SearchCursor(intersect, ['Shape@WKT', idName]) as cursor:
            x=0
            for row in cursor:
                if x%100000 == 0:
                    print('Processed {} records for duplicate shapes...'.format(x))
                segIDs[row[0]].append(row[1])
                x+=1

        #Build dictionary of all buffers overlapping each buffer
        segIDs2 = defaultdict(list)
        for v in segIDs.values():
            for segID in v:
                segIDs2[segID].extend([k for k in v if k != segID and k not in segIDs2[segID]])

        print('Assigning lines to non-overlapping sets...')
        grpdict = {}
        # Mark all non-overlapping one to group 1
        for row in arcpy.da.SearchCursor(fcbuf, [idName]):
            if row[0] in segIDs2:
                grpdict[row[0]] = None
            else:
                grpdict[row[0]] = 1

        segIDs2sort = sorted(segIDs2.items(), key=lambda x: (len(x[1]), x[0])) #Sort dictionary by number of overlapping features then by keys
        i = 2
        while None in grpdict.values(): #As long as there remain features not assigned to a group
            print(i)
            ovset = set()  # list of all features overlapping features within current group
            s_update = ovset.update
            for rec in segIDs2sort:
                if grpdict[rec[0]] is None: #If feature has not been assigned a group
                    if rec[0] not in ovset: #If does not overlap with a feature in that group
                        grpdict[rec[0]] = i  # Assign current group to feature
                        s_update(rec[1])  # Add all overlapping feature to ovList
            i += 1 #Iterate to the next group

        print('Writing out results to "expl" field in...'.format(fc))
        arcpy.AddField_management(fc, 'expl', "SHORT")
        with arcpy.da.UpdateCursor(fc,
                                   [arcpy.Describe(fc).OIDfieldName, 'expl']) as cursor:
            for row in cursor:
                if row[0] in grpdict:
                    row[1] = grpdict[row[0]]
                    cursor.updateRow(row)

        if keep == False:
            print('Deleting intermediate outputs...')
            for fc in ['intersect', "explFindID"]:
                arcpy.Delete_management(fc)
messamat
sumber
-3

Dalam hal ini saya biasanya menggunakan metode berikut:

  • Lulus kelas Fitur melalui UNION; (Ini memecah poligon di semua persimpangannya)
  • Tambahkan bidang X, Y dan Area dan hitung;
  • Larutkan hasilnya dengan X, Y, bidang Area.

Saya percaya hasilnya akan menjadi yang Anda inginkan, dan Anda bahkan dapat menghitung jumlah yang tumpang tindih. Tidak tahu apakah dalam hal kinerja akan lebih baik untuk Anda atau tidak.

Alexandre Neto
sumber
2
metode ini tidak mengarahkan Anda ke produk yang diinginkan, yang merupakan serangkaian pilihan minimal atau kelas fitur unik dari sumber asli yang tidak tumpang tindih. produk akan dimasukkan ke dalam statistik zona, dan oleh karena itu menjaga geometri asli setiap fitur sangat penting.
ndimhypervol
Anda benar, maaf. Saya tidak mengerti pertanyaan itu dengan baik. Dalam hal itu, dan tergantung pada ukuran raster, saya biasanya akan mengkonversi raster ke kelas fitur titik sementara (setiap sel titik), dan melakukan sambungan spasial antara itu dan lapisan poligon. Mungkin ini pendekatan yang sangat sederhana, dan kinerja tidak ramah tetapi berfungsi dan poligon yang tumpang tindih tidak akan memberi Anda masalah.
Alexandre Neto
Jika saya mengerti dengan benar apa yang Anda maksud dengan penggabungan spasial ini, solusi kedua Anda masih tidak akan berfungsi, Alexandre, karena ada banyak-ke-banyak hubungan antara titik dan poligon. Bagaimanapun, untuk setiap raster yang cukup besar, pendekatan berbasis vektor ini akan sangat tidak efisien dan untuk raster besar tidak mungkin dilakukan.
whuber
@whuber Anda benar tentang proses yang sangat lambat (Buat saya sekitar setengah jam dengan 4284 x 3009 raster dan 2401 poligon, dalam 2,8GHz dualcore, RAM 3Gb dengan vista). Tapi itu berhasil, karena saya sudah mengujinya. Dalam Gabung Spasial Anda harus menggunakan hubungan satu ke satu, dan mengagregasi nilai raster (seperti rata-rata, Jumlah, dll ...). Hasilnya akan menjadi layer poligon vektor yang mirip dengan aslinya tetapi dengan kolom baru dengan nilai raster agregat yang memotong setiap poligon. Bukan menjadi solusi optimal, ini bisa berguna bagi seseorang dengan keterampilan pemrograman yang lebih sedikit (seperti saya :-)).
Alexandre Neto