Générer des polygones de taille égale le long de PyQGIS?

J'aimerais créer des polygones le long d'une ligne pour les utiliser pour AtlasCreator dans une prochaine étape.

ArcMap dispose d'un outil appelé Caractéristiques de l'index de carte .

Avec cet outil, je peux choisir la hauteur et la largeur de mes polygones (disons 8 km x 4 km) et les produire / les faire pivoter automatiquement le long de la ligne.

L'un des attributs générés de chaque polygone est l'angle de rotation dont j'ai besoin pour faire pivoter ensuite mes flèches vers le nord dans Atlas Generator.

Quelqu'un a-t-il une idée de la façon de résoudre cette tâche dans QGIS / avec pyQGIS? Des algorithmes Grass ou SAGA ou un modèle prossessing-toolbox qui pourrait être utilisé dans un plugin personnalisé conviendraient également;) Edit1: J'ai besoin non seulement de l'étendue de l'impression, mais également des polygones, car je souhaite imprimer une carte avec tous les polygones / étendues comme une sorte de carte générale.

Edit2: J'offre une prime car je suis toujours à la recherche d'une solution PyQGIS pouvant être utilisée dans un plugin QGIS sans qu'il soit nécessaire d'installer un logiciel en dehors de QGIS (aucun SGBDR comme PostGIS / Oracle)

Question interessante! C'est quelque chose que j'ai voulu essayer moi-même, alors essayez-le.

Vous pouvez le faire dans PostGRES / POSTGIS avec une fonction qui génère un ensemble de polygones.

Dans mon cas, j'ai une table avec une caractéristique (un MULTILINESTRING) qui représente une ligne de chemin de fer. Il faut utiliser un CRS en mètres, j'utilise osgb (27700). J'ai fait 4 km 'pages'.

Ici, vous pouvez voir le résultat ... la substance verte est le réseau routier, attaché à un tampon de 1 km autour de la voie ferrée, ce qui correspond joliment à la hauteur des polygones.

Voici la fonction ...

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

Utiliser cette fonction

Voici un exemple. 4 km x 2 km pages, epsg: 27700 et 10% de chevauchement

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

Après avoir exécuté ceci, vous pouvez ensuite exporter de PgAdminIII dans un fichier csv. Vous pouvez importer cela dans QGIS, mais vous aurez peut-être besoin de définir manuellement le CRS pour le calque - QGIS n'utilise pas le SRID dans EWKT pour définir le CRS du calque pour vous: /

Ajout d'un attribut de relèvement

Ceci est probablement plus facile à faire dans les postgis, cela peut être fait dans les expressions QGIS mais vous devrez écrire du code. Quelque chose comme ça...

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

Mises en garde

C'est un peu bidouillé, et nous n'avons eu qu'une chance de tester sur un jeu de données.

Vous n'êtes pas sûr à 100% des deux sommets que vous devrez choisir pour cette mise à jour d'attribut de roulement query.

Je dois avouer que je ne sais pas pourquoi je dois utiliser une formule aussi compliquée pour faire pivoter le polygone afin de le faire correspondre au segment de ligne actuel. Je pensais pouvoir utiliser la sortie de ST_Azimuth () dans ST_Rotate (), mais apparemment pas.

Il y a différentes solutions. Et cela peut fonctionner avec une polyligne simple et plusieurs entités sélectionnées

diagramme:

1. Paramètres

   1. sélectionner l'orientation pour la génération et lire l'index (de gauche à droite, du nord au sud ...)
   2. définir la taille de l'objet

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

   3. définir le coef de superposition (10% par défaut?)
2. init
   1. L'ordre des polylignes (comparer les points de départ et d'arrivée) dépend de votre choix d'orientation> créer un ordre de classement de sommets avec une classe

3. boucle sur OrderNodes

   1. crée ton premier point comme ancre

   2. pour chaque sommet, ajoutez-le à dict x, y, id et calculez un vecteur

   3. générer un polygone (sur la longueur et l'orientation du vecteur) avec superposition réductrice (10% / 2)> 5% polygone gauche 5% polygone droit avec le même point d'ancrage
   4. Arrêtez-vous lorsqu'un sommet précédent est en dehors du polygone ou si le vecteur len est> pour former la longueur
   5. Générer un polygone avec la bonne solution précédente et définir le point d'ancrage avec la dernière bonne position
   6. Effectuez une nouvelle boucle et réinitialisez dict x, y, id pour générer le prochain objet polygone.

Vous pouvez changer cette proposition si elle n'est pas vraiment claire ou commentée.

Steven Kays répond en pyqgis. Il suffit de sélectionner les lignes de votre couche avant d'exécuter le script. Le script ne prend pas en charge la sortie de ligne, il ne peut donc pas fonctionner sur une couche à plusieurs chaînes.

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

Les deux réponses (au moment de la publication) sont ingénieuses et bien expliquées. Cependant, il existe aussi une solution TRÈS simple mais efficace (en supposant que vous acceptiez toutes vos cartes alignées avec le nord de la manière traditionnelle, plutôt qu’une direction aléatoire du Nord basée sur la rivière). Si vous voulez des rotations, c'est possible mais un peu plus complexe (voir en bas).

Regardez d'abord mon post ici . Cela vous explique comment créer des couvertures de carte pour Atlas. La méthode que vous souhaitez utiliser est une adaptation du "flux de travail 2" du didacticiel. Divisez votre entité linéaire par des sommets ou une longueur et mettez-la en tampon de n'importe quelle quantité. La quantité que vous mettez en mémoire tampon dicte partiellement le chevauchement (mais voir ci-dessous), mais surtout, elle crée une entité avec une zone. Vous pouvez utiliser n'importe quel nombre de plugins pour diviser les lignes, mais GRASS v.split.length et v.split.vert sont de bonnes options (disponibles dans Processing Toolbox).

Après avoir activé Atlas Generation dans Map Composer et sélectionné votre couche en mémoire tampon, revenez à l'onglet Eléments et sélectionnez votre objet cartographique. Cochez la case 'Contrôlé par Atlas', et dans votre cas d'utilisation, j'opterais pour la fonctionnalité Marge autour. Cela contrôlera votre chevauchement entre les cartes (sinon, vous préférerez peut-être une échelle fixe).

Vous pouvez prévisualiser votre Atlas à l’aide du bouton Aperçu Atlas situé dans la barre d’outils supérieure du compositeur et voir le nombre de pages qu’il produira. Notez que vous pouvez choisir d'exporter toutes les pages dans un seul PDF ou sous forme de fichiers séparés.

Pour faire pivoter la carte le long de la ligne, il existe un champ de rotation dans les propriétés de l'élément Compositeur de carte. Vous devrez définir une expression (utilisez le petit bouton à droite de la case de rotation). Sélectionnez variable comme option, puis Modifier. Un constructeur d'expression apparaîtra dans lequel vous pourrez accéder à la géométrie ou aux champs des entités de l'atlas. Vous pouvez ensuite créer un express pour faire pivoter la carte en fonction de la rotation des entités (vous pouvez calculer le relèvement en utilisant les points de début et de fin de chaque segment de ligne et un peu de trig). Répétez le même processus pour faire pivoter votre flèche du Nord (en utilisant la même expression ou une variable précalculée).