[GIS] Creating line pattern fill using geometry generator (PyQGIS)

geometry-generatorpyqgisqgis-3

I'm trying to make a line pattern fill in QGIS using the geometry generator and a Python function.

The code is inspired from this answer : Save line geometry as points geometry in QGIS (pyqgis)

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def fill_line_pattern(x, feature, parent):

#create a base line that is 45 degrees
base_xmin = 2285865.0
base_xmax = 2669093.0
base_ymin = 1017594.0
base_ymax = 1399851.0

#calculating interval between parrallels
countX = int((2800000.0-2285865.0) / xInterval)
#empty list
pt=[]
pts = []
intersected_lines = []
#raw lines creation
for xOff in range(countX+1):
    ptXmin = base_xmin + xOff*(xInterval)
    ptXmax = base_xmax + xOff*(xInterval)
    ptmin = QgsPoint(ptXmin,base_ymin)
    ptmax = QgsPoint(ptXmax,base_ymax)
    pt = ([ptmin, ptmax])
    pts.append(pt)
#intersect lines with polygons
#for f in feat :
for p in pts :
    rawlines = QgsGeometry.fromPolyline(p)
    if feature.geometry().intersects(rawlines):
        geom = feature.geometry().intersection(rawlines)
        intersected_lines.append(geom.asPolyline())
#create lines
for line_ok in intersected_lines:
    return QgsGeometry.fromPolylineXY(line_ok)

The issue is that the first line only is crossing the polygon…

How can I have all the lines crossing the polygons?

In comparison to the inspirational code above, there's no loop for the polygon features.

In the function, the feature called is a QgsFeature which can't use getQgsFeatures().

If a for loop is used :

for f in feature :
    for p in pts :
    rawlines = QgsGeometry.fromPolyline(p)
        if f.geometry().intersects(rawlines):
            geom = f.geometry().intersection(rawlines)
            intersected_lines.append(geom.asPolyline())

This error message appears:

Evaluation error: 'int' object has no attribute 'geometry'

EDIT

The code below does almost what I want, it seems a little bit buggy because not all the intersections appear… The other disadvantage is the long time QGIS takes to render it on the canvas and in the composer.

On the other hand, the pdf output size is so light and is competitive with other software.

Maps are made to be shared: QGIS biggest issue is PDF output file size

The inspiration of the source : https://www.lutraconsulting.co.uk/blog/2015/06/05/qgis-function-editor/

Lutra Consulting uses a cache but I didn't find the way to do it…

from qgis.core import *
from qgis.gui import *

def _layer_intersection(layer, xInterval):
#make 45 degrees lines from layer extent
    ext = layer.extent()
    xmax = ext.xMaximum()
    ymax = ext.yMaximum()
    xmin = ext.xMinimum()
    ymin = ext.yMinimum()

    pts_hash_box = []
    line_in_polygon = []

    x_base_min = xmin - (ymax-ymin)
    countX = int(((xmax-xmin)+(ymax-ymin)) / xInterval)

    for xOff in range(countX+1):
        ptXmin = x_base_min + xOff*(xInterval)
        ptXmax = xmin + xOff*(xInterval)

        ptmin_hash_box = QgsPointXY(ptXmin,ymin)
        ptmax_hash_box = QgsPointXY(ptXmax,ymax)

        pt_hash_box = ([ptmin_hash_box, ptmax_hash_box])
        pts_hash_box.append(pt_hash_box)

    multipolylines = QgsGeometry.fromMultiPolylineXY(pts_hash_box)

    #intersection between 45 degrees virtual lines and the polygon features
    for layer_feature in layer.getFeatures() :
        if layer_feature.geometry().intersects(multipolylines):
            pts = layer_feature.geometry().intersection(multipolylines).asMultiPolyline()
            for pt in pts :
                line_in_polygon.append(y)

return QgsGeometry.fromMultiPolylineXY(line_in_polygon)

@qgsfunction(args='auto', group='Custom')
def line_pattern_fill(layer_name, xInterval, feature, parent):
for layer in QgsProject.instance().mapLayers().values():
    if layer.name() == layer_name :
        return _layer_intersection(layer, xInterval)

Best Answer

EDIT 2

This bug is finally resolved (wait the version 3.24 (2022)) : https://github.com/qgis/QGIS/pull/45654

Here is my working solution. I took a great part of your code but rearranged and improved it.

Create a custom function in the function editor and paste the code below.

Use it as : line_pattern_fill(@layer_id, 45, 1000) where 45 is the angle in degrees and 1000 is the X interval in layer units.

I create, as Lutra Consulting blog post, a cache. I reinitialize the variable if the angle or x_interval change so that QGIS compute new pattern fills.

from qgis.core import *
from qgis.gui import *
from math import ceil, cos, radians, sin

cache = {}


@qgsfunction(args='auto', group='geom_generator')
def line_pattern_fill(layer_id, angle, x_interval, feature, parent):
    """
    Returns a MultiLineString geometry with line pattern fill, contained
    in the current feature geometry.
    <h4>Syntax</h4>
    <div class="syntax">
    <code><span class="functionname">line_pattern_fill</span>
    (<span class="argument">@layer_id</span>,
    <span class="argument">angle</span>,
    <span class="argument">x interval</span>)</code>
    <h4>Arguments</h4>
    <div class="arguments">
    <table><tr>
    <td class="argument">@layer_id</td><td>ID of current layer</td></tr><tr>
    <td class="argument">angle</td><td>angle in degrees, from 0 to 179</td></tr><tr>
    <td class="argument">x interval</td><td>spacing between lines, in layer units</td></tr>
    </table>
    </div>
    <h4>Examples</h4>
    <div class="examples">
    <ul>
    <li><code>line_pattern_fill(@layer_id, 45, 1000)</code> &rarr; 
    <code>a line pattern fill at 45 degrees, spaced by 1000 layer units</code></li>
    <li><code>line_pattern_fill(@layer_id, 0, 60)</code> &rarr; 
    <code>a vertical line pattern fill, spaced by 60 layer units</code></li>
    <li><code>line_pattern_fill(@layer_id, 90, 150)</code> &rarr; 
    <code>a horizontal line pattern fill, spaced by 150 layer units</code></li>
    </ul>
    </div>
    """
    global cache
    key1 = f"{angle}_{x_interval}"
    key2 = feature.id()
    if layer_id in cache:
        if not key1 in cache[layer_id]:
            cache[layer_id] = {key1: {}}
    else:
        cache[layer_id] = {key1: {}}

    if key2 not in cache[layer_id][key1]:
        pts_hash_box = []
        bbox = feature.geometry().boundingBox()
        xcoords = [bbox.xMinimum(), bbox.xMaximum()]
        ycoords = [bbox.yMinimum(), bbox.yMaximum()]

        xcoef = round(cos(radians(angle)), 5)
        ycoef = round(sin(radians(angle)), 5)
        fn = lambda x, y: ycoef * y - xcoef * x

        ccorners = [fn(x, y) for x in xcoords for y in ycoords]
        cmin = min(ccorners)
        cmax = max(ccorners)
        cvar = cmax - cmin

        for i in range(0, ceil(cvar), x_interval):
            c = cmin + i
            if xcoef == 0:
                xs = []
            else:
                xs = [((ycoef * y - c) / xcoef, y) for y in ycoords]

            if ycoef == 0:
                ys = []
            else:
                ys = [(x, (xcoef * x + c) / ycoef) for x in xcoords]

            pts = [
                (x, y)
                for x, y in xs + ys
                if xcoords[0] <= x <= xcoords[1] and ycoords[0] <= y <= ycoords[1]
            ]
            pt_hash_box = [QgsPointXY(x, y) for x, y in sorted(pts)]
            pts_hash_box.append(pt_hash_box)

        multipolylines = QgsGeometry.fromMultiPolylineXY(pts_hash_box)
        geom = multipolylines.intersection(feature.geometry())
        if geom.type() == QgsWkbTypes.LineGeometry:
            cache[layer_id][key1][key2] = geom
        elif geom.wkbType() == QgsWkbTypes.GeometryCollection:
            if geom.convertGeometryCollectionToSubclass(QgsWkbTypes.LineGeometry):
                cache[layer_id][key1][key2] = geom
            else:
                cache[layer_id][key1][key2] = QgsGeometry()
        else:
            cache[layer_id][key1][key2] = QgsGeometry()

    return cache[layer_id][key1][key2]

EDIT

Code improved : QgsGeometry.intersection() produce sometimes a GeometryCollection with lines and points. If the result of intersection is a collection, the code try to convert it to Line geometry type.
Angle choice added : the angle argument is now in the code. The angle is in degrees. Rewritten code for this.
Change wkbType to type : LineGeometry is more generic and works with MultiLineStringZ and all other Line types.
Formatting help : adapt function help to QGIS expression function help format

Best Answer

Related Solutions

pyqgis – Creating Centroid Point Symbol Layer via Geometry Generator in PyQGIS

Related Question