Gerando polígonos de tamanho igual ao longo da linha com PyQGIS?

Gostaria de criar polígonos ao longo de uma linha para usá-los no AtlasCreator na próxima etapa.

O ArcMap possui uma ferramenta chamada Strip Map Index Features .

Com esta ferramenta, posso escolher a altura e a largura dos meus polígonos (digamos 8 km x 4 km) e produzi-los / rotacioná-los ao longo da linha automaticamente.

Um dos atributos gerados de cada polígono é o ângulo de rotação necessário para girar minhas setas norte no Atlas Generator posteriormente.

Alguém tem uma idéia de como resolver esta tarefa no QGIS / com pyQGIS? Algum algoritmo Grass ou SAGA ou um modelo de caixa de ferramentas de prosseguimento que poderia ser usado dentro de um plug-in personalizado também seria bom;) Edit1: Eu preciso não apenas das extensões de impressão, mas também dos polígonos, pois quero imprimir um mapa com todos os polígonos / extensões como algum tipo de mapa de visão geral.

Edit2: Estou oferecendo uma recompensa, pois ainda estou procurando uma solução PyQGIS que possa ser usada em um plug-in QGIS sem a necessidade de instalar um software além do QGIS (nenhum RDBMS como PostGIS / Oracle)

Pergunta interessante! É algo que eu queria experimentar, então tentei.

Você pode fazer isso no PostGRES / POSTGIS com uma função que gera um conjunto de polígonos.

No meu caso, tenho uma tabela com um recurso (uma MULTILINESTRING) que representa uma linha ferroviária. Ele precisa usar um CRS em metros, estou usando osgb (27700). Eu fiz 4 km x 2 km 'páginas'.

Aqui, você pode ver o resultado ... o material verde é a rede de estradas, presa a um buffer de 1 km ao redor da ferrovia, o que corresponde muito bem à altura dos polígonos.

Aqui está a função ...

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' ;

Usando esta função

Aqui está um exemplo; Páginas de 4 km x 2 km, epsg: 27700 e 10% de sobreposição

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

Depois de executar isso, você poderá exportar do PgAdminIII para um arquivo csv. Você pode importar isso para o QGIS, mas pode ser necessário definir o CRS manualmente para a camada - o QGIS não usa o SRID no EWKT para definir o CRS da camada para você: /

Adicionando atributo de rolamento

Provavelmente, isso é mais fácil no postgis, pode ser feito nas expressões QGIS, mas você precisará escrever algum código. Algo assim...

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));

Ressalvas

É um pouco hackeado e só teve a chance de testar em um conjunto de dados.

Não tem 100% de certeza de quais dois vértices você precisará escolher nessa atualização de atributo do rolamento query. Talvez seja necessário experimentar.

Devo confessar que não tenho ideia de por que preciso fazer uma fórmula tão complicada para girar o polígono para corresponder ao segmento de linha atual. Eu pensei que poderia usar a saída de ST_Azimuth () em ST_Rotate (), mas aparentemente não.

Existem diferentes soluções. E isso pode funcionar com polilinha simples e várias entidades selecionadas

diagrama de bloco:

1. Parâmetros

   1. selecione a orientação para geração e índice de leitura (da esquerda para a direita, de norte para sul ...)
   2. definir tamanho do objeto

   shape = (4000,8000) # (<width>,<length>)

   3. definir coeficiente de superposição (10% por padrão?)
2. nisso
   1. A ordem da polilinha (comparar o ponto inicial e o ponto final) depende da sua opção de orientação> crie um vértice e faça a ordem da classe

3. loop em OrderNodes

   1. crie seu primeiro ponto como âncora

   2. para cada vértice, adicione-o no dict x, y, id e calcule um vetor

   3. gerar polígono (acima do comprimento e orientação do vetor) com superposição redutora (10% / 2)> 5% polígono esquerdo 5% polígono direito com o mesmo ponto de ancoragem
   4. Pare quando um ponto de vértice precedente estiver fora do polígono ou se o vetor len for> moldar o comprimento
   5. Gere polígono com boa solução anterior e defina o ponto de ancoragem com a última boa posição
   6. Execute um novo loop e redefina o dict x, y, id para gerar o próximo objeto de polígono.

Você pode alterar essa proposição se não estiver realmente claro ou comentar.

Steven Kays responde em pyqgis. Basta selecionar as linhas na sua camada antes de executar o script. O script não suporta a entrada de linha, portanto, não pode funcionar na camada com cadeia de linhas múltiplas

#!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)

As duas respostas (no momento da publicação) são engenhosas e bem explicadas. No entanto, também existe uma solução MUITO simples, mas eficaz, possível (assumindo que você aceitará todos os seus mapas alinhados com o Norte da maneira tradicional, em vez de uma direção aleatória do Norte com base no rio). Se você deseja rotações, é possível, mas um pouco mais complexo (veja abaixo).

Primeiro, dê uma olhada no meu post aqui . Isso fornece instruções sobre como criar coberturas de mapa para o Atlas. O método que você deseja é uma adaptação é 'Workflow 2' no tutorial. Divida seu recurso linear por vértices ou comprimento e armazene em buffer os recursos por qualquer valor. A quantidade que você armazena em buffer ditará parcialmente a sobreposição (mas veja abaixo), mas o mais importante é que ele cria um recurso com uma área. Você pode usar qualquer número de plug-ins para dividir as linhas, mas o GRASS v.split.length e v.split.vert são boas opções (disponíveis em Processing Toolbox).

Depois de ativar a Geração Atlas no Map Composer e selecionar sua camada em buffer, volte para a guia itens e selecione seu objeto de mapa. Marque 'Controlado pelo Atlas' e, no seu caso de uso, eu optaria pelo recurso Margin around. Isso controlará sua sobreposição entre mapas (como alternativa, você pode preferir escala fixa).

Você pode visualizar o Atlas usando o botão Visualizar Atlas na barra de ferramentas superior do compositor e ver quantas páginas ele produzirá. Observe que você pode optar por exportar todas as páginas em um único PDF ou como arquivos separados.

Para fazer o mapa girar ao longo da linha, há um campo de rotação nas propriedades do Item do compositor de mapa. Você precisará definir uma expressão (use o pequeno botão à direita da caixa de rotação). Selecione a variável como sua opção e, em seguida, edite. Um construtor de expressões será exibido e você poderá acessar a geometria ou os campos dos recursos do atlas. Você pode criar um expresso para girar o mapa de acordo com a rotação dos recursos (você pode calcular o rumo usando os pontos inicial e final de cada segmento de linha e um pouco de trigonometria). Repita o mesmo processo para girar sua seta norte (usando a mesma expressão ou variável pré-calculada).