Река-притоки

[sergiks]

Здравствуйте.
Есть линейный слой гидрография (шейп). Нужно для реки найти все притоки и притоки притоков ...
Первое, что пришло в голову: пространственный выбор "касание", который с каждым разом добавляет притоки. Работает медленно.
Может есть что более быстрое? Сетевой анализ, например. Но не пойму, как от построенного графа перейти к моей реке и ее притокам.

[Александр Мурый]

Если есть ЦМР, то логичнее всего построить бассейн и потом выделить все реки, которые внутри.

[trir]

GeoSPARQL

[Александр Мурый]

Есть такой вариант:
— строим буферные зоны фиксированной ширины по всем рекам (обязательно с объединением буферов);
— разбиваем составные объекты (буферы получаются мультиполигоном, а нужен набор полигонов);
— выделяем буферную зону нужной реки и делаем "извлечь по пространственному положению"; выбрать из слоя с реками, слой пересечения — буферные зоны, "находится внутри".

[Elf]

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

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

def inflower(a):

    for selected in a:

        for feature in layer.getFeatures():

            if feature.geometry().intersects(selected.geometry().boundingBox()):

                inflow.append(feature.id())

    layer.select(inflow)



inflow = []

layer = iface.activeLayer()

selection = layer.selectedFeatures()

inflower(selection)

Update 0:
В таком виде запускается бесконечная рекурсия, которая все же остановится через пару тысяч итераций (у меня QGIS зависает на 2-3 минуты). Цель достигнута (все притоки выделены), но способ так себе. Попробую довести до ума.

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

inflow = []



def inflower():

    layer = iface.activeLayer()

    selection = layer.selectedFeatures()

    for selected in selection:

        for feature in layer.getFeatures():

            if feature.geometry().intersects(selected.geometry().boundingBox()):

                inflow.append(feature.id())

    layer.select(inflow)

    inflower()



inflower()

[sergiks]

В таком виде запускается бесконечная рекурсия, которая все же остановится через пару тысяч итераций (у меня QGIS зависает на 2-3 минуты). Цель достигнута (все притоки выделены), но способ так себе. Попробую довести до ума.

Я тоже написал подобный скрипт. Но работает крайне медленно. Быстрее получается при использовании алгоритмов геообработки:

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

# coding: utf8

import processing

n1=1

layer = iface.activeLayer()

while True:

    riv=processing.runalg('qgis:saveselectedfeatures', layer, None)

    processing.runalg('qgis:selectbylocation', layer, riv['OUTPUT_LAYER'], ['touches'] ,1)

    fs=layer.selectedFeatures()

    print u'Выделено объектов:', len(fs)

    n2=len(fs)

    if n2==n1:

        break        

    n1=n2

Скрипт выделяет все притоки любого порядка.

[Александр Мурый]

Скрипт, написанный *sergiks*, работает лучше, если сначала прогнать слой рек через модуль GRASS <v.clean> с параметром "break" с небольшим произвольным порогом (threshold). Затем уже в полученном слое выделяем любой приток (вне зависимости от порядка) из интересующей нас реки и запускаем скрипт; выделяются вся река с притоками.
Основное и очевидное условие: все линии должны соприкасаться друг с другом.

[Elf]

Незавершенный скрипт не давал мне покоя. После консультации с ГИС-специалистом все же удалось довести дело до конца.

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

#!/usr/bin/env python

# -*- coding: utf-8 -*-

'''

Cкрипт выделяет притоки всех порядков выделенной реки.

The script selects all inflows of selected river.



Недостаток: другая река или приток другой реки также будут выделены, если они 

находятся слишком близко (пересекают минимальный прямоугольник, в который 

вписана изначально выделенная река и (или) ее приток).

'''

from qgis.core import *

from qgis.utils import *



layer = iface.activeLayer()

geom_ids = []



def inflower(line, set):

    for item in set.copy():

        if item.geometry().intersects(line.geometry().boundingBox()):

            geom_ids.append(item.id())

            set.remove(item)

            inflower(item, set)

    layer.select(geom_ids)



for selected in layer.selectedFeatures():

    river = selected



features = set()

for feature in layer.getFeatures():

    if river.id() != feature.id():

        features.add(feature)

        

inflower(river, features)



features.clear()

geom_ids = []

[sergiks]

Недостаток: другая река или приток другой реки также будут выделены, если они
находятся слишком близко (пересекают минимальный прямоугольник, в который
вписана изначально выделенная река и (или) ее приток).

А если проверять не пересечение с прямоугольником, а само касание, т.е. вместо

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

if item.geometry().intersects(line.geometry().boundingBox()):

использовать:

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

if item.geometry().touches(line.geometry()):

[Elf]

Пробовал и так. Почему-то ничего не выделялось.