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

[SMOuk96]

Вот этот python - скрипт должен сам все сделать )

Дмитрий, запустил ваш скрипт после пары import'ов, но к сожалению не увидел визуального изменения слоя, может быть что-то делаю не так? Поле id - есть, в скрипте изменил слой test на свой

Спойлер

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

from qgis.PyQt.QtCore import *
from qgis.utils import iface
from qgis.core 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 WV2_B45_20160625_C_NDVI_F165_S_Sieve_Convexhull-9 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(WV2_B45_20160625_C_NDVI_F165_S_Sieve_Convexhull-9.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)

1.png (12.36 КБ) 8027 просмотров

- таблица атрибутов после работы скрипта

2.png (22.99 КБ) 8027 просмотров

- лог скрипта, как видно отработал он успешно
Что могло пойти не так?

[SMOuk96]

полный shape-файл

WV2_B45_20160625_C_NDVI_F165_S_Sieve_Convexhull-9.shx

(420 байт) 731 скачивание

WV2_B45_20160625_C_NDVI_F165_S_Sieve_Convexhull-9.shp

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

WV2_B45_20160625_C_NDVI_F165_S_Sieve_Convexhull-9.dbf

(5.34 КБ) 553 скачивания

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

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

Традиционно: http://www.lmgtfy.com/?q=qgis%20processing%20modeler

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

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

Конечно же, я умею читать, но скрипт просто не работает. И ошибок не выдаёт.

[Dmitry Stasev]

Похоже одна из причин - длинное имя слоя и sqlite его не может обработать
тоже и с моим добавляемым полем (уменьшил его до "idSSD")
проверил на данных ТС, только имя слоя поменял на короткое, скрипт отработал

Спойлер

from qgis.PyQt.QtCore import *

srcLayer = iface.activeLayer()
dstLayerName = 'SelfSymDifference'
idName = 'idSSD'
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 %s 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(%s.geometry, CTE.geometry) FROM %s
JOIN CTE'''

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

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

не увидел визуального изменения слоя

Там создается новый виртуальный слой "dstLayerName = 'SelfSymDifference' "

[SMOuk96]

Dmitry Stasev, все замечательно получилось, добавил алгоритм в модуль к алгоритму для приведения полигонов к выпуклому виду

Спойлер

До обработки.jpg (488 КБ) 7973 просмотра

После обработки симметричной разностью.jpg (192.19 КБ) 7973 просмотра

Ваш алгоритм хорошо отработал

Симметричная разность + приведение к выпуклой форме.jpg (484.59 КБ) 7973 просмотра

А вот приведение к выпуклому виду слегка подкачало, хотя может быть потому что от "SelfSymDifference" остались остаточные линии..

[SMOuk96]

Вот такие вот линии остались почему-то

Спойлер

3.jpg (77.47 КБ) 7972 просмотра

Если интересно, вот код скрипта приведения к выпуклой форме:

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

#!/usr/bin/python
# -*- coding: UTF-8 -*-

from osgeo import ogr
import sys
import os
import math

currentPath = os.path.dirname(sys.argv[0])
# =====================================================

# *****************************************************
def convexHull(inputShapefile, outShapefile = currentPath + "/temp/convexhull.shp", inThickness='0', inArea='0'):
    # Open base SHP File
    inputData = ogr.Open(inputShapefile, False)
    if inputData is None:
        print("ERROR: open failed")
        sys.exit(1)
    # Get once layer from SHP File
    inputLayer = inputData.GetLayer(0)
    if inputLayer is None:
        print("ERROR: reading file")
        sys.exit(1)

    filter = ""
    if inThickness != '0':
        filter += "thickness <= " + inThickness
    if inArea != '0':
        filter += " and area <= " + inArea
    # "thickness <= 0.2 and area <= 50"
    if filter != "":
        print(filter)
        inputLayer.SetAttributeFilter(filter)

    # Dostup k atributam
    inFeatDef = inputLayer.GetLayerDefn()

    # Create new SHP file
    drv = ogr.GetDriverByName("ESRI Shapefile")

    # Create exit dataset
    # outShapefile = currentPath + "/temp/states_convexhull.shp"

    # Remove output shapefile if it already exists
    if os.path.exists(outShapefile):
        drv.DeleteDataSource(outShapefile)

    outputData = drv.CreateDataSource(outShapefile)
    # Create exit layer
    outputLayer = outputData.CreateLayer(outShapefile, None, ogr.wkbPolygon)

    # Add new attributes to output layer
    idField = ogr.FieldDefn("id", ogr.OFTInteger)

    areaField = ogr.FieldDefn("area", ogr.OFTReal)
    areaField.SetWidth(10)
    areaField.SetPrecision(6)

    perimeterField = ogr.FieldDefn("perimeter", ogr.OFTReal)
    perimeterField.SetWidth(10)
    perimeterField.SetPrecision(4)

    thicknessField = ogr.FieldDefn("thickness", ogr.OFTReal)
    thicknessField.SetWidth(10)
    thicknessField.SetPrecision(4)

    outputLayer.CreateField(idField)
    outputLayer.CreateField(areaField)
    outputLayer.CreateField(perimeterField)
    outputLayer.CreateField(thicknessField)

    # Create the feature and set values
    featureDefn = outputLayer.GetLayerDefn()
    feature = ogr.Feature(featureDefn)

    # ID value
    k = 0
    # chitaem postrochno file shp
    for inFeature in inputLayer:
        k = k + 1
        geom = inFeature.GetGeometryRef()
        # print(geom)
        # Out geometry type
        # print geom.GetGeometryName()
        # Calculate Convex Hull
        convexhull = geom.ConvexHull()
        # print(convexhull)
        # Area in meters kvadr
        areaImKvadMeters = convexhull.GetArea()
        # Area in kilometers kvadr
        # area = area / 1000000
        # Area in Gektar
        area = areaImKvadMeters / 10000
        # Perimeter in meters
        perimeterInMeters = convexhull.Boundary().Length()
        # Perimeter in kilometers
        perimeter = perimeterInMeters / 1000
        # Thickness in range 0..1
        thickness = (4 * 3.14) * (areaImKvadMeters / (perimeterInMeters**2))

        feature.SetField("id", k)
        # feature.SetField("CLASS", cl)
        feature.SetField("area", area)
        feature.SetField("perimeter", perimeter)
        feature.SetField("thickness", thickness)
        feature.SetGeometry(convexhull)
        outputLayer.CreateFeature(feature)

    feature = None

    # Save and close DataSource
    inputData.Destroy()
    outputData.Destroy()
# =====================================================

if __name__ == "__main__":
    if len(sys.argv) < 5:
        print ("\n\nERROR PARAMETERS. EXAMPLE: \ngeom_convexhull.py C:/NDVI/ndvi/input.shp C:/NDVI/ndvi/output.shp [thikness_value] [area_value]\n")
        sys.exit(1)
    else:
        inputSHPPath = sys.argv[1]
        outputSHPPath = sys.argv[2]
        if sys.argv[3] != None or sys.argv[4] != None:
            inThickness = sys.argv[3]
            inArea = sys.argv[4]
            convexHull(inputSHPPath, outputSHPPath, inThickness, inArea)
        else:
            convexHull(inputSHPPath, outputSHPPath)
        # print("P1: {}, P2: {}, P3: {}, P4: {}".format(sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4]))

Запускается через OSGeo4W Shell, принимает на вход (исходный слой, выходной слой, thickness, area)
thickness - толщина
area - площадь

[SMOuk96]

Там создается новый виртуальный слой "dstLayerName = 'SelfSymDifference' "

А не могли бы немного подробнее описать алгоритм? Начиная с цикла не очень понимаю что происходит?

[Dmitry Stasev]

Вот такие вот линии остались почему-то

Операции с плавающей точкой видно...
Можно добавить -буфер +буфер (при маленьком значении на геометрии почти не повлияет), время добавит, но выбросы съест:

Спойлер

# -*- coding: utf-8 -*-

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

from qgis.PyQt.QtCore import *

#работаем с активным слоем
srcLayer = iface.activeLayer()
#имя выходного слоя
dstLayerName = 'SelfSymDifference'
#имя добавленного поля id (для SQL, чтобы выбрать "неповторяющиеся" пересечения)
idName = 'idSSD'
#включаем редактирование слоя
srcLayer.startEditing()
#создаем поле (атрибут)
myField = QgsField(idName, QVariant.Int )
#добавляем атрибут к таблице
srcLayer.addAttribute( myField )
#обновляем
srcLayer.updateFields()
#смотрим индекс нашего поля
idx = srcLayer.fieldNameIndex(idName)
#создаем поле (атрибут)
#в цикле просто в записи(фичи/объекты) добавляем в наше поле порядковый № (эмулирую id)
i = 0
for f in srcLayer.getFeatures():
    f[idx] = i
    i += 1
    srcLayer.updateFeature( f )
#запоминаем изменения
srcLayer.commitChanges()

#это шаблон запроса
#если с SQL совсем никак, то тут и не обясню
#в первых 10 строках считаю объект-пересечение
#в 2 последних из каждого в исх таблице вычитаю объект-пересечение
tpl = '''
With CTE AS(
    SELECT '1' AS grp, st_union(st_intersection(a.geometry, b.geometry)) as geometry 
    FROM %s a JOIN %s 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  id, CLASS, class_mean, st_buffer(st_buffer(st_difference(%s.geometry, CTE.geometry), -0.000001), 0.000001) as geometry FROM %s
JOIN CTE'''
#подставляю данные (имена таблиц/полей) в шаблон
sql = tpl % (srcLayer.name(), srcLayer.name(), idName, idName, srcLayer.name(), srcLayer.name())
#создаю, на основе SQL, виртуальный слой
layer = QgsVectorLayer("?query="+sql, dstLayerName, "virtual")
#добавляю виртуальный слой в проект
QgsMapLayerRegistry.instance().addMapLayer(layer, True)

001.png (24.19 КБ) 7951 просмотр

Добавил комменты и вывод полей (трех как пример)
там в скрипте еще и ошибочка вкралась от копипаста,
вы, наверное слой как "test" назвали, поэтому все и прошло )

[Dmitry Stasev]

Если интересно, вот код скрипта приведения к выпуклой форме:

Посмотрел, сложно для меня. Я бы просто добавил
st_convexhull( перед st_buffer(...

[SMOuk96]

Dmitry Stasev, Буфер действительно все "съел" ) хитрый способ

там в скрипте еще и ошибочка вкралась от копипаста

Даа, в шаблоне sql запроса)

Я бы хотел все-таки вот разобраться в sql запросе, знаю sql на уровне "create table", не могли бы немного по подробнее расписать?

[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  id, CLASS, class_mean, st_buffer(st_buffer(st_difference(test.geometry, CTE.geometry), -0.000001), 0.000001) as geometry FROM test
JOIN CTE 

1.With CTE AS()- обобщенное табличное выражение - Задает временно именованный результирующий набор.
Т.е., грубо говоря обертка на запрос с именем, у меня имя - CTE
В нашем случае, внутри конструкции отрабатывает запрос, к результатам которого можно обратится по имени (СТЕ) и связать (сджойнить) с основной таблицей test и все это в одном запросе (как раз для виртуального слоя).
SELECT ... FROM test
JOIN CTE
Подробнее о нем в инете почитаете хотя бы у МС https://docs.microsoft.com/ru-ru/sql/t- ... erver-2017

2. Вложенный в CTE запрос. Он формирует объект-пересечение (мультиполигон)
SELECT '1' AS grp, st_union(st_intersection(a.geometry, b.geometry)) as geometry
- '1' AS grp - эмулирую группу, создается поле с именем grp и значением - '1'. Нужно будет для агрегации (группировки) геометрий - (в сочетании с st_union() и GROUP BY grp )
- st_intersection(a.geometry, b.geometry) - пересечение "неодинаковых" геометрий.

FROM test a JOIN test b ON (
(a.id < b.id)and
- Неодинаковые потому что, хоть и тянутся с одной таблицы с разными алиасами (test a JOIN test b), но есть условие (a.id < b.id) т.е. 1и2, 1и3, 1и4, 1и5.... 2и3, 2и4, 2и5.....

(a.geometry is NOT NULL)and(b.geometry is NOT NULL)and
- геометрии не пустые

(st_intersects(a.geometry, b.geometry))
- и пересекаются
)
GROUP BY grp
-ну и условие группировки (а группа-то одна '1', т.е. и объект один выйдет) обеспечивает склеивание геометрий

ИТАК получили на выходе из CTE - объект-пересечение, который и вычитаем из объектов test
SELECT id, CLASS, class_mean, st_buffer(st_buffer(st_difference(test.geometry, CTE.geometry), -0.000001), 0.000001) as geometry FROM test
JOIN CTE

Вот как-то так

[SMOuk96]

Доброго времени суток уважаемые форумчане Пишу модуль на Python в среде OSGeo4W используя библиотеку OGR для работы с векторными данными. Есть shape-файл в котором слой полигонов, необходимо найти симметричную разность всех полигонов. То есть в конечном итоге должен получиться слой с полигонами без пересечений.
Кажется зашел в тупик пробуя находить симметричную разность каждого полигона с каждым...
Прилагаю shape-файл, код скрипта, картинку входных и выходных данных:

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

def symDif(inputShapefile, outShapefile = currentPath + "/temp/convexhull.shp")
    # Open base SHP File
    inputData = ogr.Open(inputShapefile, False)
    inputData2 = ogr.Open(inputShapefile, False)
    
    # Get once layer from SHP File
    inputLayer = inputData.GetLayer(0)
    inputLayer2 = inputData.GetLayer(0)
        
    # Dostup k atributam
    #inFeatDef = inputLayer.GetLayerDefn()

    # Create new SHP file
    drv = ogr.GetDriverByName("ESRI Shapefile")

    # Create exit dataset
    # outShapefile = currentPath + "/temp/states_convexhull.shp"

    # Remove output shapefile if it already exists
    if os.path.exists(outShapefile):
        drv.DeleteDataSource(outShapefile)

    outputData = drv.CreateDataSource(outShapefile)
    # Create exit layer
    outputLayer = outputData.CreateLayer(outShapefile, None, ogr.wkbPolygon) 
    
    # Create the feature
    featureDefn = outputLayer.GetLayerDefn()
    feature = ogr.Feature(featureDefn)
    
    i = 0
    j = 0
    
    feat = inputLayer.GetNextFeature()
    feat2 = inputLayer2.GetNextFeature()
    
    if feat is not None and feat2 is not None:
        print "not None"
        
    while feat is not None:
        count = 0
        feat = inputLayer.GetNextFeature()
        geom = feat.GetGeometryRef()
        while feat2 is not None:
            feat2 = inputLayer2.GetNextFeature()
            geom2 = feat2.GetGeometryRef()
            geom = geom.SymDifference(geom2)
            feat2 = inputLayer2.GetNextFeature()
            j+=1
        feature.SetGeometry(geom)
        outputLayer.CreateFeature(feature)
        feat = inputLayer.GetNextFeature()
        i+=1

    feature = None

    # Save and close DataSource
    inputData.Destroy()
    outputData.Destroy()

Моя идея в том чтобы находить SymDif сначала первого полигона со 2,3,4 и т. д., затем...
второго полигона с 1, 3,4 и т. д.
и каждый раз когда SymDif будет находиться перезаписывать полигон

Кажется идея верная, но реализация ужасная, помогите :
UPD 16:34 МСК (Код исправил на нужный)

[SMOuk96]

Если бы к полигонам можно было бы обращаться по индексам, то вообще задача была бы проста, но сейчас понял, что второй вложенный цикл исполняет только один раз. Да - есть функция ReturnReading() но она начнет чтение слоя с самого перового полигона...
Кстати при работе пользовался вот этой темой - http://gis-lab.info/qa/ogr-python.html

https://gis.stackexchange.com/questions ... -in-python
Нашел похожую проблему на ин. ресурсе, там предлагается решение библиотекой Fiona. Кто-то работал с ней?
Похожее решение как мне кажется:

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

with fiona.open('symdiff.shp','w',...) as output1:
   with fiona.open('normal.shp', 'w',...) as output2:
     for i in geoms1:
         for j in geoms2:
             if i.overlaps(j):
                # write the symdifference shapefile
             else:
                #write the normal shapefile

[SMOuk96]

Dmitry Stasev, Спасибо вам огромное - разобрался. Оказалось все не очень сложно) Теперь пытаюсь реализовать эту же задачу без qgis'овской библиотеки, только GDAL/OGR