Поиск пересечений двух shape (python)

Вопросы общего характера по ГИС и дистанционному зондированию, не связанные с конкретным ПО.
Elena92
Новоприбывший
Сообщения: 9
Зарегистрирован: 09 мар 2019, 23:18
Репутация: 0
Откуда: Мск

Поиск пересечений двух shape (python)

Сообщение Elena92 » 13 апр 2019, 19:39

Коллеги, добрый день!

Задача стоит простая, но никак не могу найти решения.
Исходные данные: 2 векторных слоя в формате shp:
- лесные кварталы (полигоны);
- контуры вырубок (линии).
Во вложении наглядное изображение двух наложенных друг на друга shp.

Задача состоит в следующем:
Необходимо написать python-скрипт, который выполняет некоторое преобразование.
На входе: номер вырубки (соответствующий одному из объектов в файле «контуры вырубок»).
На выходе: shp-файл, содержащий один полигональный объект – квартал, в котором находится указанная на входе вырубка.

По сути, необходимо найти пересечение обоих shp и выгрузить их поочерёдно. Возникла проблема с поиском пересечений. Я пыталась найти решение и на этом форуме в том числе, но, скорее всего, допускаю ошибки.
На данный момент есть вот это (опущу import):

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

# Драйвер
driver = ogr.GetDriverByName('ESRI Shapefile')

# WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# Открываем SHP-файл №1
# shp1 = ogr.Open("files/kv_AO/kv_AO.shp", 1)
os.chdir('files/kv_AO')
inputData = driver.Open('kv_AO.shp', 0)  # 0 - read, 1 - write
# Проверяем, что он не пустой
if inputData is None:
    print("Ошибка: файл №1 не найден.")
    sys.exit(1)

# Указываем рабочий слой, Записываем 0-ой слой в кортеж layer, обращение к слою по индексу
inputLayer = inputData.GetLayer(0)
if inputLayer is None:
    print("Ошибка чтения слоя в shp #1")
    sys.exit(1)

# Открываем SHP-файл №2
# shp2 = ogr.Open("files/test/test.shp", 1)
os.chdir('../test')
inputData2 = driver.Open('test.shp', 0)
# Проверяем, что он не пустой
if inputData2 is None:
    print("Ошибка: файл №2 не найден.")
    sys.exit(1)

# Указываем рабочий слой, Записываем 0-ой слой в кортеж layer, обращение к слою по индексу
inputLayer2 = inputData2.GetLayer(0)
inputLayer2.ResetReading()
if inputLayer2 is None:
    print("Ошибка чтения слоя в shp #2")
    sys.exit(1)

# подготавливаем все для записи результата
# выбираем драйвер
outputData = ".."
driverName = "ESRI Shapefile"
drv = ogr.GetDriverByName(driverName)
if drv is None:
    print("% драйвер недоступен." % driverName)
    sys.exit(1)
# создаем выходной набор данных
outputData = drv.CreateDataSource('SomeFilename.shp')
if outputData is None:
    print("Ошибка создания выходного файла.")
    sys.exit(1)
# и выходной слой
layerName = os.path.basename('SomeFilename').strip(".shp")
outputLayer = outputData.CreateLayer(layerName, None, geom_type=ogr.wkbLineString)
if outputLayer is None:
    print("Ошибка создания слоя.")
    sys.exit(1)
# Создаем схему (таблица атрибутов) слоя
# featureDefn = outputLayer.GetLayerDefn()
# feature = ogr.Feature(featureDefn)
feature = ogr.Feature(outputLayer.GetLayerDefn())

# начинаем просматривать объекты в исходном слое
inFeat = inputLayer2.GetNextFeature()
# и сразу же формируем таблицу атрибутов результирущего слоя
inFeatDef = inputLayer2.GetLayerDefn()

# делаем копию всех полей
for i in range(inFeatDef.GetFieldCount()):
    fieldDef = inFeatDef.GetFieldDefn(i)
    if outputLayer.CreateField(fieldDef) != 0:
        print("Ошибка создания поля %" % fieldDef.GetNameRef())
        sys.exit(1)
# добавлем поле для площади
fieldDef = ogr.FieldDefn("area_kv.m.", ogr.OFTReal)
fieldDef.SetPrecision(2)
if outputLayer.CreateField(fieldDef) != 0:
    print("Ошибка создания поля %" % fieldDef.GetNameRef())
    sys.exit(1)

# Пересечение
k = 1  # номер объекта
for i in inputLayer:
    count = 0  # счет кол-ва пересечений
    geom = i.GetGeometryRef()  # геометрия объекта
    for j in inputLayer2:
        geom2 = j.GetGeometryRef()  # геометрия объекта
        intersection = geom.Intersection(geom2) 
        # если пересечение найдено и это не сам объект
        if not intersection.IsEmpty():
            # записываем полученные пересечения в файл
            feature.SetGeometry(intersection)
            feature.SetField("area_kv.m.")
            outputLayer.CreateFeature(feature)
            outFeat = ogr.Feature(outputLayer.GetLayerDefn())
            count += 1
            # убираем мусор и переходим к следующему объекту
            outFeat.Destroy()
            inFeat = inputLayer2.GetNextFeature()
    inputLayer2.ResetReading()
    print("\nObj {} intersect {} objs\n".format(k, count))
    k += 1

# все объекты обработаны, закрываем наборы данных
inputData.Destroy()
outputData.Destroy()
Итоговый shp не заполняется, нет даже 1 объекта, хотя print исправно прописывает их нумерацию. Прошу вас подсказать, в чём моя ошибка? Может быть, есть другие пути решений?
Вложения
Снимок экрана от 2019-04-13 19-15-15.png
Снимок экрана от 2019-04-13 19-15-15.png (62.85 КБ) 11338 просмотров

trir
Гуру
Сообщения: 5278
Зарегистрирован: 09 апр 2010, 19:30
Репутация: 1014
Ваше звание: просто мимо прохожу
Откуда: Ё-бург

Re: Поиск пересечений двух shape (python)

Сообщение trir » 13 апр 2019, 20:28

(опущу import)
зря
вот по этому я всегда агитирую за БД и SQL - там всё гораздо проще!

gamm
Гуру
Сообщения: 4049
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 1050
Ваше звание: программист
Откуда: Казань

Re: Поиск пересечений двух shape (python)

Сообщение gamm » 13 апр 2019, 21:27

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

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

        # create feature
        feature = ogr.Feature(layer.GetLayerDefn())
        feature.SetField("ImageId", building['ImageId'])
        feature.SetField("BuildingId", building['BuildingId'])
        feature.SetGeometry(building['polyPix'])

        # Create the feature in the layer (geojson)
        layer.CreateFeature(feature)
        # Destroy the feature to free resources
        feature.Destroy()

Elena92
Новоприбывший
Сообщения: 9
Зарегистрирован: 09 мар 2019, 23:18
Репутация: 0
Откуда: Мск

Re: Поиск пересечений двух shape (python)

Сообщение Elena92 » 13 апр 2019, 22:39

outFit вылез из цикла, который уже не актуален. Удалила, спасибо. Порядок восстановила по схеме, но пока всё та же проблема.

gamm
Гуру
Сообщения: 4049
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 1050
Ваше звание: программист
Откуда: Казань

Re: Поиск пересечений двух shape (python)

Сообщение gamm » 14 апр 2019, 09:53

у вас в коде много подобного, обо что глаз спотыкается. Я конкретных функций не знаю, но
Elena92 писал(а):
13 апр 2019, 19:39
feature.SetField("area_kv.m.")
- ???, и скорее всего и
Elena92 писал(а):
13 апр 2019, 19:39
feature.SetGeometry(intersection)
тоже (не то вы выводите) ... вы бы попробовали для начала просто слой входной скопировать на выход в цикле.

Elena92
Новоприбывший
Сообщения: 9
Зарегистрирован: 09 мар 2019, 23:18
Репутация: 0
Откуда: Мск

Re: Поиск пересечений двух shape (python)

Сообщение Elena92 » 14 апр 2019, 19:22

gamm, добрый день!
1. С area я пыталась вывести площадь полигонов в новый слой (и, кстати, это получилось). Столбец не был переименован в связи с тем, что я плюнула на красоту и оставила его так х)
2. Да, я выводила входной слой. И он сохраняется в новом shp. Но при попытке сравнить 2 shape начинается какая-то вакханалия.

Возможно, я косячу на какой-нибудь ерунде.
Может быть, стоит использовать не intertools? Может быть, кто-то знает, с помощью каких фреймворков можно добиться лучшего результата?

trir
Гуру
Сообщения: 5278
Зарегистрирован: 09 апр 2010, 19:30
Репутация: 1014
Ваше звание: просто мимо прохожу
Откуда: Ё-бург

Re: Поиск пересечений двух shape (python)

Сообщение trir » 14 апр 2019, 20:04

Может быть, кто-то знает, с помощью каких фреймворков можно добиться лучшего результата?
БД гораздо удобней

STIntersection может вернуть коллекцию, а shp кажись её не понимает...

Elena92
Новоприбывший
Сообщения: 9
Зарегистрирован: 09 мар 2019, 23:18
Репутация: 0
Откуда: Мск

Re: Поиск пересечений двух shape (python)

Сообщение Elena92 » 14 апр 2019, 21:33

trir, добрый вечер.
К сожалению, ссылка не рабочая. Не совсем поняла, что Вы имеете ввиду. Какую БД Вы предлагаете использовать?

trir
Гуру
Сообщения: 5278
Зарегистрирован: 09 апр 2010, 19:30
Репутация: 1014
Ваше звание: просто мимо прохожу
Откуда: Ё-бург

Re: Поиск пересечений двух shape (python)

Сообщение trir » 14 апр 2019, 22:05


gamm
Гуру
Сообщения: 4049
Зарегистрирован: 15 окт 2010, 08:33
Репутация: 1050
Ваше звание: программист
Откуда: Казань

Re: Поиск пересечений двух shape (python)

Сообщение gamm » 15 апр 2019, 04:40

Elena92 писал(а):
14 апр 2019, 19:22
1. С area я пыталась вывести площадь полигонов в новый слой (и, кстати, это получилось). Столбец не был переименован в связи с тем, что я плюнула на красоту и оставила его так х)
2. Да, я выводила входной слой. И он сохраняется в новом shp. Но при попытке сравнить 2 shape начинается какая-то вакханалия.
Я в современных чудесах не силен, привык все руками делать, просто скажу, обо что глаз спотыкается:
1. Имя поля видно, присваиваемого значения не видно (может современные фреймворки его сами чудесным образом находят)
2. Вместо объекта из входного слоя выводится пересечение (там другая геометрия, а не исходный полигон). Т.е. проверили пересечение, если оно не пусто, то нужно выводить объект. А чтобы дублей не плодить, сначала накопить уникальные ID объектов, для которых есть пересечение (или поставить отметку таким объектам), а потом вторым циклом их вывести. Это, собственно, развертка SQL запроса, о котором говорит уважаемый trir

trir
Гуру
Сообщения: 5278
Зарегистрирован: 09 апр 2010, 19:30
Репутация: 1014
Ваше звание: просто мимо прохожу
Откуда: Ё-бург

Re: Поиск пересечений двух shape (python)

Сообщение trir » 15 апр 2019, 06:42


Dmitry Stasev
Участник
Сообщения: 67
Зарегистрирован: 13 мар 2018, 08:59
Репутация: 22
Откуда: MO

Re: Поиск пересечений двух shape (python)

Сообщение Dmitry Stasev » 15 апр 2019, 15:31

Да, код тяжеловат...
1. если нужен полигональный слой, тип надо другой geom_type
outputLayer = outputData.CreateLayer(layerName, None, geom_type=ogr.wkbLineString)
2. если открыли/создали датасет, желательно его освободить
3. srs создан, но не пользуется
4. Поля создаются, а данные не копируются...
и как у Вас получилось площадь записать без площади
feature.SetField("area_kv.m.")
?
5. Каша с циклами
Собирали из двух или более примеров?

попробуйте удалить всё ниже строки
# Создаем схему (таблица атрибутов) слоя
и добавьте этот код

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

inFeatDef = inputLayer2.GetLayerDefn()
# делаем копию всех полей
for i in range(inFeatDef.GetFieldCount()):
    fieldDef = inFeatDef.GetFieldDefn(i)
    if outputLayer.CreateField(fieldDef) != 0:
        print("Ошибка создания поля %" % fieldDef.GetNameRef())
        sys.exit(1)
# добавлем поле для площади
fieldDef = ogr.FieldDefn("area_kv.m.", ogr.OFTReal)
fieldDef.SetPrecision(2)
if outputLayer.CreateField(fieldDef) != 0:
    print("Ошибка создания поля %" % fieldDef.GetNameRef())
    sys.exit(1)

outFeatDef = outputLayer.GetLayerDefn()

# Пересечение
k = 1  # номер объекта
for i in inputLayer:
    count = 0  # счет кол-ва пересечений
    geom = i.GetGeometryRef()  # геометрия объекта
    for j in inputLayer2:
        geom2 = j.GetGeometryRef()  # геометрия объекта
        intersection = geom.Intersection(geom2) 
        # если пересечение найдено и это не сам объект
        if not intersection.IsEmpty():
            # записываем полученные пересечения в файл
            feature = ogr.Feature(outFeatDef)
            feature.SetGeometry(intersection)
            for f in range(0, outFeatDef.GetFieldCount()):
                feature.SetField(outFeatDef.GetFieldDefn(f).GetNameRef(), j.GetField(f))
            feature.SetField("area_kv.m.", intersection.Area())
            outputLayer.CreateFeature(feature)
            feature = None
            count += 1

    print("\nObj {} intersect {} objs\n".format(k, count))
    k += 1

# все объекты обработаны, закрываем наборы данных
inputData.Destroy()
inputData2.Destroy()
outputData.Destroy()
Последний раз редактировалось Dmitry Stasev 15 апр 2019, 18:07, всего редактировалось 1 раз.

Dmitry Stasev
Участник
Сообщения: 67
Зарегистрирован: 13 мар 2018, 08:59
Репутация: 22
Откуда: MO

Re: Поиск пересечений двух shape (python)

Сообщение Dmitry Stasev » 15 апр 2019, 15:57

Да, и если
На выходе: shp-файл, содержащий один полигональный объект – квартал, в котором находится указанная на входе вырубка.
То записывать надо не пересечение, а квартал.

Ivor
Завсегдатай
Сообщения: 345
Зарегистрирован: 11 дек 2006, 09:46
Репутация: 102
Откуда: Иркутск

Re: Поиск пересечений двух shape (python)

Сообщение Ivor » 16 апр 2019, 00:26

А почему не воспользоваться штатным инструментом Spatial Join (Пространственное соединение)? Получите список вырубок с присоединёнными id кварталов, а дальше элементарно выбрать нужный квартал под нужную вырубку. Или даже сразу список кварталов с вырубками, как подойти. Можно всё это в model builder запихнуть для автоматизации - получите инструмент, не привязанный к конкретным шейпам

Elena92
Новоприбывший
Сообщения: 9
Зарегистрирован: 09 мар 2019, 23:18
Репутация: 0
Откуда: Мск

Re: Поиск пересечений двух shape (python)

Сообщение Elena92 » 18 апр 2019, 20:54

Всем добрый день!

trir, интересная тема, я обязательно её поизучаю накануне :)

gamm, Вы абсолютно правы. Я переписала цикл и присвоила значения полям. После этого цикл заработал. С пересечением пока разбираюсь.

Dmitry Stasev, замечания в работу приняла, исправила недочёты) Код, к сожалению, пока не получилось применить, не заработал.

Ivor, К сожалению, мне не подходит готовое решение вроде qgis/arcgis/и т. д. Но, конечно, они решают мой вопрос)

Ответить

Вернуться в «Общие вопросы»

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

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