¿Generando polígonos de igual tamaño a lo largo de la línea con PyQGIS?

Me gustaría crear polígonos a lo largo de una línea para usarlos para AtlasCreator en el siguiente paso.

ArcMap tiene una herramienta llamada Strip Index Index Features .

Con esta herramienta puedo elegir la altura y el ancho de mis polígonos (por ejemplo, 8 km x 4 km) y producirlos / rotarlos a lo largo de la línea automáticamente.

Uno de los atributos generados de cada polígono es el ángulo de rotación que necesito para rotar mis flechas norte en Atlas Generator después.

¿Alguien tiene una idea de cómo resolver esta tarea en QGIS / con pyQGIS? Los algoritmos Grass o SAGA o un modelo de caja de herramientas de prossessing que podría usarse dentro de un complemento personalizado también estaría bien;) Edit1: necesito no solo las extensiones de impresión sino también los polígonos en sí, ya que quiero imprimir un mapa con todos los polígonos / extensiones como algún tipo de mapa general.

Edit2: estoy ofreciendo una recompensa ya que todavía estoy buscando una solución PyQGIS que pueda usarse en un complemento QGIS sin la necesidad de instalar un software aparte de QGIS (sin RDBMS como PostGIS / Oracle)

¡Interesante pregunta! Es algo que quería probar yo mismo, así que probé.

Puede hacer esto en PostGRES / POSTGIS con una función que genera un conjunto de polígonos.

En mi caso, tengo una tabla con una característica (una MULTILINESTRING) que representa una línea de ferrocarril. Necesita usar un CRS en metros, estoy usando osgb (27700). He hecho 4km x 2km 'páginas'.

Aquí puede ver el resultado ... lo verde es la red de carreteras, sujeta a un búfer de 1 km alrededor del ferrocarril, que corresponde a la altura de los polígonos muy bien.

Aquí está la función ...

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 función

Aquí hay un ejemplo; Páginas de 4 km x 2 km, epsg: 27700 y superposición del 10%

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

Después de ejecutar esto, puede exportar desde PgAdminIII a un archivo csv. Puede importar esto a QGIS, pero es posible que necesite configurar el CRS manualmente para la capa: QGIS no usa el SRID en EWKT para configurar el CRS de capa para usted: /

Agregar atributo de demora

Probablemente esto sea más fácil de hacer en postgis, se puede hacer en expresiones QGIS pero necesitará escribir algo de código. Algo como esto...

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

Advertencias

Está un poco pirateado, y solo tuvo la oportunidad de probar en un conjunto de datos.

No estoy 100% seguro de qué dos vértices necesitará elegir en la actualización de ese atributo de demora query... podría necesitar experimentar.

Debo confesar que no tengo idea de por qué necesito hacer una fórmula tan complicada para rotar el polígono para que coincida con el segmento de línea actual. Pensé que podría usar la salida de ST_Azimuth () en ST_Rotate (), pero aparentemente no.

Hay diferentes soluciones. Y esto puede funcionar con polilíneas simples y múltiples entidades seleccionadas

diagrama de bloques:

1. Parámetros

   1. seleccione la orientación para el índice de generación y lectura (de izquierda a derecha, de norte a sur ...)
   2. establecer el tamaño del objeto

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

   3. definir coeficiente de superposición (10% por defecto?)
2. en eso
   1. Los pedidos de polilínea (comparar el punto inicial y final) dependen de su elección de orientación> crear una clase de entidad de orden de vértices OrderNodes

3. bucle en OrderNodes

   1. crea tu primer punto como ancla

   2. para cada vértice agréguelo en dict x, y, id y calcule un vector

   3. generar polígono (sobre la longitud y orientación del vector) con superposición reductora (10% / 2)> 5% polígono izquierdo 5% polígono derecho con el mismo punto de anclaje
   4. Deténgase cuando un punto de vértice precedente esté fuera del polígono o si el vector len tiene> la longitud de la forma
   5. Genere un polígono con una buena solución anterior y establezca el punto de anclaje con la última buena posición
   6. Realice un nuevo bucle y restablezca dict x, y, id para generar el siguiente objeto de polígono.

Puede cambiar esta propuesta si no es realmente clara o no comenta.

Respuesta de Steven Kays en pyqgis. Simplemente seleccione las líneas en su capa antes de ejecutar el script. La secuencia de comandos no admite la combinación de líneas, por lo que no puede funcionar en la capa con multilínea

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

Las dos respuestas (en el momento de la publicación) son ingeniosas y están bien explicadas. Sin embargo, también existe una solución MUY simple pero efectiva para esto (suponiendo que acepte todos sus mapas alineados con el Norte hacia arriba de la manera tradicional, en lugar de una dirección aleatoria del Norte basada en el río). Si quieres rotaciones, es posible pero un poco más complejo (ver abajo).

Primero eche un vistazo a mi publicación aquí . Esto le proporciona una guía práctica para crear coberturas de mapas para Atlas. El método que desea es una adaptación de 'Workflow 2' en el tutorial. Divide tu entidad lineal por vértices o longitud y amortigua las entidades por cualquier cantidad. La cantidad de búfer dictará parcialmente la superposición (pero ver más abajo), pero lo más importante, crea una función con un área. Puede usar cualquier cantidad de complementos para dividir las líneas, pero GRASS v.split.length y v.split.vert son buenas opciones (disponibles en Processing Toolbox).

Una vez que haya habilitado Atlas Generation en Map Composer y haya seleccionado su capa almacenada, vuelva a la pestaña de elementos y seleccione su objeto de mapa. Marque 'Controlado por Atlas', y en su caso de uso, optaría por la función Margen alrededor. Esto controlará su superposición entre mapas (alternativamente, puede preferir la escala fija).

Puede obtener una vista previa de su Atlas utilizando el botón Vista previa del Atlas en la barra de herramientas superior del compositor y ver cuántas páginas producirá. Tenga en cuenta que puede elegir exportar todas las páginas en un solo PDF o como archivos separados.

Para que el mapa gire a lo largo de la línea, hay un campo de rotación en las propiedades del elemento del Compositor de mapas. Deberá establecer una expresión (use el pequeño botón a la derecha del cuadro de rotación). Seleccione variable como su opción y luego Editar. Aparecerá un generador de expresiones y allí podrá acceder a la geometría o los campos de las características del atlas. Luego puede construir un expreso para rotar el mapa de acuerdo con la rotación de las entidades (puede calcular la demora utilizando los puntos de inicio y final de cada segmento de línea y un poco de trigonometría). Repita el mismo proceso para rotar su flecha Norte (usando la misma expresión o variable calculada previamente).