Gleichgroße Polygone nach PyQGIS generieren?

42

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.

Bildbeschreibung hier eingeben

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)

Berlinmapper
quelle
4
Das sieht nach einer lustigen Idee für ein Plugin aus.
Alphabetasoup
1
Als wilder Gedanke denke ich, dass etwas, das auf der Peucker-Douglas-Verallgemeinerung basiert, funktionieren könnte
plablo09
1
vielleicht v.split.length, dann ziehen Sie eine gerade Linie zwischen Start- und Endpunkt der Segmente und dann v.buffer mit der Option "Keine Kappen an den Enden von Polylinien machen"
Thomas B
1
Ich würde gerne ein Kopfgeld für diese Frage
auslösen,
2
In Implementierungen mit "Label-Follow-Line" könnte es wiederverwendbaren Code geben. Ihre Rechtecke sind wie Fußabdrücke von Glyphen einer monospaced Schriftart.
user30184

Antworten:

29

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.

postgis generierte streifen karte

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.

Steven Kay
quelle
Ihre Antwort ist wirklich großartig und ich werde es auf jeden Fall versuchen. Eine Einschränkung für mich ist, dass ich für das Projekt, an dem ich arbeite, keine Postgres verwenden kann und etwas brauche, das nicht von einer Serverseite abhängt. Aber vielleicht kann ich Ihre verwenden Tolle Logik, um so etwas mit pyQGIS zu reproduzieren.
Berlinmapper
2
Wenn dies der Fall ist, werfen Sie einen Blick auf die QgsGeometry-Klasse . Es enthält eine Teilmenge der Geometrieoperationen von PostGIS und ist ein guter Ausgangspunkt, wenn Sie die pyQGIS-Route einschlagen möchten. Der Algorithmus sollte auf pyQGIS portierbar sein.
Steven Kay
3
Ich denke, für Postgis wäre ein Ansatz, der ST_Simplify verwendet , um Referenzlinien zu generieren und die Linie in Segmente zu unterteilen, und dann ST_Buffer und ST_Envelope verwendet , kürzer und effizienter.
Matthias Kuhn
@Matthias Kuhn: Wenn ich die Linie in Segmente zerlege, könnte ich gleich große Linien bekommen, aber nicht unbedingt auch gleich große Polygone. Wenn die Linie beispielsweise ziemlich kurvig ist, ist das Polygon wahrscheinlich kürzer, nicht wahr?
Berlinmapper
2
Ich habe Ihre Lösung und die PyQGIS-Version Ihres Skripts getestet .
Berlinmapper
12

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>)
    1. Ü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.

GeoStoneMarten
quelle
das klingt raffiniert, aber ich muss zugeben, dass ich noch nicht weiß, wie ich das für den Modeler oder PyQGIS verwenden soll. Übrigens: Was ist ein Überlagerungskoeffizient?
Berlinmapper
@Berlinmapper, der in diesem Fall Teil eines Polygons mit der Superposition 8000 x 10% ist. Sie können eine andere auswählen oder eine feste Distanz zwischen Polygonen erstellen. Sie können das in jedem Atlas sehen, um die nächste
Kachelseite
Soll Ihre Lösung mit pyQGIS oder der Verarbeitungs-Toolbox verwendet werden? es hört sich toll an, aber ich weiß immer noch nicht, wie ich vorgehen soll
Berlinmapper
1
@Berlinmapper Ich denke, Sie müssen pyQGIS verwenden, um das Prozessskript zu erstellen und die Eingabe- und Ausgabeparameter in der Verarbeitungs-Toolbox oder im QGIS-Plugin festzulegen. Wie arcgistoolbox. Ich habe eigentlich keine Zeit, das zu tun und es zu testen.
GeoStoneMarten
12

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)
lejedi76
quelle
1
Toll. Ich habe die Lösung getestet. Irgendeine Idee, wie man diese Probleme löst, die Lösung hat noch: bit.ly/1KL7JHn ?
Berlinmapper
Vielleicht gibt es hier etwas "Inspiration": github.com/maphew/arcmapbook/blob/master/Visual_Basic/…
Thomas B
thanks.great resource, um zu verstehen, wie das ArcMap-Tool funktioniert. Leider bin ich nicht an VB gewöhnt, aber vielleicht kann jemand anderes damit eine Antwort / einen Kommentar
posten
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).

MappaGnosis
quelle
danke für diese lösung. aber ich denke, auf diese Weise würde ich keine Polygone der einzelnen Bereiche erhalten, die ich drucken möchte. Ich brauche sie, um auch eine "Übersichtskarte" mit allen Druckausmaßen zu erstellen.
Berlinmapper