Создание полигонов одинакового размера по линии PyQGIS?

42

Я хотел бы создать полигоны вдоль линии, чтобы использовать их для AtlasCreator на следующем этапе.

ArcMap имеет инструмент под названием Strip Map Index Features .

С помощью этого инструмента я могу выбрать высоту и ширину моих полигонов (скажем, 8 км х 4 км) и автоматически производить / вращать их вдоль линии.

Один из сгенерированных атрибутов каждого многоугольника - это угол поворота, который мне нужен, чтобы потом повернуть стрелки на север в Генераторе Атласа.

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

У кого-нибудь есть идеи, как решить эту задачу в QGIS / с pyQGIS? Алгоритмы Grass- или SAGA-алгоритмы или модель набора инструментов-инструментов, которые можно использовать внутри собственного плагина, тоже подойдут;) Edit1: мне нужны не только экстенты печати, но и сами полигоны, так как я хочу напечатать карту с все полигоны / экстенты как своего рода обзорная карта.

Edit2: я предлагаю щедрость, поскольку я все еще ищу решение PyQGIS, которое можно использовать в плагине QGIS без необходимости установки программного обеспечения помимо QGIS (без RDBMS, такой как PostGIS / Oracle)

Berlinmapper
источник
4
Это похоже на забавную идею для плагина.
alphabetasoup
1
Как дикая мысль, я думаю, что что-то, основанное на обобщении
Пекера
1
возможно v.split.length, затем нарисуйте прямую линию между начальной и конечной точкой сегментов, а затем v.buffer с опцией «Не делать заглушки на концах полилиний»
Томас Б.
1
Я хотел бы получить награду за этот вопрос, но у меня пока недостаточно репутации; (
Berlinmapper 12.12.15
2
В реализациях "label-follow line" может быть многократно используемый код. Ваши прямоугольники похожи на следы глифов какого-то моноширинного шрифта.
user30184 15.12.15

Ответы:

29

Интересный вопрос! Это то, что я хотел попробовать сам, так что попробуй.

Вы можете сделать это в PostGRES / POSTGIS с помощью функции, которая генерирует набор полигонов.

В моем случае у меня есть таблица с одной функцией (MULTILINESTRING), которая представляет железнодорожную линию. Нужно использовать CRS в метрах, я использую osgb (27700). Я сделал "страницы" 4 км х 2 км.

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

сгенерированная Postgis карта

Вот функция ...

CREATE OR REPLACE FUNCTION getAllPages(wid float, hite float, srid integer, overlap float) RETURNS SETOF geometry AS
$BODY$
DECLARE
    page geometry; -- holds each page as it is generated
    myline geometry; -- holds the line geometry
    startpoint geometry;
    endpoint geometry;
    azimuth float; -- angle of rotation
    curs float := 0.0 ; -- how far along line left edge is
    step float;
    stepnudge float;
    currpoly geometry; -- used to make pages
    currline geometry;
    currangle float;
    numpages float;
BEGIN
    -- drop ST_LineMerge call if using LineString 
    -- replace this with your table.
    SELECT ST_LineMerge(geom) INTO myline from traced_osgb; 
    numpages := ST_Length(myline)/wid;

    step := 1.0/numpages;
    stepnudge := (1.0-overlap) * step; 
    FOR r in 1..cast (numpages as integer)
    LOOP
        -- work out current line segment

        startpoint :=  ST_SetSRID(ST_Line_Interpolate_Point(myline,curs),srid);
        endpoint :=  ST_SetSRID(ST_Line_Interpolate_Point(myline,curs+step),srid);
        currline := ST_SetSRID(ST_MakeLine(startpoint,endpoint),srid);

        -- make a polygon of appropriate size at origin of CRS
        currpoly := ST_SetSRID(ST_Extent(ST_MakeLine(ST_MakePoint(0.0,0.0),ST_MakePoint(wid,hite))),srid);

        -- then nudge downwards so the midline matches the current line segment
        currpoly := ST_Translate(currpoly,0.0,-hite/2.0);

        -- Rotate to match angle
        -- I have absolutely no idea how this bit works. 
        currangle := -ST_Azimuth(startpoint,endpoint) - (PI()/2.0) + PI();
        currpoly := ST_Rotate(currpoly, currangle);

        -- then move to start of current segment
        currpoly := ST_Translate(currpoly,ST_X(startpoint),ST_Y(startpoint));

        page := currpoly;

        RETURN NEXT page as geom; -- yield next result
        curs := curs + stepnudge;
    END LOOP;
    RETURN;
END
$BODY$
LANGUAGE 'plpgsql' ;

Используя эту функцию

Вот пример; Страницы 4 x 2 км, epsg: 27700 и перекрытие 10%

select st_asEwkt(getallpages) from getAllPages(4000.0, 2000.0, 27700, 0.1);

После запуска вы можете экспортировать из PgAdminIII в CSV-файл. Вы можете импортировать это в QGIS, но вам может потребоваться установить CRS вручную для слоя - QGIS не использует SRID в EWKT, чтобы установить CRS слоя для вас: /

Добавление атрибута подшипника

Это, вероятно, легче сделать в postgis, это можно сделать в выражениях QGIS, но вам нужно написать некоторый код. Что-то вроде этого...

create table pages as (
    select getallpages from getAllPages(4000.0, 2000.0, 27700, 0.1)
);

alter table pages add column bearing float;

update pages set bearing=ST_Azimuth(ST_PointN(getallpages,1),ST_PointN(getallpages,2));

Предостережения

Он немного взломан, и у него была возможность протестировать только один набор данных.

Не уверен на 100%, какие две вершины вам нужно выбрать при обновлении атрибута подшипника query. Возможно, придется поэкспериментировать.

Должен признаться, я понятия не имею, почему мне нужно сделать такую ​​сложную формулу, чтобы повернуть многоугольник, чтобы соответствовать текущему отрезку. Я думал, что мог бы использовать вывод из ST_Azimuth () в ST_Rotate (), но, похоже, нет.

Стивен Кей
источник
Ваш ответ действительно хорош, и я обязательно попробую кое-что. Одно из ограничений для меня заключается в том, что я не могу использовать postgres для проекта, над которым я работаю, и мне нужно что-то, что не зависит от серверной стороны. Но, возможно, я могу использовать ваш отличная логика для воспроизведения чего-то подобного с pyQGIS.
Берлинмаппер
2
если это так, взгляните на класс QgsGeometry . Он имеет подмножество геометрических операций PostGIS и будет хорошей отправной точкой, если вы хотите пойти по пути pyQGIS. Алгоритм должен быть переносимым на pyQGIS ..
Стивен Кей
3
Я думаю, что для postgis подход, использующий ST_Simplify для генерации опорных линий и разбиения линии на сегменты, а затем использующий ST_Buffer и ST_Envelope , будет короче и эффективнее.
Матиас Кун
@Matthias Kuhn: если я разделю линию на сегменты, я могу получить линии одинакового размера, но необязательно также многоугольники одинакового размера. например, если линия довольно «соблазнительная», полигон, вероятно, будет короче, не так ли?
Berlinmapper
2
Я проверил ваше решение и PyQGIS-версию вашего скрипта. Осталась какая-то идея, как решить некоторые незначительные проблемы: bit.ly/1KL7JHn ?
Berlinmapper
12

Есть разные решения. И это может работать с простой ломаной линией и несколькими выбранными объектами

блок-схема:

  1. параметры

    1. выберите ориентацию для генерации и считайте индекс (слева направо, с севера на юг ...)
    2. установить размер объекта

    shape = (4000,8000) # (<width>,<length>)
    1. определить суперпозицию coef (10% по умолчанию?)
  2. в этом
    1. Упорядочение полилинии (сравнение начальной и конечной точек) упорядочение зависит от вашего выбора ориентации> создать упорядочение вершин, класс объектов OrderNodes
  3. цикл на узлах заказов

    1. создать первую точку в качестве якоря

    2. для каждой вершины добавьте ее в dict x, y, id и вычислите вектор

    3. создать полигон (по длине и векторной ориентации) с уменьшением суперпозиции (10% / 2)> 5% левого многоугольника 5% правого многоугольника с той же точкой привязки
    4. Остановиться, когда предыдущая точка вершины выходит за пределы многоугольника или если вектор len> для формирования длины
    5. Создать полигон с предыдущим хорошим решением и установить точку привязки с последней хорошей позицией
    6. Выполните новый цикл и сбросьте dict x, y, id для создания следующего объекта многоугольника.

Вы можете изменить это предложение, если оно не совсем понятно или прокомментировано.

GeoStoneMarten
источник
это звучит сложно, но я должен признать, что я еще не знаю, как использовать это для моделиста или PyQGIS. кстати: что такое коэффициент суперпозиции?
Berlinmapper
@Berlinmapper, который является частью многоугольника с суперпозицией 8000 x 10% в этом случае. Вы можете выбрать другой или создать фиксированное расстояние между полигонами. Вы можете видеть это во всех атласах, чтобы указать следующую страницу плиток в углу
GeoStoneMarten
Ваше решение предназначено для использования с pyQGIS или набором инструментов для обработки? это звучит здорово, но я все еще не знаю, как поступить
Berlinmapper
1
@Berlinmapper Я думаю, вам нужно использовать pyQGIS для создания сценария процесса и установки входных и выходных параметров в наборе инструментов обработки или плагине QGIS. То же, что и arcgistoolbox. У меня нет времени на самом деле, чтобы сделать это и проверить это.
GeoStoneMarten
12

Стивен Кейс отвечает в Pyqgis. Просто выберите линии в вашем слое перед запуском скрипта. Скрипт не поддерживает линейное объединение, поэтому он не может работать на слое с множеством строк

#!python
# coding: utf-8

# https://gis.stackexchange.com/questions/173127/generating-equal-sized-polygons-along-line-with-pyqgis
from qgis.core import QgsMapLayerRegistry, QgsGeometry, QgsField, QgsFeature, QgsPoint
from PyQt4.QtCore import QVariant


def getAllPages(layer, width, height, srid, overlap):
    for feature in layer.selectedFeatures():
        geom = feature.geometry()
        if geom.type() <> QGis.Line:
            print "Geometry type should be a LineString"
            return 2
        pages = QgsVectorLayer("Polygon?crs=epsg:"+str(srid), 
                      layer.name()+'_id_'+str(feature.id())+'_pages', 
                      "memory")
        fid = QgsField("fid", QVariant.Int, "int")
        angle = QgsField("angle", QVariant.Double, "double")
        attributes = [fid, angle]
        pages.startEditing()
        pagesProvider = pages.dataProvider()
        pagesProvider.addAttributes(attributes)
        curs = 0
        numpages = geom.length()/(width)
        step = 1.0/numpages
        stepnudge = (1.0-overlap) * step
        pageFeatures = []
        r = 1
        currangle = 0
        while curs <= 1:
            # print 'r =' + str(r)
            # print 'curs = ' + str(curs)
            startpoint =  geom.interpolate(curs*geom.length())
            endpoint = geom.interpolate((curs+step)*geom.length())
            x_start = startpoint.asPoint().x()
            y_start = startpoint.asPoint().y()
            x_end = endpoint.asPoint().x()
            y_end = endpoint.asPoint().y()
            # print 'x_start :' + str(x_start)
            # print 'y_start :' + str(y_start)
            currline = QgsGeometry().fromWkt('LINESTRING({} {}, {} {})'.format(x_start, y_start, x_end, y_end))
            currpoly = QgsGeometry().fromWkt(
                'POLYGON((0 0, 0 {height},{width} {height}, {width} 0, 0 0))'.format(height=height, width=width))
            currpoly.translate(0,-height/2)
            azimuth = startpoint.asPoint().azimuth(endpoint.asPoint())
            currangle = (startpoint.asPoint().azimuth(endpoint.asPoint())+270)%360
            # print 'azimuth :' + str(azimuth)
            # print 'currangle : ' +  str(currangle)

            currpoly.rotate(currangle, QgsPoint(0,0))
            currpoly.translate(x_start, y_start)
            currpoly.asPolygon()
            page = currpoly
            curs = curs + stepnudge
            feat = QgsFeature()
            feat.setAttributes([r, currangle])
            feat.setGeometry(page)
            pageFeatures.append(feat)
            r = r + 1

        pagesProvider.addFeatures(pageFeatures)
        pages.commitChanges()
        QgsMapLayerRegistry.instance().addMapLayer(pages)
    return 0

layer = iface.activeLayer()
getAllPages(layer, 500, 200, 2154, 0.4)
lejedi76
источник
1
Отлично. Я проверил решение. Любая идея, как решить эти проблемы, решение все еще имеет: bit.ly/1KL7JHn ?
Berlinmapper
возможно, здесь есть какое-то «вдохновение»: github.com/maphew/arcmapbook/blob/master/Visual_Basic/…
Томас Б.
Спасибо. отличный ресурс, чтобы понять, как работает инструмент ArcMap. К сожалению, я не привык к VB, но, возможно, кто-то другой может использовать его для публикации ответа / комментария;)
Berlinmapper
4

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

Сначала взгляните на мой пост здесь . Это дает вам инструкции по созданию покрытий карты для Atlas. Метод, который вы хотите адаптировать - это «Рабочий процесс 2» в руководстве. Разделите ваш линейный объект по вершинам или длине и буферизируйте объекты по любому количеству Величина, по которой вы буферизуете, частично будет определять перекрытие (но см. Ниже), но, что более важно, оно создает элемент с областью. Вы можете использовать любое количество плагинов для разделения строк, но GRASS v.split.length и v.split.vert являются хорошими вариантами (доступны в Processing Toolbox).

Включив Генерацию атласа в Map Composer и выбрав буферный слой, вернитесь на вкладку элементов и выберите объект карты. Отметьте «Контролируемый Atlas», и в вашем случае я бы выбрал функцию Margin вокруг. Это будет контролировать ваше перекрытие между картами (в качестве альтернативы вы можете предпочесть фиксированный масштаб).

Вы можете предварительно просмотреть свой Atlas с помощью кнопки Preview Atlas на верхней панели инструментов композитора и посмотреть, сколько страниц он будет создавать. Обратите внимание, что вы можете выбрать экспорт всех страниц в один PDF или в виде отдельных файлов.

Чтобы заставить карту вращаться вдоль линии, в свойствах элемента «Компоновщик карты» имеется поле поворота. Вам нужно будет установить выражение (используйте маленькую кнопку справа от поля поворота). Выберите переменную в качестве вашего варианта и затем Изменить. Появится конструктор выражений, и там вы сможете получить доступ к геометрии или полям объектов атласа. Затем вы можете построить экспресс для поворота карты в соответствии с поворотом объектов (вы можете рассчитать направление, используя начальную и конечную точки каждого отрезка и немного триггера). Повторите тот же процесс, чтобы повернуть стрелку севера (используя то же выражение или предварительно рассчитанную переменную).

MappaGnosis
источник
спасибо за это решение. но я думаю, что таким образом я бы не получил полигоны с единичными экстентами, которые я хочу напечатать. Мне нужно, чтобы они также создали «обзорную карту» со всеми экстентами печати.
Berlinmapper