使用PyQGIS沿线生成大小相等的多边形？

#1 楼

有趣的问题！这是我想尝试的事情，所以就去尝试一下。就我而言，我有一张表，其中包含一个代表铁路线的要素（MULTILINESTRING）。它需要以米为单位使用CRS，我正在使用osgb（27700）。我已经完成了4公里x 2公里的“页面”。多边形的高度很好。



这里是函数...

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

使用此函数

这是一个例子； 4公里x 2公里的页面，epsg：27700，重叠10％
您可以将其导入QGIS，但是您可能需要手动为图层设置CRS-QGIS不会在EWKT中使用SRID为您设置图层CRS：/

添加方位角属性

这可能在postgis中更容易完成，可以在QGIS表达式中完成，但是您需要编写一些代码。像这样的东西...

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

注意事项

有点黑，并且只有一次可以对一个数据集进行测试。

不是100％不确定需要在方位属性更新query上选择哪两个顶点。..可能需要试验。

我必须承认我不知道为什么我需要做一个复杂的公式来旋转多边形以匹配当前线段。我以为我可以在ST_Rotate（）中使用ST_Azimuth（）的输出，但似乎不能。

#2 楼

Steven Kays用pyqgis回答。
在运行脚本之前，只需在图层中选择行。

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

#3 楼

有不同的解决方案。这可以用于简单的折线和多个选定的实体

框图：



参数


>选择生成方向并读取索引（从左到右，从北到南...）
设置对象大小

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

/>

定义叠加系数（默认情况下为10％？）


init



订购折线（比较开始和端点）的排序取决于您的方向选择>创建顶点排序要素类OrderNodes



在OrderNodes上循环


创建您第一个点作为锚点
每个顶点将其添加到dict x，y，id上并计算向量
通过减少重叠（10％/ 2）> 5％生成多边形（在长度和方向上）左多边形5％右多边形，具有相同的锚点
当先验顶点不在多边形范围内或矢量len>达到形状长度时停止
用先前的好解生成多边形nd设置锚点的最后一个好位置
执行新循环并重置字典x，y，id以生成下一个多边形对象。



您可以更改此命题如果不清楚或没有评论。

#4 楼

这两个答案（在发布时）非常巧妙，并得到了很好的解释。但是，还有一种非常简单但有效的解决方案（假设您将接受所有以传统方式与北向上对齐的地图，而不是基于河流的随机北向）。如果您要轮换，则可以，但是稍微复杂一点（请参阅底部）。

首先在这里看看我的帖子。这为您提供了创建Atlas地图覆盖的方法。您想要的方法是方法指南中“工作流程2”的改编。通过顶点或长度分割线性特征，并以任意数量缓冲特征。缓冲的量将部分决定重叠（但请参见下文），但更重要的是，它创建具有面积的要素。您可以使用任意数量的插件来分割行，但是GRASS v.split.length和v.split.vert是不错的选择（在Processing Toolbox中可用）。

在Map Composer中启用了Atlas Generation并选择了缓冲层，然后切换回项目标签并选择地图对象。选中“由Atlas控制”，在您的用例中，我会选择围绕功能的边距。这将控制您在地图之间的重叠（或者，您可能更喜欢固定比例的地图）。请注意，您可以选择将所有页面导出为单个PDF或单独的文件。

为了使地图沿直线旋转，在Map Composer项属性中有一个旋转字段。您将需要设置一个表达式（使用旋转框右侧的小按钮）。选择变量作为选项，然后选择编辑。将弹出一个表达式构建器，您可以在其中访问图集要素的几何或字段。然后，您可以构建一个Express，以根据要素的旋转来旋转地图（可以使用每个线段的起点和终点以及一点点Trig来计算方位角）。重复相同的过程以旋转北箭头（使用相同的表达式或预先计算的变量）。

#5 楼

我知道这是一个老问题，但是已经创建了一个新插件来解决此问题。 https://plugins.qgis.org/plugins/polystrip/

我不是创建者，也不为插件效劳。