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

40

Me gustaría crear polígonos a lo largo de una línea para usarlos para AtlasCreator en un próximo paso.

ArcMap tiene una herramienta llamada Características del índice de mapa de franjas .

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 hacia el norte en Atlas Generator después.

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

Edit2: ofrezco una recompensa porque sigo buscando una solución PyQGIS que se puede utilizar en un complemento de QGIS sin la necesidad de instalar un software aparte de QGIS (no RDBMS como PostGIS / Oracle)

    
pregunta Berlinmapper 08.12.2015 - 22:36

4 respuestas

28

Pregunta interesante! Es algo que he querido probar yo mismo, así que decidí intentarlo.

Puedes 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 (un MULTILINESTRING) que representa una línea ferroviaria. Necesita usar un CRS en metros, estoy usando osgb (27700). He hecho 'páginas' de 4 km x 2 km.

Aquí, puedes ver el resultado ... lo verde es la red de carreteras, recortada a un búfer de 1 km alrededor del ferrocarril, que se corresponde con la altura de los polígonos.

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

Utilizando esta función

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

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 en QGIS, pero es posible que deba 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 rumbo

Esto es probablemente más fácil de hacer en Postgis, se puede hacer en expresiones QGIS pero necesitarás 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));

Caveats

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

No estoy 100% seguro de qué dos vértices deberás elegir en esa actualización de atributo de rumbo query .. podría necesitar experimentar.

Debo confesar que no tengo ni 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 podía usar la salida de ST_Azimuth () en ST_Rotate (), pero aparentemente no.

    
respondido por el Steven Kay 15.12.2015 - 17:57
13

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 la generación y el índice de lectura (de izquierda a derecha, de norte a sur ...)
    2. establecer tamaño de objeto

    shape = (4000,8000) # (<width>,<length>)
    
    1. defina la superposición de coef (10% por defecto?)
  2. init
    1. El orden de la polilínea (comparar los puntos de inicio y final) depende de su elección de orientación > crear vértices ordenando featureclass 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 la orientación del vector) con una superposición reductora (10% / 2) > Polígono izquierdo del 5% Polígono derecho del 5% 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 es > para dar forma a la longitud
    5. Generar polígono con una buena solución anterior y establecer un punto de anclaje con la última posición buena
    6. Realizar un nuevo bucle y restablecer dict x, y, id para generar el siguiente objeto poligonal.

Puede cambiar esta propuesta si no está realmente clara o comentar.

    
respondido por el GeoStoneMarten 13.12.2015 - 15:00
11

Steven Kays responde en pyqgis. Simplemente seleccione las líneas en su capa antes de ejecutar el script. La secuencia de comandos no es compatible con el linemerging, 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)
    
respondido por el lejedi76 25.01.2016 - 02:39
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 (asumiendo que aceptará todos los mapas alineados con el Norte de la manera tradicional, en lugar de una dirección Norte aleatoria basada en el río). Si quieres rotaciones, es posible pero un poco más complejo (ver la parte inferior).

Primero eche un vistazo a mi publicación aquí . Esto le brinda instrucciones para crear coberturas de mapas para Atlas. El método que desea es una adaptación de es 'Workflow 2' en el manual. Divide tu entidad lineal por vértices o longitud y almacena en búfer las entidades por cualquier cantidad. La cantidad por la que se almacenará en el búfer determinará parcialmente la superposición (pero consulte a continuación), pero lo que es más importante, crea una característica 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 la Caja de herramientas de procesamiento).

Habiendo habilitado la Generación Atlas en Map Composer y seleccionando tu capa amortiguada, regresa a la pestaña de elementos y selecciona tu 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 los mapas (alternativamente, puede preferir una escala fija).

Puedes previsualizar tu Atlas usando 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 hacer que el mapa gire a lo largo de la línea, hay un campo de rotación en las propiedades del elemento Map Composer. Necesitará establecer una expresión (use el botón pequeño a la derecha del cuadro de rotación). Seleccione la variable como su opción y luego Editar. Aparecerá un generador de expresiones y allí podrá acceder a la geometría o campos de las características del atlas. Luego puede crear un expreso para rotar el mapa según la rotación de las entidades (puede calcular el rumbo 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).

    
respondido por el MappaGnosis 15.12.2015 - 20:07

Lea otras preguntas en las etiquetas