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

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

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

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

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

Задача стоит простая, но никак не могу найти решения.
Исходные данные: 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 КБ) 5469 просмотров

wladfm
Участник
Сообщения: 69
Зарегистрирован: 04 июл 2016, 16:02
Репутация: 8

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

Сообщение wladfm » 17 апр 2019, 16:49

Elena92 писал(а):
13 апр 2019, 19:36
По сути, необходимо найти пересечение обоих shp и выгрузить их поочерёдно
Не совсем понятна задача к сожалению. Вы данную задачу можете решить с помощью ГИС-приложения?

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

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

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

wladfm, добрый день!
Да, могу.
P.S. Тема была перемещена администраторами и теперь, после перемещения темы, все ответы и обсуждение остались в старом разделе: viewtopic.php?f=1&t=25359

Ответить

Вернуться в «Я новичок!»

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

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