Gleichgroße Polygone nach PyQGIS generieren?

Ich möchte Polygone entlang einer Linie erstellen, um sie in einem nächsten Schritt für AtlasCreator zu verwenden.

ArcMap verfügt über ein Tool mit dem Namen " Strip Map Index Features" .

Mit diesem Werkzeug kann ich die Höhe und Breite meiner Polygone (z. B. 8 km x 4 km) auswählen und sie automatisch entlang der Linie erzeugen / drehen.

Eines der generierten Attribute jedes Polygons ist der Drehwinkel, den ich benötige, um meine Nordpfeile anschließend in Atlas Generator zu drehen.

Hat jemand eine Idee, wie man diese Aufgabe in QGIS / mit pyQGIS löst? Gras- oder SAGA-Algorithmen oder ein Prossessing-Toolbox-Modell, das in einem benutzerdefinierten Plugin verwendet werden könnte, wären auch in Ordnung Alle Polygone / Bereiche als eine Art Übersichtskarte.

Edit2: Ich biete ein Kopfgeld an, da ich noch nach einer PyQGIS-Lösung suche, die in einem QGIS-Plugin verwendet werden kann, ohne dass eine Software außerhalb von QGIS installiert werden muss (kein RDBMS wie PostGIS / Oracle)

Interessante Frage! Es ist etwas, das ich selbst ausprobieren wollte, also habe ich es ausprobiert.

Sie können dies in PostGRES / POSTGIS mit einer Funktion tun, die eine Menge von Polygonen erzeugt.

In meinem Fall habe ich eine Tabelle mit einem Merkmal (MULTILINESTRING), das eine Eisenbahnlinie darstellt. Es muss ein CRS in Metern verwendet werden, ich verwende osgb (27700). Ich habe 4 km x 2 km Seiten erstellt.

Hier können Sie das Ergebnis sehen ... das Grüne ist das Straßennetz, das an einem 1 km langen Puffer um die Eisenbahn befestigt ist, der der Höhe der Polygone gut entspricht.

Hier ist die Funktion ...

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

Verwendung dieser Funktion

Hier ist ein Beispiel; 4 km x 2 km große Seiten, epsg: 27700 und 10% Überlappung

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

Nachdem Sie dies ausgeführt haben, können Sie von PgAdminIII in eine CSV-Datei exportieren. Sie können dies in QGIS importieren, müssen das CRS jedoch möglicherweise manuell für den Layer festlegen. QGIS verwendet die SRID in EWKT nicht, um das Layer-CRS für Sie festzulegen: /

Lagerattribut hinzufügen

Dies ist wahrscheinlich einfacher in Postgis, es kann in QGIS-Ausdrücken gemacht werden, aber Sie müssen etwas Code schreiben. Etwas wie das...

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

Vorbehalte

Es ist ein bisschen zusammen gehackt und hatte nur die Chance, an einem Datensatz zu testen.

Nicht 100% sicher, welche zwei Eckpunkte Sie für die Aktualisierung der queryLagerattribute auswählen müssen. Möglicherweise müssen Sie experimentieren.

Ich muss zugeben, dass ich keine Ahnung habe, warum ich eine derart verschlungene Formel ausführen muss, um das Polygon so zu drehen, dass es mit dem aktuellen Liniensegment übereinstimmt. Ich dachte, ich könnte die Ausgabe von ST_Azimuth () in ST_Rotate () verwenden, aber anscheinend nicht.

Es gibt verschiedene Lösungen. Dies funktioniert mit einfachen Polylinien und mehreren ausgewählten Objekten

Blockdiagramm:

1. Parameter

   1. Wählen Sie die Ausrichtung für die Erstellung und lesen Sie den Index (von links nach rechts, von Nord nach Süd ...).
   2. Objektgröße einstellen

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

   3. Überlagerungskoeffizient definieren (standardmäßig 10%?)
2. drin
   1. Die Reihenfolge der Polylinien (Start- und Endpunkt vergleichen) hängt von Ihrer Wahl der Ausrichtung ab

3. Schleife auf OrderNodes

   1. erstelle deinen ersten punkt als anker

   2. Fügen Sie für jeden Scheitelpunkt das Diktat x, y, id hinzu und berechnen Sie einen Vektor

   3. Erzeugen eines Polygons (über die Länge und Vektororientierung) mit reduzierter Überlagerung (10% / 2)> 5% linkes Polygon 5% rechtes Polygon mit demselben Ankerpunkt
   4. Halten Sie an, wenn ein vorhergehender Scheitelpunkt außerhalb des Polygons liegt oder wenn der Vektor len> die Länge formen soll
   5. Erzeugen Sie ein Polygon mit der vorherigen guten Lösung und setzen Sie den Ankerpunkt mit der letzten guten Position
   6. Führen Sie eine neue Schleife durch und setzen Sie dict x, y, id zurück, um das nächste Polygonobjekt zu generieren.

Sie können diesen Satz ändern, wenn er nicht wirklich klar ist oder nicht kommentiert wird.

Steven Kays antwortet in Pyqgis. Wählen Sie einfach die Linien in Ihrer Ebene aus, bevor Sie das Skript ausführen. Das Skript unterstützt das Zusammenführen von Linien nicht und kann daher nicht auf Ebenen mit Mehrfachlinienfolgen arbeiten

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

Die beiden Antworten (zum Zeitpunkt der Veröffentlichung) sind genial und gut erklärt. Es gibt jedoch auch eine SEHR einfache, aber effektive Lösung dafür (vorausgesetzt, Sie akzeptieren alle Karten, die nach Norden ausgerichtet sind, auf herkömmliche Weise und nicht in zufälliger Richtung nach Norden, basierend auf dem Fluss). Wenn Sie Rotationen wollen, ist es möglich, aber etwas komplexer (siehe unten).

Schauen Sie sich zuerst meinen Beitrag hier an . Auf diese Weise erhalten Sie eine Anleitung zum Erstellen von Kartenabdeckungen für Atlas. Die gewünschte Methode ist eine Anpassung von 'Workflow 2' in der Anleitung. Teilen Sie Ihr lineares Feature nach Scheitelpunkten oder Länge und puffern Sie die Features um einen beliebigen Betrag. Der Betrag, um den Sie puffern, bestimmt teilweise die Überlappung (siehe unten), aber was noch wichtiger ist, es wird ein Feature mit einer Fläche erstellt. Sie können eine beliebige Anzahl von Plugins verwenden, um die Zeilen zu teilen. Die Optionen GRASS v.split.length und v.split.vert sind jedoch gut geeignet (verfügbar in der Processing Toolbox).

Nachdem Sie die Atlas-Generierung in Map Composer aktiviert und Ihren gepufferten Layer ausgewählt haben, wechseln Sie zurück zur Registerkarte Elemente und wählen Sie Ihr Kartenobjekt aus. Aktivieren Sie das Kontrollkästchen "Von Atlas gesteuert", und in Ihrem Anwendungsfall würde ich mich für die Funktion "Margin around" entscheiden. Dadurch wird die Überlappung zwischen den Karten gesteuert (alternativ können Sie einen festen Maßstab bevorzugen).

Sie können Ihren Atlas in der Vorschau anzeigen, indem Sie auf die Schaltfläche Atlas in der oberen Symbolleiste des Komponisten klicken und sehen, wie viele Seiten erzeugt werden. Beachten Sie, dass Sie auswählen können, ob alle Seiten in einer einzelnen PDF-Datei oder als separate Dateien exportiert werden sollen.

In den Eigenschaften des Map Composer-Elements befindet sich ein Drehfeld, mit dem die Karte entlang der Linie gedreht werden kann. Sie müssen einen Ausdruck festlegen (verwenden Sie die kleine Schaltfläche rechts neben dem Drehfeld). Wählen Sie Variable als Ihre Option und dann Bearbeiten. Ein Ausdrucksgenerator öffnet sich und Sie können dort auf die Geometrie oder Felder der Atlas-Features zugreifen. Sie können dann einen Express erstellen, um die Karte entsprechend der Drehung der Features zu drehen (Sie können die Peilung anhand der Start- und Endpunkte jedes Liniensegments und einer Triggerung berechnen). Wiederholen Sie denselben Vorgang, um den Nordpfeil zu drehen (mit demselben Ausdruck oder derselben vorberechneten Variablen).