Симметрическая разность полигонов

[SMOuk96]

Доброго дня уважаемые форумчане! Столкнулся с задачей по созданию модуля на Python для Qgis. Необходимо обойти все полигональные объекты shape-файла и применить к ним операцию симметрической разности. Пример такого shape-файла приложил

WV2_B45_20160625_C_NDVI_F165_S_Sieve_Convexhull-9.shp

(12.58 КБ) 926 скачиваний

.
Нашел в qgis в модуле ftools функцию Симметрической разность, она работает с исходным слоем и слоем разности. Можно ли как-то посмотреть код данной функции, чтобы переделать ее под работу с одним слоем И двигаюсь ли я в нужно направлении. Прошу вашей помощи потому что новичок в этой теме и зашел в тупик. Help me, plz
На данный момент сам сделал (по примерам из интернета) загрузку слоев в comboBox.

[Александр Мурый]

Вы выложили только часть шейп-файла (сам .shp). Ещё надо как минимум .shx и .dbf.

Необходимо обойти все полигональные объекты shape-файла и применить к ним операцию симметрической разности

Вам нужно получить пересечения всех полигонов в одном слое? Если так, используйте модуль GRASS <v.clean> из панели инструментов ("Анализ данных"/Processing). Опция "bpol".

[SMOuk96]

Результат работы скрипта должен быть примерно таким (сделал вручную)

Screenshot_2.png (109.71 КБ) 9766 просмотров

Если так, используйте модуль GRASS <v.clean> из панели инструментов ("Анализ данных"/Processing). Опция "bpol".

Результат работы bpol:

Screenshot_3.png (50.42 КБ) 9766 просмотров

Думал он оставит только пересечения, но нет.
Александр Мурый, возможно ли получить пересечения полигонов и потом от исходного слоя отнять данное пересечение, чтобы и получить симметрическую разность?
Т. к. делаю скрипт, то вызывать команды надо оттуда, а не из интерфейса, а как вызвать Опцию "blop" модуля GRASS <v.clean>?
Возможно что-то такое? processing.runalg('qgis:blop',исходный файл,True,id,'файл для записи') ?

[Александр Мурый]

Сделал простую модель для Processing - "Polygon self symmetrical difference". Используется алгоритм SAGA "Polygon self intersection" и извлечение по атрибутам.
Экспорт модели в Python-скрипт:

Код: Выделить всё

##Polygon self symmetrical difference=name
##poly=vector
##polygonselfsymmetricaldifference=output vector
outputs_QGISADDUNIQUEVALUEINDEXFIELD_1=processing.runalg('qgis:adduniquevalueindexfield', poly,'diff_id',None)
outputs_SAGAPOLYGONSELFINTERSECTION_1=processing.runalg('saga:polygonselfintersection', outputs_QGISADDUNIQUEVALUEINDEXFIELD_1['OUTPUT'],'diff_id',None)
outputs_QGISEXTRACTBYEXPRESSION_1=processing.runalg('qgis:extractbyexpression', outputs_SAGAPOLYGONSELFINTERSECTION_1['INTERSECT'],'"ID" not like '%|%'',polygonselfsymmetricaldifference)

Вот файл модели:

[rhot]

Я вообще начал замечать, что SAGA работает в разы быстрее QGIS, потому что последний обрастает ненужным функционалом и плюшками.

[Александр Мурый]

Я вообще начал замечать, что SAGA работает в разы быстрее QGIS, потому что последний обрастает ненужным функционалом и плюшками.

Ну не знаю, алгоритм SAGA "Polygon self intersection" тупит нещадно даже на нескольких полигонах.

[rhot]

Ну не знаю, алгоритм SAGA "Polygon self intersection" тупит нещадно даже на нескольких полигонах.

Линукс или винда? Какая версия саги?

[Александр Мурый]

Линукс или винда? Какая версия саги?

В линуксах 2.3.1, в виндах 2.3.2, и везде проблемы с Сагой. Модель делалась на винде, работает нормально только там пока.

[rhot]

везде проблемы с Сагой. Модель делалась на винде, работает нормально только там пока.

Хмм,.. выложите, пожалуйста, вектор потестить.

[Dmitry Stasev]

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

Можно сиквелами, только не уверен в стабильности на больших объемах

Вот попробовал:

Код: Выделить всё

With CTE AS(
	SELECT '1' AS grp, st_union(st_intersection(a.geometry, b.geometry)) as geometry 
	FROM test a JOIN test b ON (
		(a.id < b.id)and
		(a.geometry is NOT NULL)and
		(b.geometry is NOT NULL)and
		(st_intersects(a.geometry, b.geometry))
	) 
	GROUP BY grp
) 
SELECT st_difference(test.geometry, CTE.geometry) FROM test
JOIN CTE 

1. Тут сам слой - "test"
2. Для корректного поиска пересечений добавил поле "id" - int
3. В обобщенном выражении считаю разницу и потом вычитаю ее из всех объектов
Для экспериментов достаточно вставить SQL в виртуальный слой и поправить под себя "test" и "id"

PS
Сделал для эксперимента (захотелось попробовать одним запросом), поэтому что там с производительностью ...

[Александр Мурый]

Хмм,.. выложите, пожалуйста, вектор потестить.

Тестовый слой рисуется за несколько секунд в кугисе. Дело не в векторе, дело, видимо, в самом модуле Processing. Например, в линуксовой версии QGIS 2.18.* алгоритм SAGA "Polygon self intersection" не работает в принципе почему-то.
А вообще, мы отклонились от темы.

[Dmitry Stasev]

Тестовый слой рисуется за несколько секунд в кугисе

Интересно как мой скриптик отработает, не посмотрите?

[Александр Мурый]

Интересно как мой скриптик отработает, не посмотрите?

На виртуальном слое, сделанном на базе memory-слоя, не взлетело. Не понимает поле "id".

[Dmitry Stasev]

Не понимает поле "id"

конечно я же выше писал, что надо добавить или поменять на свое.
Вот этот python - скрипт должен сам все сделать )
Сделает "SelfSymDifference" слой из активного

Код: Выделить всё

from qgis.PyQt.QtCore import *

srcLayer = iface.activeLayer()
dstLayerName = 'SelfSymDifference'
idName = 'idSelfSymDifference'
srcLayer.startEditing()
myField = QgsField(idName, QVariant.Int )
srcLayer.addAttribute( myField )
srcLayer.updateFields()

idx = srcLayer.fieldNameIndex(idName)
i = 0
for f in srcLayer.getFeatures():
    f[idx] = i
    i += 1
    srcLayer.updateFeature( f )
srcLayer.commitChanges()

tpl = '''
With CTE AS(
    SELECT '1' AS grp, st_union(st_intersection(a.geometry, b.geometry)) as geometry 
    FROM %s a JOIN test b ON (
        (a.%s < b.%s)and
        (a.geometry is NOT NULL)and
        (b.geometry is NOT NULL)and
        (st_intersects(a.geometry, b.geometry))
    ) 
    GROUP BY grp
) 
SELECT st_difference(test.geometry, CTE.geometry) FROM %s
JOIN CTE'''

sql = tpl % (srcLayer.name(), idName, idName, srcLayer.name())

layer = QgsVectorLayer("?query="+sql, dstLayerName, "virtual")
QgsMapLayerRegistry.instance().addMapLayer(layer, True)

ЗЫ
без проверок на валидность и типы слоев.

[SMOuk96]

Сделал простую модель для Processing - "Polygon self symmetrical difference". Используется алгоритм SAGA "Polygon self intersection" и извлечение по атрибутам.

Запускаю ваш скрипт, но к сожалению у меня почему-то нет второго алгоритма SAGAовского, версия QGIS 2.18.9

1.png (13.11 КБ) 9614 просмотров

Вот файл модели:

Впервые сталкиваюсь с таким, попробовал открыть в QT Creator, но особо не понял что куда. Не посоветуете что-то почитать чтобы понимать хотя бы что за файл с расширение .model или как-то "на пальцах"..