3

I wanted in QGIS to be able to number (whole number: 1,2,3,4...) multiple polygons in a sequence according to the the direction of a line intersecting through them (so that I have control over the numbering sequence, unlike @row_number), as shown in the picture bellow:

enter image description here

From the start point of the line each polygon intersected by it will receive a number in an already existing field (number), starting from one (1) until the end of the line as well as the polygons.

Is it possible to do in Field calculator or Virtual layer (preferably in field calculator)?

4

Let's assume there area two layers, a polygon layer 'layer' and a polyline layer 'lines_test', see the image below

input

In the polygon layer create a new field "NewId" using Field calculator with the following expression:

with_variable(
    'lines',
    'lines_test', -- name of the polyline layer
    array_find(
        array_agg($id,
            group_by:=overlay_intersects(
                            layer:=@lines,
                            expression:="id"
                        )[0],
            filter:=line_locate_point(
                            overlay_intersects(
                                layer:=@lines,
                                expression:=$geometry
                            )[0],
                            point_on_surface($geometry)
                        ) IS NOT NULL,
            order_by:=line_locate_point(
                            overlay_intersects(
                                layer:=@lines,
                                expression:=$geometry
                            )[0],
                            point_on_surface($geometry)
                        )
            ),
            $id
    ) + 1
)
    

and get the output

result

Note: 0 is the value for no overlaps

4
  • 1
    Thank you very much. It works perfectly!
    – Felippe M.
    Nov 23 '21 at 19:42
  • 1
    @Taras - nice solution. However, in some cases it does not produce the desired result - see screenshot, where I applied your solution: i.stack.imgur.com/5WYxx.png - I added another solution, similar to yours, that avoids these kind of problems. Alos using point_on_surface takes some time to calculate. I took the centroid of the interesction between line and polygon: this is faster to calculate.
    – Babel
    Dec 18 '21 at 19:46
  • 1
    Thank you so much for your comment!
    – Taras
    Dec 18 '21 at 20:04
  • It takes especially long with layers with many polygons (the houses of a town, a neighborhood...) : even if you have just one single line that crosses a few houses, it still calculates the array_agg($id,group_by:=overlay_intersects()... part for every polygon feature again and again. In such cases, it takes a very long time to finally calculate just a few numbers. I added a way to overcome this and use a possibility to do the count for several lines at once. Still, your solution is great, OP didn't ask for a layer with larger number of features - it was just my curiosity how to improve that.
    – Babel
    Dec 18 '21 at 23:24
3

Here is an example of how this can be done. Replace the values of the first four variables to match your use case:

    # -*- coding: utf8 -*-
    from qgis.PyQt.QtCore import *
    from qgis.PyQt.QtGui import *
    from qgis.PyQt.QtWidgets import *
    from qgis.utils import iface
    from qgis.core import *
    import processing

    layers = [layer for layer in \
             QgsProject.instance().mapLayers().values()]
    # change name to layername you wish to number
    layer_to_number = "parcels2021"
    # change to match name of linelayer
    linelayer = "lines"
    # change to match name of field
    numfield = "lotnumber"
    # change to match unique field in polygon layer
    uniquefield = "capakey"

    layer = None
    line_layer = None

    for lyr in list(QgsProject.instance().mapLayers().values()):
        if lyr.name() == layer_to_number:
            layer = lyr
        if lyr.name() == linelayer:
            line_layer = lyr
    
    geom = None        
    
    if line_layer:
        linefeatlist = list(line_layer.getFeatures())
        geom = linefeatlist[0].geometry()


    result = processing.run("native:pointsalonglines",{'INPUT':           line_layer,'DISTANCE':2,'START_OFFSET':1,'END_OFFSET':0,'OUTPUT':'memory:'})

    QgsProject.instance().addMapLayer(result['OUTPUT'])

    outputlayer =  QgsProject.instance().mapLayersByName("output")[0]

    layer_provider=layer.dataProvider()
    layer_provider.addAttributes([QgsField(numfield,QVariant.Int)])
    layer.updateFields()

    #put parcel features in array and fill it with parcels that intersect with the line geometry
    parcels = []
    featlist = list(layer.getFeatures())
    for feat in featlist:
        if feat.geometry().intersects(geom):
            parcels.append(feat)

    field_idx = layer.fields().indexOf(numfield)
    
    pointlist = list(outputlayer.getFeatures())
    counter = 1
    #make list to store unique field for parcellayer
    capakeys = []
    for point in pointlist:
        pgeom = point.geometry()
        for parcel in parcels:
            capa = parcel[uniquefield]
            if capa in capakeys:
                continue
            if pgeom.intersects(parcel.geometry()):
                capakeys.append(capa)
                fid = parcel.id()
                layer.startEditing()
                layer.changeAttributeValue(fid, field_idx, counter)               
                layer.commitChanges()
                counter += 1
            
    print("ready") 
1

The principles

This answer is similar to the one posted by @Taras, but avoids some of the problems it creates (see below). It improves efficiency for larger layers and allows for the use of different lines for different countings as can be seen on this screenshot: enter image description here

How it works

The basic idea works as follows: on the polygon layer, get for each polygon the line that crosses the polygon and create the centroid of the intersecting part (centroid of that part of the crossing line that is inside the polygon). Then measure how far away this centroid is from the line's start point. Order all these distances in an array and then for each point (again) measure the distance from the start point and look up in the array at what position (index) it can be found. As array indices start with 0 for the first element, add 1.

The solution

For optimizing with several lines see below.

Use this expression:

array_find ( 
   array_sort( 
        array_filter(  
            array_agg(  
                line_locate_point( 
                    overlay_crosses(
                        'line', 
                        $geometry
                    )[0],
                    array_foreach (
                        overlay_crosses(
                            'line', 
                            $geometry
                        ),
                        centroid (
                            intersection ($geometry, @element)
                        )
                    )[0]
                )
            ), 
            @element > 0
        )
    ),
    line_locate_point( 
        overlay_crosses(
            'line', 
            $geometry
        )[0],
        array_foreach (
            overlay_crosses('line', $geometry),
            centroid (intersection ($geometry, @element)
        )
    )[0]
    )
)+1

Here, the expression is used as to label the polygon layer. The expression ignores to polygon that is not crossed by the line (red arrow) and correctly numbers polygons that are crossed by the line more then once: polygons 8 (3 times), 10, 11 as well as 1 (two times each): enter image description here

Optimizing the expression for several lines

If you have different lines you want to use for successiveley number the polygons they cross, use an improved version of the expression.

  1. First, however, on the polygon layer create a field lineid that gets the unique id value of the line that crosses each polygon. Use this expression to create the field:

    overlay_crosses('line', $id)[0]

  2. Now, you can use this expression for the numbering:

     array_find (
         array_sort(
             array_filter(  
                 array_agg(
                     line_locate_point(
                         geometry ( get_feature_by_id ('line', lineid)),
                         centroid (
                             intersection ($geometry, geometry (get_feature_by_id ('line', lineid)))
                         )
                     ),
                     group_by:=lineid
                 ),
                 @element > 0
             )
         ),
         line_locate_point( 
             geometry (get_feature_by_id ('line', lineid)),
             centroid (
                 intersection ($geometry, geometry (get_feature_by_id ('line', lineid)))
             )
         )
     )+1
    

The expression in use with real world data from OSM. The numbers are created as dynamic label with the expression from above. Based on this value, the color is set to a color ramp (Viridis) in increasing values from the first to the last number:

enter image description here

The problem

The solution by @Taras works, but in some cases does not create the correct output, see the same data as before with the expression by @Taras, even using a straight line without crossing the same polygon more than once. Also, using point_on_surface quite takes a while to calculate, so the solution above avoids that and is faster.

enter image description here

0

I think the best way to do this is to add vertices to your line (for example every 2 meter) and then with python you can loop over the vertices and check for intersection with the parcellayer. If there is intersection, then you write the value of the counter (starting with 1) to a field of the intersected parcel. Next you check for the following vertex and if it intersects the same parcel as before you omit the editing for that parcel and continue. Each time the intersecting parcel is different you increase the counter.

2
  • If a person does not know how to loop over the vertices and check for intersection with the parcel layer with python, what then?
    – Taras
    Nov 23 '21 at 8:45
  • I will post an example of how this can be done Nov 27 '21 at 12:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.