Вложенный цикл по слою Shape-файла (перебор объектов)

Вопросы по нескольким пакетам сразу, или вопросы, которые непонятно к какой ГИС отнести
Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 17 июн 2018, 21:48

Доброго времени суток уважаемые форумчане! Столкнулся с проблемой перебора всех объектов слоя shape-файла... имеется следующий код:

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

    inputData = ogr.Open(inputShapefile, False)
    inputData2 = ogr.Open(inputShapefile, False)
    inputLayer = inputData.GetLayer(0)
    inputLayer2 = inputData2.GetLayer(0)
    k = 0
    for i in inputLayer:
        count = 0
        geom = i.GetGeometryRef()
        for j in inputLayer2:
            geom2 = j.GetGeometryRef()
            if geom.Intersection(geom2):
                count+=1
        print "Obj {} intersect {} objs".format(k, count)
        k +=1
Считаю количество пересечений объекта с другими объектами для примера, но проблема в другом. Вложенный цикл исполняется только 1 раз.

Может кто-то знает иной способ, как перебрать объекты каждый с каждым в одном слое?

Результат работы скрипта на shape-файле из 40 объектов во вложении (считает неправильно даже для первого объекта)
Вложения
WV2_B45_20160625_C_NDVI_F165_S_Sieve_Convexhull-9.shp
(12.58 КБ) 753 скачивания
1.jpg
1.jpg (25.6 КБ) 10310 просмотров

doujin
Активный участник
Сообщения: 163
Зарегистрирован: 28 июн 2012, 01:02
Репутация: 84
Откуда: Vladivostok

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение doujin » 18 июн 2018, 05:45

Смотрите в документацию: Iterate over Features.
А именно:
You must call ResetReading if you want to start iterating over the layer again.
Перебор объектов идет в итераторе. Дойдя до конца в первый раз, итератор там в конце и остается, поэтому у вас последующие переборы не выполняются. Перед каждым новым проходом итератор надо обновить, вызвав ResetReading. Смотрите пример в документации по ссылке.

Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 18 июн 2018, 13:55

doujin, Спасибо! Добавил ResetReading после вложенного цикла, все работает, спасибо)
Количество объектов пересечения считает верно, однако теперь записываю пересечения в выходной файл и столкнулся с другой проблемой. А именно, что пересечения записываются дважды или трижды..

Простой пример того как работает алгоритм:
Есть 2 объекта, которые пересекаются. Обходим их: Ищем пересечение Объект[1] с Объект[2] и Объект[2] с Объект[1].

Как исключить поиск пересечения таких пар? Есть идея записывать пару индексов объектов в список и затем каждый раз при поиске пересечения искать реверс текущей пары. Пример:

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

list[[0,1],[0,2]] #сюда записываем пары
list2[2,0]
for i in list:
	if i == list2:
		#тогда эту пару исключаем
А вот код, который есть сейчас:

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

f1 = 0 #индекс 1ого объекта
    k = 1 #индекс для наглядности (номер объекта)
    for i in inputLayer:
        count = 0 # переменная для счета кол-ва пересечений
        geom = i.GetGeometryRef() # геометрия объекта
        f2 = 0 #индекс 2ого объекта
        for j in inputLayer2:
            geom2 = j.GetGeometryRef() # геометрия объекта
            inters = geom.Intersection(geom2) #пересечение 1ого объекта со 2ым
            #если пересечение найдено и это не сам объект
            if inters.GetGeometryType() == ogr.wkbPolygon and f1 != f2:
            	#записываем полученные пересечения в файл
                feature.SetGeometry(inters)
                outputLayer.CreateFeature(feature)
                count+=1
            f2+=1 #увеличиваем индекс 2ого объекта
        inputLayer2.ResetReading()
        print "\nObj {} intersect {} objs\n".format(k, count)
        k +=1
        f1+=1 #увеличиваем индекс 1ого объекта
На выходе получается слой на котором по несколько объектов-пересечения друг на друге

Искал в сети решение по лучше, но не могу найти. Уверен есть метод проще и легче - исключить такие пары
Вложения
2.jpg
2.jpg (120.92 КБ) 10238 просмотров

doujin
Активный участник
Сообщения: 163
Зарегистрирован: 28 июн 2012, 01:02
Репутация: 84
Откуда: Vladivostok

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение doujin » 18 июн 2018, 16:34

Не совсем понимаю, почему у вас дублируются пересечения... inputLayer и inputLayer2 это один и тот же файл? То есть задача найти все пересечения объектов внутри одного шейп-файла?
Если так, то мое предложение такое:

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

import itertools

layer = tuple(inputData.GetLayer(0))
intersections_count = 0
for i, j in itertools.combinations(range(len(layer)), 2):
    geom_i = layer[i].GetGeometryRef()
    geom_j = layer[j].GetGeometryRef()
    intersection = geom_i.Intersection(geom_j)
    if intersection:
        feature.SetGeometry(intersection)
        outputLayer.CreateFeature(feature)
        intersections_count += 1
Т.е не делать лишних проходов. Составить все комбинации индексов без повторов. Потом дергаем объекты по их индексам за один проход. Если пересечение чего-нибудь вернуло, увеличиваем счетчик и пишем в выходной файл.
Код я не тестировал, просто хотел показать идею, но, вроде, должно сработать.

Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 18 июн 2018, 19:27

doujin, циклы itertools куда удобнее и функциональнее стандартных питоновских, не пробовал их еще использовать) Спасибо) Работает, немного изменил, чтобы не создавало объекты без геометрии :)

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

    id = 0
    for i, j in itertools.combinations(range(len(layer)), 2):
        geom_i = layer[i].GetGeometryRef()
        geom_j = layer[j].GetGeometryRef()
        intersection = geom_i.Intersection(geom_j)
        Check_Intersection = intersection.GetGeometryType()
        if Check_Intersection == ogr.wkbPolygon:
            feature.SetGeometry(intersection)
            outputLayer.CreateFeature(feature)
            feature.SetField("id", id)
            id+=1
Одна непонятность еще - в таблице атрибутов почему-то со второго элемента id только начинает заполнять.

А вообще моя цель найти не пересечения, а симметричную разность. И вот думаю что выйдет же все отлично, если из каждого объекта исходного файла таким же способом отнимать (методом Difference) полученные объекты пересечения.
Спойлер

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

    out_layer = tuple(outputLayer)
    for i, j in itertools.combinations(range(len(out_layer)), 2):
        geom_i = layer[i].GetGeometryRef()
        geom_j = out_layer[j].GetGeometryRef()
        geom_i = geom_i.Difference(geom_j)
        Check_Difference = geom_i.GetGeometryType()
        #??? как не записывать все разности по отдельности
        if Check_Difference == ogr.wkbPolygon:
            feature.SetGeometry(geom_i)
            outputLayer.CreateFeature(feature)
Наверное можно найти Difference и в предыдущем цикле, но не знаю. А здесь получается так, что очень много разностей друг на друга накладываются.

Еще читал что есть метод SymmetricDifference, пробовал его применять, но не могу додумать каким образом не накладывать результаты функции друг на друга :cry:

P.S.: Начал с пересечений (Intersection) т.к. проще разобраться с этим)
Вложения
3.jpg
3.jpg (47.59 КБ) 10203 просмотра

doujin
Активный участник
Сообщения: 163
Зарегистрирован: 28 июн 2012, 01:02
Репутация: 84
Откуда: Vladivostok

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение doujin » 19 июн 2018, 04:44

Я полагал, что Intersection возвращает None в случае, если пересечение не найдено. Но на самом деле он возвращает объект геометрии, но пустой. Проверить его лучше методом IsEmpty, а в вашем коде проверка какая-то слишком замороченная.

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

intersection = geom_i.Intersection(geom_j)
if not intersection.IsEmpty():
    # запись в выходной файл
В отношении симметричной разности.
Согласно документации есть специальный метод SymDifference. В плане алгоритма нам надо:
1. взять первые два объекта из слоя
2. получить их разность
3. взять разность и следующий объект
4. получить уже их разность
5. повторить 3. и 4. пока объекты не закончатся
6. последняя разность и будет искомым результатом

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

layer = inputData.GetLayer(0)
first_geom = next(layer).GetGeometryRef()
second_geom = next(layer).GetGeometryRef()
diff = first_geom.SymDifference(second_geom)
for goem in layer:
    # удобно, что layer это итератор, можно продолжить брать из него объекты
    diff = diff.SymDifference(geom.GetGeometryRef())
# В итоге в diff получаем _мультиполигон_, который и представляет собой
# конечный результат симметричной разности. Тут его можно записать в
# выходной файл
Вроде все просто, но могут быть какие-то нюансы.
По поводу пустого id.
Вы сначала в feature устанавливаете геометрию, потом записываете feature в файл, потом для feature устанавливаете значение поля id... но стоп! Вы же уже записали его в файл!
В следующей итерации цикла у вас feature уже с id, который остался с прошлого раза. Вы устанавливаете ему геометрию и пишете в файл... В общем, сначала установите все настройки для feature, а только потом создавайте из него объект в файле.

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

feature.SetGeometry(intersection)
feature.SetField("id", id)
outputLayer.CreateFeature(feature)
Не уверен, что такое подробное объяснение требовалось. Мне просто показалось это забавным.

Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 19 июн 2018, 07:17

doujin, пробовал реализовать нечто подобное, но более сложным методом через два цикла :)
Скрипт вылетал без ошибки. Ваш код красивее конечно, но к сожалению тоже вылетает :o

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

    layer = inputData.GetLayer(0)
    print "line 1"
    #first_geom = next(layer).GetGeometryRef()
    #second_geom = next(layer).GetGeometryRef()
    first_geom = layer.GetNextFeature().GetGeometryRef()
    print "line 2"
    second_geom = layer.GetNextFeature().GetGeometryRef()
    print "line 3"
    diff = first_geom.SymDifference(second_geom)
    print "line 4"
    for geom in layer:
        diff = diff.SymDifference(geom.GetGeometryRef())
    feature.SetGeometry(diff)
    outputLayer.CreateFeature(diff)
Пробовал по-другому получать следующий объект из слоя, но нет, проблема не в этом. А скорее всего мне кажется на этапе нахождения "diff". Не знаю как дебажить python код в OSGeo4W, поэтому понатыкал принтов после каждой строчки) вот что выдает:
4.jpg
4.jpg (31.13 КБ) 10143 просмотра

Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 19 июн 2018, 13:11

Вот полностью готовый скрипт, который ищет все пересечений полигонов в shape-файле и сохраняет в выходной shape-файл.

Входные данные: input.shp
Выходные данные: output.shp

Скрипт запускается из OSGeo4W и использует библиотеку osgeo.ogr, а также:
  • Os
  • Math
  • Sys
  • itertools
Спойлер

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

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

from osgeo import ogr
import sys
import os
import math
import itertools

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

def intersections(inputShapefile, outShapefile = currentPath + "/temp/intersections.shp"):
    # Открываем SHP-файл
    inputData = ogr.Open(inputShapefile, False)
    
    # Проверяем, что он не пустой
    if inputData is None:
        print("ERROR: open failed")
        sys.exit(1)
        
    # Записываем 0-ой слой в кортеж layer
    layer = tuple(inputData.GetLayer(0))
    
    # Проверяем, что он не пустой
    if layer is None:
        print("ERROR: reading file")
        sys.exit(1)
        
    # Создаем новый shp-файл
    drv = ogr.GetDriverByName("ESRI Shapefile")
    
    # Удаляем shp-файл, если он уже создан
    if os.path.exists(outShapefile):
        drv.DeleteDataSource(outShapefile)
        
    # Создаем DataSet
    outputData = drv.CreateDataSource(outShapefile)
    
    # Создаем outputLayer - выходной слой
    try:
        outputLayer = outputData.CreateLayer(outShapefile, None, ogr.wkbPolygon)
    except AttributeError:
        print ("Log:\tCan not delete old output.shp")
        sys.exit(1)
    
    # Добавляем Поле "id" типа Integer
    idField = ogr.FieldDefn("id", ogr.OFTInteger)
    outputLayer.CreateField(idField)
    
    # Создаем схему (таблица атрибутов) слоя
    featureDefn = outputLayer.GetLayerDefn()
    feature = ogr.Feature(featureDefn)
    
    id = 0 # переменная для записи id объекта
    
    # интератор без повторяющихся элементов длиной 2 из layer
    for i, j in itertools.combinations(range(len(layer)), 2):
    
        # записываем геометрию i, j объектов в geom_i и geom_j
        geom_i = layer[i].GetGeometryRef()
        geom_j = layer[j].GetGeometryRef()
        
        # записываем в intersection пересечение geom_i и geom_j геометрий
        intersection = geom_i.Intersection(geom_j)
        
        # если intersection не пустое (т.е. пересечение geom_i и geom_j найдено)
        if not intersection.IsEmpty():
        
            # В значение поля текущего объекта записываем id
            feature.SetField("id", id)
            
            # Записываем геометрию пересечения в feature
            feature.SetGeometry(intersection)
            
            # записываем feature в выходной файл
            outputLayer.CreateFeature(feature)
            # увеличиваем id на единицу
            id+=1 
            
    # очищаем feature         
    feature = None

    # Сохраняем и закрываем DataSet
    inputData.Destroy()
    outputData.Destroy()
    
if __name__ == "__main__":
    if len(sys.argv) != 3:
        print ("\n\nERROR PARAMETERS. EXAMPLE: \nintersections.py C:/NDVI/ndvi/input.shp C:/NDVI/ndvi/output.shp\n")
        sys.exit(1)
    else:
        inputSHPPath = sys.argv[1]
        outputSHPPath = sys.argv[2]
        intersections(inputSHPPath, outputSHPPath)
        print ("\n\tReady!\n\n\tSee result in {}\output.shp".format(os.path.dirname(os.path.abspath(__file__))))
Результат работы скрипта:
Синие полигоны: исходный слой
Зеленые полигоны: выходной слой
Спойлер
5.jpg
5.jpg (78.85 КБ) 10106 просмотров
Надеюсь кому-то это пригодиться, по крайней мере поможет разобраться в OGR :wink:
И найдем решение для Симметричной разности :cry:

Аватара пользователя
rhot
Гуру
Сообщения: 1727
Зарегистрирован: 25 янв 2011, 17:50
Репутация: 194
Ваше звание: доктор
Откуда: Архангельск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение rhot » 19 июн 2018, 13:21

SMOuk96 писал(а):
19 июн 2018, 07:17
Не знаю как дебажить python код в OSGeo4W
Также как и везде:

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

try:
    func(x) # кусок кода, который надо отладить
exception Exception as e:
    print(e)
___________(¯`·.¸(¯`·.¸ Scientia potentia est _/ {SILVA}:::{FOSS}:::{GIS} \_ Знание сила ¸.·´¯)¸.·´¯)___________

Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 19 июн 2018, 13:34

rhot, спасибо, попробовал использовать данную конструкцию, но к сожалению все равно просто также зависает и «Прекращена работа python.exe», никакой ошибки не пишет :shock:

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

    layer = inputData.GetLayer(0)
    first_geom = layer.GetNextFeature().GetGeometryRef()
    second_geom = layer.GetNextFeature().GetGeometryRef()
    try:
        diff = first_geom.SymDifference(second_geom)
    except Exception as e:
        print(e)
    for geom in layer:
        diff = diff.SymDifference(geom.GetGeometryRef())
    feature.SetGeometry(diff)
    outputLayer.CreateFeature(diff)
Не знаете в чем может быть ошибка :?:

Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 19 июн 2018, 13:38

doujin, с id очевидна ошибка, сам понял насколько это было глупо сначала записывать объект в файл, а потом устанавливать значение полю id этому объекта :)

Аватара пользователя
rhot
Гуру
Сообщения: 1727
Зарегистрирован: 25 янв 2011, 17:50
Репутация: 194
Ваше звание: доктор
Откуда: Архангельск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение rhot » 19 июн 2018, 14:13

А если явно попросить выводить

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

ogr.UseExceptions()
___________(¯`·.¸(¯`·.¸ Scientia potentia est _/ {SILVA}:::{FOSS}:::{GIS} \_ Знание сила ¸.·´¯)¸.·´¯)___________

Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 19 июн 2018, 15:10

rhot писал(а):
19 июн 2018, 14:13
А если явно попросить выводить
Тоже ничего не выдает. Решил проблему случайно.

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

first_geom = layer.GetNextFeature().GetGeometryRef()
Не знаю с чем это связано, но получать следующий объект и сразу же получать его геометрию нельзя :o

Рабочий код, который выдает мультиполигон симметричной разности:

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

    inputData = ogr.Open(inputShapefile, False)
    layer = inputData.GetLayer(0)
    first = next(layer)
    second = next(layer)
    first_geom = first.GetGeometryRef()
    second_geom = second.GetGeometryRef()
    diff = first_geom.SymDifference(second_geom)
    for geom in layer:
        diff = diff.SymDifference(geom.GetGeometryRef())
    feature.SetGeometry(diff)
    outputLayer.CreateFeature(feature)
Теперь только осталось решить проблему того, как разделить мульти полигон на отдельные полигоны... Знаю, что в QGIS есть в модуле fTools функция "Разбить составную геометрию (multi_to_single)"
А в OGR есть что-то такое же?

doujin
Активный участник
Сообщения: 163
Зарегистрирован: 28 июн 2012, 01:02
Репутация: 84
Откуда: Vladivostok

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение doujin » 19 июн 2018, 15:58

Все еще проще. Просто перебрать отдельные части в цикле.

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

for part in diff:
    feature.SetGeometry(part)
    outputLayer.CreateFeature(feature)

Аватара пользователя
SMOuk96
Интересующийся
Сообщения: 29
Зарегистрирован: 30 июн 2017, 17:07
Репутация: 2
Откуда: Красноярск

Re: Вложенный цикл по слою Shape-файла (перебор объектов)

Сообщение SMOuk96 » 19 июн 2018, 18:41

doujin писал(а):
19 июн 2018, 15:58
Все еще проще. Просто перебрать отдельные части в цикле.
Ох, отлично. То есть мультиполигон как список, можно обращаться к каждому объекту внутри него :)

Ответить

Вернуться в «Общий - ПО»

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 6 гостей