[GIS] Generating equal sized polygons along line with PyQGIS

grasspyqgisqgis-processingsagastrip-map

I would like to create polygons along a line to use them for AtlasCreator in a next step.

ArcMap has a tool called Strip Map Index Features.

With this tool I can choose the height and width of my polygons (say 8km x 4km) and produce/rotate them along the line automatically.

One of the generated attributes of each polygon is the rotation angle that I need to rotate my north arrows in Atlas Generator afterwards.

enter image description here

Does anyone have an idea how to solve this task in QGIS / with PyQGIS?

GRASS- or SAGA-algorithms or a prossessing-toolbox-model which could be used inside a custom plugin would be fine, too.

I need not only the print extents but also the polygons itself as I want to print a map with all polygons/extents as some sort of overview map. I am looking for a PyQGIS-solution that can be used in a QGIS-Plugin without the need for installing a software aside from QGIS (no RDBMS like PostGIS / Oracle).

Best Answer

Interesting question! It's something I've wanted to try myself, so gave it a go.

You can do this in PostGRES/POSTGIS with a function which generates a set of polygons.

In my case, I have a table with one feature (a MULTILINESTRING) which represents a railway line. It needs to use a CRS in meters, I'm using osgb (27700). I've done 4km x 2km 'pages'.

Here, you can see the result...the green stuff is the road network, clipped to a 1km buffer around the railway, which corresponds to the height of the polygons nicely.

postgis generated strip map

Here's the function...

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

Using this function

Here's an example; 4km x 2km pages, epsg:27700 and 10% overlap

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

After running this you can then export from PgAdminIII into a csv file. You can import this into QGIS, but you might well need to set the CRS manually for the layer - QGIS doesn't use the SRID in EWKT to set the layer CRS for you :/

Adding bearing attribute

This is probably easier done in postgis, it can be done in QGIS expressions but you'll need to write some code. Something like this...

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

It's a bit hacked-together, and only had a chance to test on one dataset.

Not 100% sure which two vertices you'll need to choose on that bearing attribute update query.. might need to experiment.

I must confess I have no idea why I need to do such a convoluted formula to rotate the polygon to match the current line segment. I thought I could use the output from ST_Azimuth() in ST_Rotate(), but seemingly not.

Related Question