Как замкнуть соседние линии?

[Denis Rykov]

Предположим, что у нас есть несколько линий. И нам необходимо соединить их концы в случае, если они лежат в пределах определенного расстояния. Я попробовал использовать для этого v.clean:

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

v.clean input=input_line@PERMANENT output=output_line@PERMANENT type=line tool=snap thres=200

но полученный результат не устраивает по той причине, что он меняет ориентацию исходных геометрий (см. вложение). Можно ли каким-то образом замкнуть линии путём их продолжения до места пересечения?

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

Вот так сходу этого сделать нельзя. Надо подумать

[Denis Rykov]

Может быть стоит перенести тему в Общие? Просто я подумал, что в GRASS точно должен быть нужный инструмент.

[trir]

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

FUNCTION [dbo].[LS_AddPoint]
(
	@Ls geometry,
	@np geometry
)
RETURNS geometry
AS
BEGIN
  Declare @I Integer =1, @ICount Integer;
  Declare @nStr nvarchar(max);
  Declare @crP geometry;
  Set @ICount = @Ls.STNumPoints() + 1;
  Set @nStr = 'LineString(';
  while @I < @ICount 
   begin
	Set @crP = @Ls.STPointN(@I);
	Set @nStr = Concat(@nStr, SUBSTRING(@crP.ToString(), 8, Len(@crP.ToString())-8), ',');
	--Concat(@nStr, ltrim(STR(@crP.STX, 7, 11)), ' ', ltrim(STR(@crP.STY, 99, 11)), ',');
	set @I = @I + 1;
   end;
  Set @nStr = Concat(@nStr, SUBSTRING(@nP.ToString(), 8, Len(@nP.ToString())-8), ')');
  --Concat(@nStr, ltrim(STR(@nP.STX, 99, 11)), ' ', ltrim(STR(@nP.STY, 99, 11)), ')');
  Return geometry::STGeomFromText(@nStr, 0);
END

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

FUNCTION [dbo].[LS_Union]
(
@sls geometry,
@els geometry
)
RETURNS geometry
AS
BEGIN
	-- Declare the return variable here
	DECLARE @i int = 1, @iCount int;
	DECLARE @wP geometry;
	--
	Set @ICount = @eLs.STNumPoints();
	while @I < @ICount + 1
	 begin
	   Set @wP = @eLs.STPointN(@I);
	   Set @sls = dbo.LS_AddPoint(@sls, @wP);
	   set @I = @I + 1;
	 end;
	--
	RETURN @sls;
END

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

Может быть стоит перенести тему в Общие? Просто я подумал, что в GRASS точно должен быть нужный инструмент.

Пусть это висит здесь, в "Общем" лучше отдельную тему создать. Кстати, очень распространённое заблуждение насчёт GRASS.

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

А если надо по-быстрому просто соединить концы линий, то вот велосипедик:

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

# получаем вершины линий
v.to.points -v in=lns out=lns_vert
# делаем текстовый дамп вершин, фильтруем координаты и загоняем их назад в единую линию
v.out.ascii in=lns_vert | cut -d'|' -f1,2 | v.in.lines in=- out=lns_new fs='|'

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

trir, как я могу применить это в GRASS? А так да, выглядит очень круто, я просто сражён.

[Denis Rykov]

Написал небольшой скрипт, который можно запустить в консоли QGIS и продолжения линий замкнуться (пока без их физического объединения). Одно ограничение: каждый сегмент должен быть пересекаем только с одним другим, то есть ситуация, когда, например, в порог попадают три и более концов разных линий не учитывается, да и непонятно в этом случае как поступать. Попробуйте, пожалуйста, кто-нибудь на своих тестовых данных. Да, и тему, думаю, нужно куда-нибудь переместить. Вот что у меня получилось в итоге:

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

FirstVertex = 0

SecondVertex = 1

LastVertex = -1

PenultVertex = -2



def calculateNodes(features, thresh):

    nodes = []

    for i in features[:-1]:

        igeom = i.geometry().asPolyline()

        for j in features[features.index(i)+1:]:

            jgeom = j.geometry().asPolyline()

            distances = {

                'ff': pow(igeom[FirstVertex].sqrDist(jgeom[FirstVertex]), 0.5),

                'fl': pow(igeom[FirstVertex].sqrDist(jgeom[LastVertex]), 0.5),

                'lf': pow(igeom[LastVertex].sqrDist(jgeom[FirstVertex]), 0.5),

                'll': pow(igeom[LastVertex].sqrDist(jgeom[LastVertex]), 0.5)

            }

            min_key = min(distances, key=distances.get)

            if distances[min_key] <= thresh:

                nodes.append(dict(source=i, target=j, direction=min_key))

    return nodes



def getCross(source, target, direction):

    geom_source = source.geometry().asPolyline()

    geom_target = target.geometry().asPolyline()

    if direction == 'ff':

        p1,p2 = geom_source[FirstVertex], geom_source[SecondVertex]

        p3,p4 = geom_target[FirstVertex], geom_target[SecondVertex]

    elif direction == 'fl':

        p1,p2 = geom_source[FirstVertex], geom_source[SecondVertex]

        p3,p4 = geom_target[LastVertex], geom_target[PenultVertex]

    elif direction == 'lf':

        p1,p2 = geom_source[LastVertex], geom_source[PenultVertex]

        p3,p4 = geom_target[FirstVertex], geom_target[SecondVertex]

    elif direction == 'll':

        p1,p2 = geom_source[LastVertex], geom_source[PenultVertex]

        p3,p4 = geom_target[LastVertex], geom_target[PenultVertex]

    k1 = (p2.y()-p1.y()) / (p2.x()-p1.x())

    k2 = (p4.y()-p3.y()) / (p4.x()-p3.x())

    b1 = p2.y() - k1*p2.x()

    b2 = p4.y() - k2*p4.x()       

    return ((b1-b2)/(k2-k1), (k2*b1 - k1*b2)/(k2-k1)) if k1 != k2 else None



def getVertexId(feature, position):

    if position == 'f':

        vertexId = 0

    elif position == 'l':

        vertexId = len(feature.geometry().asPolyline()) - 1

    return vertexId



def main():

    canvas = qgis.utils.iface.mapCanvas()

    layer = qgis.utils.iface.activeLayer()

    features = [f for f in layer.getFeatures()]

    thresh = 200   

    layer.startEditing()

    eu = QgsVectorLayerEditUtils(layer)

    canvas.refresh()

    for node in calculateNodes(features, thresh):

        cross = getCross(source=node['source'], target=node['target'], direction=node['direction'])

        source_position = node['direction'][0]

        target_position = node['direction'][1]

        if cross is not None:

            eu.moveVertex(cross[0], cross[1], node['source'].id(), getVertexId(node['source'], source_position))

            eu.moveVertex(cross[0], cross[1], node['target'].id(), getVertexId(node['target'], target_position))

    canvas.refresh()



main()

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

Денис, у меня скрипт отрабатывает без ошибок, но линии не объединяет. Просто делает текущий слой редактируемым. Порог менял, не помогло.

[Denis Rykov]

Скрипт не объединяет линии, он просто удлинняет каждую до точки пересечения, объединение - это отдельная задача. И так не делает? А данные в какой СК?

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

Скрипт не объединяет линии, он просто удлинняет каждую до точки пересечения, объединение - это отдельная задача. И так не делает? А данные в какой СК?

Я просто неверно выразился, имел в виду удлинение, конечно. Всё, что он делает, я написал выше, больше ничего. Пока пробовал в местной СК на базе Гаусса-Крюгера / Пулково 42. Проверю ещё и в других СК (хотя если там те же метры, то какая разница?).

[Denis Rykov]

А можешь данные приложить?

[ Сообщение с мобильного устройства ]

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

А можешь данные приложить?

Да там просто несколько линий, специально нарисованных разомкнутыми (т.е. это никакие не данные). А может быть дело в способе запуска скрипта? Я делаю в Питон-консоли " execfile '~/путь/к/скрипту' ".

[Denis Rykov]

Попробовал и таким способом запустить - работает.

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

execfile("/home/rykov/test.py")

Я сейчас в скрипте заметил, что пропущена одна строка (определение переменной canvas), добавил ее, но без нее вообще не работало - выводилось сообщение об ошибке. Попробуй еще раз с новой версией, может в этом дело было.

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

Теперь заработало (правда, в другом проекте с другой СК). Денис, спасибо
Переменную "canvas" я сам добавлял, но, видимо, как-то не так.