Generujesz wielokąty tej samej wielkości zgodnie z PyQGIS?

Chciałbym utworzyć wielokąty wzdłuż linii, aby użyć ich w AtlasCreator w następnym kroku.

ArcMap ma narzędzie o nazwie Funkcje indeksu map paskowych .

Za pomocą tego narzędzia mogę wybrać wysokość i szerokość moich wielokątów (powiedzmy 8 km x 4 km) i automatycznie produkować / obracać je wzdłuż linii.

Jednym z generowanych atrybutów każdego wielokąta jest kąt obrotu, którego potrzebuję, aby później obrócić moje strzałki północy w Atlas Generatorze.

Czy ktoś ma pomysł, jak rozwiązać to zadanie w QGIS / z pyQGIS? Równie dobrze byłyby algorytmy trawiaste lub SAGA albo profesjonalny zestaw narzędzi, który można by użyć w niestandardowej wtyczce;) Edytuj1: Potrzebuję nie tylko zakresu drukowania, ale także samych wielokątów, ponieważ chcę wydrukować mapę za pomocą wszystkie wielokąty / zakresy jako rodzaj mapy przeglądowej.

Edycja2: Oferuję nagrodę, ponieważ wciąż szukam rozwiązania PyQGIS, którego można użyć w wtyczce QGIS bez potrzeby instalowania oprogramowania oprócz QGIS (brak RDBMS jak PostGIS / Oracle)

Interesujące pytanie! Chciałem spróbować, więc spróbowałem.

Możesz to zrobić w PostGRES / POSTGIS za pomocą funkcji, która generuje zestaw wielokątów.

W moim przypadku mam tabelę z jedną funkcją (MULTILINESTRING), która reprezentuje linię kolejową. Musi używać CRS w metrach, używam osgb (27700). Zrobiłem 4 strony x 2 km „stron”.

Tutaj możesz zobaczyć wynik ... zielony jest siecią drogową, przyciętą do 1 km bufora wokół linii kolejowej, co ładnie odpowiada wysokości wielokątów.

Oto funkcja ...

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

Korzystanie z tej funkcji

Oto przykład; Strony 4 km x 2 km, epsg: 27700 i 10% nakładania się

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

Po uruchomieniu możesz następnie wyeksportować z PgAdminIII do pliku csv. Możesz zaimportować to do QGIS, ale może być konieczne ręczne ustawienie CRS dla warstwy - QGIS nie używa SRID w EWKT, aby ustawić dla ciebie CRS warstwy: /

Dodanie atrybutu namiaru

Prawdopodobnie łatwiej to zrobić w postgis, można to zrobić za pomocą wyrażeń QGIS, ale trzeba będzie napisać trochę kodu. Coś takiego...

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

Ostrzeżenia

Jest trochę zhakowany i miał okazję przetestować tylko jeden zestaw danych.

Nie jestem w 100% pewien, które dwa wierzchołki musisz wybrać w tej aktualizacji atrybutu namiaru query. Może być konieczne eksperymentowanie.

Muszę wyznać, że nie mam pojęcia, dlaczego muszę wykonać tak skomplikowaną formułę, aby obrócić wielokąt w celu dopasowania do bieżącego segmentu linii. Myślałem, że mogę użyć danych wyjściowych z ST_Azimuth () w ST_Rotate (), ale najwyraźniej nie.

Istnieją różne rozwiązania. I to może działać z prostą polilinią i wieloma wybranymi elementami

Schemat blokowy:

1. Parametry

   1. wybierz orientację dla generowania i odczytu indeksu (od lewej do prawej, z północy na południe ...)
   2. ustaw rozmiar obiektu

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

   3. zdefiniuj współczynnik superpozycji (domyślnie 10%?)
2. w tym
   1. Porządkowanie polilinii (porównaj punkt początkowy i końcowy) porządkowanie zależy od wybranego ustawienia orientacji> utwórz klasę funkcji zamawiania wierzchołków OrderNodes

3. zapętlić na węzłach OrderNodes

   1. stwórz pierwszy punkt jako kotwicę

   2. dla każdego wierzchołka dodaj go na dykta x, y, id i oblicz wektor

   3. generuj wielokąt (na całej długości i orientacji wektorowej) ze zmniejszającą się superpozycją (10% / 2)> 5% lewy wielokąt 5% prawy wielokąt z tym samym punktem kotwiczenia
   4. Zatrzymaj się, gdy poprzedzający punkt wierzchołka jest poza wielokątem lub jeśli długość wektora ma> kształtować długość
   5. Generuj wielokąt z poprzednim dobrym rozwiązaniem i ustaw punkt zaczepienia na ostatniej dobrej pozycji
   6. Wykonaj nową pętlę i zresetuj dict x, y, id, aby wygenerować następny obiekt wielokąta.

Możesz zmienić tę propozycję, jeśli nie jest tak naprawdę jasna lub komentowana.

Steven Kays odpowiada w pyqgis. Po prostu wybierz linie w warstwie przed uruchomieniem skryptu. Skrypt nie obsługuje łączenia linem, więc nie może działać na warstwie z wielowarstwowym testowaniem

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

Dwie odpowiedzi (w momencie publikacji) są genialne i dobrze wyjaśnione. Jednak istnieje również BARDZO proste, ale skuteczne rozwiązanie (zakładając, że zaakceptujesz wszystkie mapy ustawione w kierunku północnym do góry w tradycyjny sposób, a nie losowy kierunek północny oparty na rzece). Jeśli chcesz rotacji, jest to możliwe, ale nieco bardziej złożone (patrz na dole).

Najpierw spójrz na mój post tutaj . Daje to instrukcje tworzenia map dla Atlasu. Metodą, którą chcesz zastosować, jest adaptacja „Workflow 2” w poradniku. Podziel element liniowy na wierzchołki lub długość i buforuj elementy o dowolną wartość. Ilość buforowana przez ciebie częściowo narzuci nakładanie się (ale patrz poniżej), ale co ważniejsze, tworzy funkcję z obszarem. Możesz użyć dowolnej liczby wtyczek do podziału linii, ale GRASS v.split.length i v.split.vert to dobre opcje (dostępne w Processing Toolbox).

Po włączeniu generowania atlasu w programie Map Composer i wybraniu buforowanej warstwy, wróć do karty elementów i wybierz obiekt mapy. Zaznacz opcję „Kontrolowane przez Atlas”, a w twoim przypadku wybrałbym opcję Margin wokół. To kontroluje nakładanie się map (alternatywnie możesz preferować stałą skalę).

Możesz wyświetlić podgląd Atlasa za pomocą przycisku Podgląd atlasu na górnym pasku narzędzi kompozytora i zobaczyć, ile stron on wyprodukuje. Uwaga: możesz wyeksportować wszystkie strony w jednym pliku PDF lub jako osobne pliki.

Aby mapa obracała się wzdłuż linii, we właściwościach elementu do tworzenia mapy znajduje się pole rotacji. Musisz ustawić wyrażenie (użyj małego przycisku po prawej stronie pola rotacji). Wybierz zmienną jako opcję, a następnie Edytuj. Pojawi się okno konstruktora wyrażeń, w którym można uzyskać dostęp do geometrii lub pól elementów atlasu. Następnie możesz zbudować ekspres, aby obrócić mapę zgodnie z obrotem obiektów (możesz obliczyć namiar przy użyciu punktów początkowych i końcowych każdego segmentu linii oraz odrobiny wyzwalacza). Powtórz ten sam proces, aby obrócić strzałkę Północ (używając tego samego wyrażenia lub wstępnie obliczonej zmiennej).