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