Взрывается наложение на новые неперекрывающиеся полигоны?

10

Учитывая несколько полигонов, которые перекрываются несколькими способами, я хотел бы экспортировать из этих объектов все полигоны, которые не перекрываются с другими, итеративно.

Продукт будет иметь ряд функций без каких-либо совпадений, которые при суммировании образуют оригинал.

Затем эти продукты можно использовать в качестве входных данных для зональной статистики, и это будет намного быстрее, чем итерация зональной статистики по каждому полигону.

Я пытался закодировать это в ArcPy без успеха.

Код для этого уже существует?

ndimhypervol
источник
Вы хотите сказать, что хотите «сгладить» данные в топологически правильный набор?
nagytech
@ Geoist ZonalStats требует многоугольников, которые не перекрываются. Когда у вас есть перекрывающаяся коллекция, очевидное, но неэффективное решение состоит в том, чтобы циклически проходить полисы и вычислять зональную статистику одну за другой. Было бы более эффективно выбрать подмножество непересекающихся полисов, применить к ним зональные статистические данные и выполнить итерацию. Вопрос состоит в том, как сделать такой выбор эффективно.
whuber
whuber - я думаю, что @Geoist предлагает создать набор непересекающихся многоугольников из пересечений входных многоугольников. Посмотрите на это изображение - (не можете публиковать изображения в комментариях?). Вход слева. Весь регион покрыт тремя полигонами, каждый из которых пересекает оба других. Единственными неперекрывающимися подмножествами являются синглтоны, и они не удовлетворяют требованию Готануки о том, что объединение заполняет пространство. Я думаю, что Geoist предлагает создать набор непересекающихся областей справа, который действителен для zonalstats
Llaves
Я думаю, что есть некоторая путаница относительно того, каким должен быть конечный продукт. Не могли бы вы привести пример? Моя интерпретация заключается в том, что вы хотели бы, чтобы на выходе был выбор полигонов, которые не перекрываются - при отбрасывании или растворении оставшихся полигонов. Вы работаете с одним или несколькими классами объектов?
Аарон
1
Мне кажется, что @gotanuki хочет создать минимальное количество классов объектов, которые содержат только неперекрывающиеся полигоны из класса полигонов с перекрывающимися полигонами
PolyGeo

Ответы:

14

Это проблема раскраски графа .

Напомним, что раскраска графа - это присвоение цвета вершинам графа таким образом, что никакие две вершины, которые имеют общее ребро, также не будут иметь одинаковый цвет. В частности, (абстрактными) вершинами графа являются многоугольники. Две вершины связаны с (ненаправленным) ребром, когда они пересекаются (как многоугольники). Если мы возьмем какое-либо решение проблемы - которое представляет собой последовательность (скажем, k ) непересекающихся наборов многоугольников - и назначим уникальный цвет каждому набору в последовательности, то мы получим k- раскраску графа. , Желательно найти маленький к .

Эта проблема довольно сложна и остается нерешенной для произвольных графов. Рассмотрим примерное решение, которое просто закодировать. Последовательный алгоритм должен сделать. Алгоритм Уэлша-Пауэлла является жадным решением, основанным на нисходящем упорядочении вершин по степени. В переводе на язык исходных полигонов, сначала отсортируйте полигоны в порядке убывания числа других полигонов, которые они перекрывают. Работая по порядку, дайте первому многоугольнику начальный цвет. На каждом последующем шаге попробуйте закрасить следующий многоугольник существующим цветом: выберите цвет, который не соответствуетуже используется любым из соседей этого многоугольника. (Есть много способов выбрать один из доступных цветов; попробуйте тот, который был наименее использован, или выберите один случайным образом.) Если следующий полигон не может быть окрашен существующим цветом, создайте новый цвет и закрасьте его этим.

Как только вы добились раскраски с небольшим количеством цветов, выполните зональную статистику цвет за цветом: по построению вы гарантируете, что никакие два многоугольника данного цвета не перекрываются.


Вот пример кода в R. (Код Python не будет сильно отличаться.) Во-первых, мы описываем перекрытия между семью показанными полигонами.

Карта семи полигонов

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)

То есть полигоны 1 и 2 перекрываются, как и полигоны 2 и 3, 3 и 4, ..., 1 и 7.

Сортировка вершин по убыванию:

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))]

Алгоритм последовательной раскраски (грубый) использует самый ранний доступный цвет, еще не использованный ни одним перекрывающимся многоугольником:

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
}

Инициализируйте структуры данных ( colorsи color.next) и примените алгоритм:

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

Разделите полигоны на группы по цвету:

split(vertices, colors)

Вывод в этом примере использует четыре цвета:

$`1`
[1] 2 4

$`2`
[1] 3 6 7

$`3`
[1] 5

$`4`
[1] 1

Четырехцветные многоугольники

Он разделил полигоны на четыре непересекающиеся группы. В этом случае решение не является оптимальным ({{3,6,5}, {2,4}, {1,7}} является трехцветной для этого графа). В целом, решение, которое оно получает, не должно быть слишком плохим.

Whuber
источник
Я не уверен, отвечает ли это на вопрос или каков вопрос, но тем не менее это хороший ответ.
nagytech
@ Geoist Есть ли способ, как я мог бы проиллюстрировать иллюстрацию или лучше объяснить проблему?
whuber
6

Методология, рекомендованная #whuber, вдохновила меня на новое направление, и вот мое загадочное решение, состоящее из двух функций. Первый, называемый countOverlaps, создает два поля, «overlaps» и «ovlpCount», чтобы записывать для каждого poly, который перекрывается с polys, и сколько произошло перекрытий. Вторая функция, explodeOverlaps, создает третье поле «expl», которое дает уникальное целое число каждой группе непересекающихся полисов. Пользователь может затем экспортировать новые фк на основе этого поля. Процесс разбит на две функции, потому что я думаю, что инструмент countOverlaps может оказаться полезным сам по себе. Пожалуйста, извините за неряшливость кода (и небрежное соглашение об именах), так как он довольно предварительный, но он работает. Также убедитесь, что «idName» field представляет поле с уникальными идентификаторами (проверяется только с целочисленными идентификаторами). Спасибо, что предоставили мне основу, необходимую для решения этой проблемы!

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
источник
2
Чтобы связать это с моим решением: вы countOverlapsсоответствует двум строкам nbrhoods <- sapply(vertices, neighbors); degrees <- sapply(nbrhoods, length)в моем коде: degreesэто количество совпадений. Конечно, ваш код длиннее, потому что он отражает большую часть анализа ГИС, который считается само собой разумеющимся в моем решении: а именно, что вы сначала идентифицируете, какие полигоны перекрываются, и что в конце вы используете решение для вывода наборов данных полигонов. Было бы неплохо инкапсулировать теоретико-графовые вычисления, поэтому, если вы когда-нибудь найдете лучший алгоритм раскраски, его будет легко подключить.
whuber
1

Это было некоторое время, но я использовал этот код для своего собственного приложения, и он работал отлично - спасибо. Я переписал его часть, чтобы обновить, применить к строкам (с допуском) и значительно ускорить его (ниже - я запускаю его на 50 миллионах пересекающихся буферов, и это занимает всего пару часов).

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
источник
-3

В этом случае я обычно использую следующий метод:

  • Передайте класс объектов через UNION; (Это ломает многоугольники во всех его пересечениях)
  • Добавьте поля X, Y и Area и рассчитайте их;
  • Растворить результат по полям X, Y, Area.

Я считаю, что результатом будет именно тот результат, который вы хотели получить, и вы даже можете подсчитать количество совпадений. Не знаю, будет ли это с точки зрения производительности лучше для вас или нет.

Александр Нето
источник
2
этот метод не позволяет получить нужный продукт, представляющий собой минимальную серию выборок или уникальных классов объектов оригинала, которые не перекрываются. продукты будут введены в зональную статистику, и поэтому важно сохранить исходную геометрию каждого элемента.
ndimhypervol
Вы правы, извините. Я не очень хорошо понял вопрос. В этом случае, в зависимости от размера растра, я обычно преобразовываю растр во временный класс точечных объектов (каждая ячейка - точка) и выполняю пространственное соединение между ним и слоем многоугольника. Может быть, это очень упрощенный и недружественный к производительности подход, но он работает, и перекрывающиеся полигоны не доставят вам никаких проблем.
Александр Нето
Если я правильно понимаю, что вы подразумеваете под этим пространственным объединением, ваше второе решение все равно не сработает, Александр, потому что между точками и полигонами существует отношение многих ко многим. Независимо от этого, для любого значительного растра этот векторный подход будет крайне неэффективным, а для больших растров его будет невозможно реализовать.
whuber
@whuber Вы правы в том, что очень медленный процесс (мне понадобилось около получаса с 4284 x 3009 растром и 2401 полигоном в двухъядерном 2,8 ГГц, 3 Гб ОЗУ с Vista). Но это работает, как я уже проверял. В Пространственном соединении вы должны использовать отношение один к одному и агрегировать растровые значения (как среднее, сумма и т. Д.). Результатом будет слой векторного многоугольника, аналогичный исходному, но с новым столбцом с агрегированными растровыми значениями, которые пересекают каждый многоугольник. Не являясь оптимальным решением, это может быть полезно для тех, у кого меньше навыков программирования (таких как я :-)).
Александр Нето