Программный поиск полигонов, которые> 90% перекрываются другим слоем векторного полигона, используя QGIS?

9

Примеры слоев

Я пытаюсь выяснить, как использовать Python для извлечения полигонов в одном векторе, которые перекрываются> 90% другим вектором. Затем я хотел бы иметь вектор / карту, которая будет показывать только эти полигоны. Пример изображения показывает мои слои. Я хочу, чтобы все серые многоугольники были> 90% красного цвета.

Мне нужно сделать все это через Python (или аналогичным образом автоматизированные методы). У меня есть ~ 1000 карт для обработки таким же образом.

dnormous
источник
Вы хотите сделать наложение 'union' (см. Infogeoblog.wordpress.com/2013/01/08/geo-processing-in-qgis для некоторых основ), а затем для каждого исходного полигона вычислить статистику 'in' и статистику 'out' gis.stackexchange.com/questions/43037/… чтобы определить процент оверлея ... подсказка: вам нужно провести измерение площади gis.stackexchange.com/questions/23355/…
Майкл Стимсон
Спасибо за советы. Это тот же подход, который я только что попробовал. Я могу сделать объединение через консоль Python достаточно легко. Уже добавлены в область значения атрибутов. Это следующий шаг, в котором я не уверен. Как я могу использовать python для вычисления статистики «in» и «out», чтобы я мог идентифицировать / выбрать / обрезать и т. Д.> 90% полигонов?
dnormous
Я думаю, что это возможно без питона. Вам нужен Python или решение с виртуальными слоями подходит для вас?
Пьерма
Области «in» будут иметь атрибуты из обоих полигонов, области «out» имеют атрибуты только из одного набора полигонов. Получите оба набора статистики области и присоединитесь к исходным многоугольникам, добавьте поле для «in», «out» и покрытия, вычислите значения для «in» и «out» из суммы площадей, затем разделите «in» на исходная область (или «в» + «вне») для расчета процентов.
Майкл Стимсон
1
Пьерма - мне просто нужен автоматизированный метод, чтобы найти полигоны.
августа

Ответы:

3

Следующий код работает в моей консоли Python QGIS. Он создает слой памяти с полигонами, которые> 90% перекрываются красными областями.

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for polygon_intersects
feats_lyr1 = [ feat for feat in layers[0].getFeatures() ]

#for xwRcl
feats_lyr2 = [ feat for feat in layers[1].getFeatures() ]

selected_feats = []

for i, feat1 in enumerate(feats_lyr1):
    area1 = 0
    area2 = 0
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area = feat1.geometry().intersection(feat2.geometry()).area()
            print i, j, area, feat2.attribute('class')
            if feat2.attribute('class') == 1:
                area1 += area
            else:
                area2 += area
    crit = area1/(area1 + area2)
    print crit
    if crit > 0.9:
        selected_feats.append(feat1)

epsg = layers[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "mem_layer",
                           "memory")

prov = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Я опробовал код с этими двумя векторными слоями:

введите описание изображения здесь

После запуска кода на консоли Python QGIS для подтверждения результатов были напечатаны индексы i, j задействованных объектов, областей пересечения, атрибута поля в polygons_intersects (1 для красных областей и 2 для серых областей) и критерий перекрытия ,

0 0 9454207.56892 1
0 1 17429206.7906 2
0 2 10326705.2376 2
0 4 40775341.6814 1
0 5 26342803.0964 2
0 7 11875753.3216 2
0.432253120382
1 6 1198411.02558 2
1 7 1545489.96614 2
1 10 27511427.9909 1
0.90930850584
2 7 750262.940888 2
2 8 12012343.5859 1
0.941213972294
3 6 23321277.5158 2
0.0

Созданный слой памяти (зеленые объекты) можно наблюдать на следующем изображении. Это было, как и ожидалось.

введите описание изображения здесь

xunilk
источник
6

Здесь решение, которое не требует Python.

Добавьте новый виртуальный слой с запросом вроде:

WITH r AS (
SELECT 
    Basins800.rowid AS idGray, 
    area(Basins800.geometry) AS areaGray, 
    area(Intersection(Basins800.geometry, Severity.geometry)) AS aeraInter, 
    Basins800.geometry AS geomGray 
  FROM Basins800, Severity
)

SELECT *, areaInterSum/areaGray  AS overlap , geomGray 
    FROM (
        SELECT 
           idGray, 
           areaGray, 
           sum(areaInter) AS areaInterSum, 
           geomGray 
        FROM r 
        GROUP BY idGray) 
     WHERE areaInterSum/areaGray > 0.9

С :

  • Basins800 в качестве вашего слоя вы хотите фильтровать с серыми полигонами

  • Серьезность: ваш красный слой перекрывается.

Результатом будет новый слой, в котором только все серые полигоны> 90% будут перекрыты красными полигонами, а новое поле будет содержать процент перекрытия.

введите описание изображения здесь

Надеюсь, это работает. Я могу добавить больше деталей по запросу, если это необходимо.

Примечание: ваши данные содержат очень маленькие полигоны (исходя из обработки растра и соответствующие растровому пикселю (на рисунке мы видим 4 полигона, но есть еще 25 небольших полигонов). Это делает запрос очень медленным для выполнения (функция пересечения) генерирует один элемент для каждой пары элементов из двух слоев).

Pierma
источник
Я получаю сообщение об ошибке при выполнении запроса с помощью кнопки «создать виртуальный слой». "Ошибка выполнения запроса при CREATE TEMP VIEW _tview AS WITH r AS (" ... остальная часть кода здесь ... с последующим: "1 - рядом" WITH ": синтаксическая ошибка" Я довольно новичок в QGIS. Могу ли я создать этот виртуальный слой программно? Спасибо за вашу помощь!
dnormous
Вот ссылка для скачивания шейп-файлов: ссылка
dnormous
Извините, плохая копия между серым и серым (извините за мой примерный английский). Я отредактировал запрос. Это должно работать сейчас. Почему вы хотите создать слой программно? Преимущества виртуального слоя в том, что он неразрушающий, и если вы редактируете свои данные (серый или красный многоугольники), виртуальный слой автоматически обновляется.
Пьерма
Это только один маленький кусочек процесса. У меня есть ~ 1000 таких карт, поэтому автоматизация процесса будет чрезвычайно полезна.
dnormous
Я все еще получаю ту же ошибку -> "1 - рядом с" WITH ": синтаксическая ошибка". Я включил локальные имена для каждого слоя в GrayLayer и RedLayer. Это местное имя, что я должен использовать? то есть: серый слой помечен как «Basins_800», поэтому у меня есть код, такой как «Basins_800.geometry»
dnormous
2

Увидев ссылку на шейп- файлы Severity и Basins800 , я смог понять необходимый геопроцесс. Я изменил код в:

Программный поиск полигонов, которые> 90% перекрываются другим слоем векторного полигона, используя QGIS?

для получения этого:

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

#for Severity
feats_lyr1 = [ feat for feat in layers[0].getFeatures() ]

#for Basins800
feats_lyr2 = [ feat for feat in layers[1].getFeatures() ]

selected_feats = []

print "processing..."

for i, feat1 in enumerate(feats_lyr1):
    for j, feat2 in enumerate(feats_lyr2):
        if feat1.geometry().intersects(feat2.geometry()):
            area1 = feat1.geometry().intersection(feat2.geometry()).area()
            area2 = feat1.geometry().area()
            print i, j, area1, area2
    crit = area1/area2
    print crit
    if crit > 0.9:
        selected_feats.append(feat1)

epsg = layers[0].crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           "mem_layer",
                           "memory")

prov = mem_layer.dataProvider()

for i, feat in enumerate(selected_feats):
    feat.setAttributes([i])

prov.addFeatures(selected_feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Через несколько минут после запуска кода с этими шейп-файлами на Python Console в QGIS я получил результат, аналогичный Pierma ; где слой памяти имел 31 объект (разные из 29 полученных им полигонов).

введите описание изображения здесь

Я не собираюсь отлаживать результаты, потому что существует 1901 * 3528 = 6706728 взаимодействий для функций. Тем не менее, код выглядит многообещающе.

xunilk
источник