Учитывая дыры / ограничения в создании многоугольника Вороного в QGIS?

Я пытаюсь создать в QGIS полигоны вороной, которые бы рассматривали «дыры» в общей области. Примером может быть:

На самом деле я создал Вороного на этом изображении, используя QGIS через команду GRASS, а затем с помощью инструмента «Разница» для создания отверстий. В качестве слоя «Разница» использовался отдельный шейп-файл многоугольника, который содержит размеры отверстий. Примером приложения будет создание полигонов вокруг точек отбора проб, которые были собраны между структурами, которые должны быть исключены из анализа.

Здесь возникают две проблемы:

1. Функция «разности» не работает должным образом на 100%, при этом некоторые границы многоугольника простираются в «дыры». Это можно исправить, найдя строку в таблице атрибутов, у которой нет номера ID полигона (или ID «0»).

2. Этот тип "пробивания отверстий" после факта может привести к прерывистым многоугольникам, как показано красной стрелкой на изображении.

У меня вопрос: есть ли инструмент или плагин Вороного, который может рассматривать наличие «дырок» в центре области как одношаговый процесс, а также устранять образование разрывных многоугольников? Я предполагаю, что такой инструмент будет расширять границу многоугольника до ближайшего пересечения с другой границей, если только эта другая граница не ударяет сначала о границу «дыры».

Это может быть возможно с использованием растров. Сначала преобразуйте точки и полигоны в растр высокого разрешения. Установите маску для ваших границ, используя r.mask. Затем бегите r.grow.distanceв GRASS и используйте Value= output. Это даст вам за каждый пиксель, который является ближайшей точкой. Преобразуйте это обратно в векторные многоугольники. Там могут быть дополнительные шаги, необходимые для избавления от осколков полигонов.

Это, конечно, возможно с растрами.

Надеемся, что этот скриншот показывает проблему более четко. Часть В вороноя «ближе к лету» ближе к первоначальному центру вороноя, но это не принимает во внимание тот факт, что для обхода здания потребуется больше времени. Мое понимание вопроса ОП заключается в том, что вороной должен учитывать это дополнительное расстояние, чтобы пройтись вокруг здания.

Мне нравится предложение от @Guillaume. Однако, когда я попробовал это, у меня были проблемы с получением r.grow.distanceмаски (см. Ниже. Рябь не должна проходить через здания).

Мои знания о траве не так сильны, как могли бы, поэтому, возможно, я делаю что-то глупое. Определенно, сначала проверьте это предложение, так как оно будет намного меньше работы, чем мое ;-)

Шаг 1 - Создайте поверхность затрат

Первый шаг - создать поверхность затрат. Это нужно сделать только один раз.

• создать редактируемый слой, отверстия и все.
• добавьте поле с именем 'unit', установите его в 1.
• используя полигон-растр на «вырубленном» векторном слое (который имеет отверстия), используя поле «unit». Теперь у вас есть слой «маска», где 1 - свободное пространство, а 0 - здание.

• используйте растровый калькулятор, чтобы превратить это в поверхность затрат. Я установлю «на улице» на 1 и «в помещении» на 9999. Это затруднит передвижение по зданиям.

  (( "Маска @ 1" = 1) * 1) + (( "маска @ 1" = 0) * 9999)

Вы можете получить более «естественные» результаты, добавив немного шума на поверхность затрат (например, используйте случайное число от 1 до 3, а не просто 1 для наружных пикселей).

Шаг 2. Создайте растры совокупной стоимости для каждого вороного центра

Теперь мы можем запустить (для одной вороной ячейки за раз) алгоритм GRASS r.cost.coordinatesпротив нашего поверхностного слоя затрат.

Для начальной координаты используйте центр vornoi. В качестве конечной координаты выберите один из углов вашей области. Я предлагаю использовать «Knights Tour», поскольку это дает более плавные результаты.

Результат показывает линии равного времени прохождения от одного вороного центра. Обратите внимание, как полосы обвивают здания.

Не уверен, как лучше всего это автоматизировать. Может быть, обработка в пакетном режиме или в Pyqgis.

Шаг 3. Слить растры

Это, вероятно, потребуется код. Алгоритм будет

create a raster 'A' to match the size of your cumulative cost images
fill raster 'A' with a suitably high number e.g. 9999
create an array of the same size as the raster.
for each cumulative cost raster number 1..N
    for each cell in image
        if cell < value in raster 'A'
            set value in raster 'A' to cell value
            set corresponding cell in array to cum. cost image number
write out array as a raster

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

Затем вы можете использовать растр-полигон. Затем вы можете использовать плагин Generalize для удаления артефактов эффекта «step» из растра.

Приносим извинения за неопределенность на шагах 2 и 3 ... Я надеюсь, что кто-то присоединится к более элегантному решению :)

Примечание № 1 : Я не смог воспроизвести предложенную проблему, потому что инструмент Разница хорошо работал для меня в нескольких тестах, которые я выполнил (возможно, это было связано с простой геометрией проблемы или потому что инструмент был улучшен, так как вопрос был спросил 1 год назад).

Тем не менее, я предлагаю обходной путь в PyQGIS, чтобы избежать использования инструмента Разница . Все основано на локальном пересечении двух входных слоев (см. Рисунок ниже):

1. векторный слой многоугольника, представляющий многоугольники Вороного;
2. векторный слой многоугольника, представляющий дыры / ограничения, которые необходимо исключить из анализа.

Примечание № 2 : Поскольку я не хочу использовать инструмент « Разница» , я не могу избежать создания «осколков» (см. Затем), поэтому мне нужно было запустить v.cleanинструмент для их устранения. Кроме того, как сказал @Chris W,

[...] но второе, приводящее к осколкам, не является совершенно неожиданным - в конце концов, эта область все равно будет выделена той же точке, даже если дыра присутствует. Вы можете использовать этот метод, а затем включить ленты в их соседей с помощью метода очистки .

После этих необходимых условий я выкладываю свой код:

##Voronoi_Polygons=vector polygon
##Constraints=vector polygon
##Voronoi_Cleaned=output vector

from qgis.core import *

voronoi = processing.getObject(Voronoi_Polygons)
crs = voronoi.crs().toWkt()
ex = voronoi.extent()
extent = '%f,%f,%f,%f' % (ex.xMinimum(), ex.xMaximum(), ex.yMinimum(), ex.yMaximum())

constraints = processing.getObject(Constraints)

# Create the output layer
voronoi_mod = QgsVectorLayer('Polygon?crs='+ crs, 'voronoi' , 'memory')
prov = voronoi_mod.dataProvider()
fields = voronoi.pendingFields() # Fields from the input layer
prov.addAttributes(fields) # Add input layer fields to the outLayer
voronoi_mod.updateFields()

# Spatial index containing all the 'constraints'
index_builds = QgsSpatialIndex()
for feat in constraints.getFeatures():
    index_builds.insertFeature(feat)

final_geoms = {}
final_attrs = {}

for feat in voronoi.getFeatures():
    input_geom = feat.geometry()
    input_attrs = feat.attributes()
    final_geom = []
    multi_geom = input_geom.asPolygon()
    input_geoms = [] # edges of the input geometry
    for k in multi_geom:
        input_geoms.extend(k)
    final_geom.append(input_geoms)
    idsList = index_builds.intersects(input_geom.boundingBox())
    mid_geom = [] # edges of the holes/constraints
    if len(idsList) > 0:
        req = QgsFeatureRequest().setFilterFids(idsList)
        for ft in constraints.getFeatures(req):
            geom = ft.geometry()
            hole = []
            res = geom.intersection(input_geom)
            res_geom = res.asPolygon()
            for i in res_geom:
                hole.extend(i)
                mid_geom.append(hole)
        final_geom.extend(mid_geom)
    final_geoms[feat.id()] = final_geom
    final_attrs[feat.id()] = input_attrs

# Add the features to the output layer
outGeom = QgsFeature()
for key, value in final_geoms.iteritems():
    outGeom.setGeometry(QgsGeometry.fromPolygon(value))
    outGeom.setAttributes(final_attrs[key])
    prov.addFeatures([outGeom])

# Add 'voronoi_mod' to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(voronoi_mod)

# Run 'v.clean'
processing.runalg("grass7:v.clean",voronoi_mod, 2, 0.1, extent, -1, 0.0001, Voronoi_Cleaned, None)

# Remove 'voronoi_mod' to the Layers panel
QgsMapLayerRegistry.instance().removeMapLayer(voronoi_mod)

что приводит к такому результату:

Просто для ясности, это будет результат без использования v.cleanинструмента:

Разница с результатом @LeaningCactus заключается в том, что к настоящему времени геометрии не нарушены и их можно «очистить» без ошибок .