6

The problem

In QGIS, it's easy to calculate areas of polygons or lengths of lines using QGIS expressions with area ($geometry) / $area or length($geoemtry) / $length respectively. However, there is no expression to measure angles based on the three points that define it.

The question

Given a polygon of any shape in QGIS, how can interior angle be measured apart from the manual measuring tool (toolbar icon)?

What I tried

What sounds quite easy at the first glance is more complicated. QGIS expressions offer a few angle-functions as:

  • line_interpolate_angle
  • angle_at_vertex
  • main_angle
  • azimuth

However, there is no get_angle_at_vertex(geometry,index) function so that you could calculate the angle of any polygon- or line-geometry at a certain vertex and neither an angle (point1,point2, point3) function.

Given the points 1, 2 and 3: how to get the angle at point 2, formed by the line connecting these three points on the left side in drawing-order?

enter image description here

3 Answers 3

17

The short answer:

To get the left-hand angle of any line or polygon at a vertex n, use this expression and replace 20 on line 3 with the index of the vertex at which you want to calculate the angle.

with_variable (
    'vertex',
    20,
    with_variable (
        'azimuth1',
        degrees (
            azimuth(
                point_n($geometry,@vertex-1),
                point_n($geometry,@vertex)
            )
        ),
        with_variable (
            'azimuth2',
            degrees (azimuth(
                point_n($geometry,@vertex),
                point_n($geometry,@vertex+1)
            )
        ),
        case
        when (@azimuth1 > @azimuth2) and (@azimuth1 > @azimuth2+180) then 540-@azimuth1+@azimuth2
        when (@azimuth1 > @azimuth2) then 180-@azimuth1+@azimuth2
        when (@azimuth1 < @azimuth2) and (@azimuth1+180>@azimuth2) then 180 + @azimuth2-@azimuth1
        when (@azimuth1 < @azimuth2) then @azimuth2-@azimuth1-180 
    end
        )
    )
)

The solution in more details:

In what follows, I will present different variants of this solution, based on what you need it for. The solution contains:

  • Get the angle at vertex n for a line or polygon (see above)
  • A. Geometric principles that this solution is based on
  • B. Get list of all the angles of a line or polygon
  • C. Create a separate value for each angle of a line or polygon

A. The geometric principles of the calculation of interior angles of a polygon and/or line

Let's consider an irregularily shaped polygon with all types of angles as in the screenshot below: weird shape (a monster?) with extreme angles. Let's first consider the geometric solution, than come to it's implementation in QGIS.

Screenshot: how to calculate all interior angles of this polygon: enter image description here

The angle values are distributes as shown in the following histogram, with min/max values at 5 and 354 degrees, respectively:

enter image description here

1. Get azimuth and what that is

For any line-segment, we can easily get the azimuth with azimuth(point_a,point_b). Azimuth means the north-based angle measured clockwise from the vertical for each line (start- to end-point). Consider it as the angle the minute hand of a clock forms: it is the angle formed by the line starting at 12 ("north"), going to the center of the clock (this vertical line remains) and from there to the head of the minute hand - the start- and end-points of this second part of the line are the arguments for the azimuth function. So at 5 past full hour (x:05 h), the azimuth is 30°, at a quarter past (x:15), it is 90°, at half past 180° (x:30), at 20 to (x:40) 240°

2. Formula for getting interior angle based an azimuths

Thus we can create for each line-segment the azimuth. Together with the azimuth of the next line, we can calculate the angle these two lines enclose. We want to get the angle at the left-hand side of the line (for right-hand-side: simply substract the result from 360). So for a polygon, we get the interior angle if their outer ring is drawn counter-clockwise (inside of the polygon is at left hand). The formula to calculate the angle depends on which of the azimuths (first or second line) is bigger and if the difference between both is more than 180° or not. The formula is as follows, whereby the azimuth of the first line is called a1, the azimuth of the second line a2:

  • if (a1 > a2) and (a1 > a2 + 180) then interior_angle = 540 -a1 + a2
  • if (a1 > a2) and (a1 < a2 + 180) then interior_angle = 180 - a1 + a2
  • if (a1 < a2) and (a1 + 180 > @azimuth2) then interior_angle = 180 + a2 - a1
  • if (a1 < a2) and (a1 + 180 < @azimuth2) then interior_angle = a2 - a1 - 180

Implementation: how to calculate interior angles of polygons with QGIS expressions:

B. Get list of all the angles of a line or polygon

A polygon or a line-string can only have as much different attributes as it contains features. So you can't directly create a different attribute for each angle (for this, see below). But you can easily create a list of all left-hand angles of a polygon or line using the expression that follows.

It creates a list of values as an array with as many values as there are vertices (for polygons: one angle at each vertex) respectively with 2 less values than there are vertices for lines: there are no angles at the start- and end-point of the line.

If you want the list in a string (text) format, simply enclose the expression in an array_to_string() statement. From this, you could also convert to a key/value map.

array_foreach (
    generate_series (
        1, 
         num_points( $geometry)-1 - 
        if (
            num_rings( $geometry) is null, 1,0
        )
    ),
    with_variable (
        'vertex',
        @element,
        with_variable (
            'azimuth1',
            degrees (
                azimuth(
                    point_n($geometry,@vertex),
                    point_n($geometry,@vertex+1)
                )
            ),
            with_variable (
                'azimuth2',
                degrees (azimuth(
                    point_n($geometry,@vertex+1),
                    point_n($geometry, if (@vertex+2>num_geometries( nodes_to_points( $geometry,true))+1,2,@vertex+2))
                )
            ),
            case
            when (@azimuth1 > @azimuth2) and (@azimuth1 > @azimuth2+180) then 540-@azimuth1+@azimuth2
            when (@azimuth1 > @azimuth2) then 180-@azimuth1+@azimuth2
            when (@azimuth1 < @azimuth2) and (@azimuth1+180>@azimuth2) then 180 + @azimuth2-@azimuth1
            when (@azimuth1 < @azimuth2) then @azimuth2-@azimuth1-180 
            end
            )
        )
    ) 
)

Screenshot: based on the expression above + a few additions for styling, I created a list of the interior angles of the polygon at each vertex and used this as a label. The vertices layer is only for visualization purpose (red dots, identifying the vertices to compare with the list at the right): enter image description here


C. Create a separate value for each angle

  1. If using polygons: Polygons to lines; if using lines: start with step 2.

  2. Explode lines. From now on, we work on this layer only.

You have now two possibilities: either jump directly to step 5 or take two intermediary steps ( 3 and 4) to first calculate the azimuth for the own and the next lines so that you could use them with the formula explained above. Step 5 however includes the calculation of these two azimuths, but because of this looks a bit more complex.

3. Get the azimuth for each line. As azimuth returns results in radians, convert to degrees. Create a new attribute field called azimuth_own with field calculator and this expression: degrees(azimuth (start_point ($geometry), end_point ($geometry))).

4. Get the azimuth for the next line: those that forms the angle with the current line. Again create a new filed called azimuth_next with this expression:

degrees (
    azimuth (
        start_point (
            geometry(
                get_feature_by_id(
                    @layer, 
                    $id+1
                )
            )
        ),
        end_point (
            geometry(
                get_feature_by_id(
                    @layer, 
                    $id+1
                )
            )
        )
    )
)
  1. Now, based on these two azimuths, you can calculate the angle at the left side of the line. For each line, the angle at it's end point (where it meets the start-point of the next line) is calculated. In the following expression, we (again) calculate azimuth_own and azimuth_next as variables, here called azimuth1 and azimuth2 according to the formula above. If you let out step 3 and 4, you can directly use this expression in Field Calculator to create a field inner_angle:
with_variable (
    'azimuth1',
    degrees (
        azimuth(
            start_point (geometry (get_feature_by_id(@layer,$id))),
            end_point (geometry (get_feature_by_id(@layer,$id)))
        )
    ),
    with_variable (
        'azimuth2',
        degrees (azimuth(
            start_point (geometry (get_feature_by_id(@layer,$id+1))),
            end_point (geometry (get_feature_by_id(@layer,$id+1)))
        )
    ),
    case
        when (@azimuth1 > @azimuth2) and (@azimuth1 > @azimuth2+180) then 540-@azimuth1+@azimuth2
        when (@azimuth1 > @azimuth2) then 180-@azimuth1+@azimuth2
        when (@azimuth1 < @azimuth2) and (@azimuth1+180>@azimuth2) then 180 + @azimuth2-@azimuth1
        when (@azimuth1 < @azimuth2) then @azimuth2-@azimuth1-180 
    end
    )
)

Visualization and control:

  1. To visualize (and control!) the results, you can create circle-segments to symbolize the angle based on the inner_angle value calculated. Use Geometry generator or Geometry by expression with this expression:
wedge_buffer (
    end_point($geometry),
    degrees (
        azimuth (
            end_point( $geometry),
            project (
                end_point ( $geometry),
                1000,
                radians (
                    angle_at_vertex( 
                        make_line (
                            start_point($geometry),
                            end_point ($geometry), 
                            end_point (
                                geometry (
                                    get_feature_by_id(
                                        @layer,
                                        $id+1
                                    )
                                )
                            )
                        ),
                        1
                    )-90
                )
            )
        )
    ),
    "inner_angle"
    ,
    200
)

Screenshot: blue circle segments show the angle that is additionally labeled with it's size. The polygon is draw counter clockwise: enter image description here

6
  • When I try your process using steps C5 & C6 there are two problems: 1. The angles calculated and displayed are the outer angles, and 2. The last record results in a NULL value for an angle. Thoughts?
    – Stu Smith
    Sep 21, 2021 at 3:54
  • I wrote the expression a few months before, so I'm not quite sure, but I guess 1. has to do with the direction the polygons are drawn. As stated, this solution is valid for "interior angle if their outer ring is drawn counter-clockweise (inside of the polygon is at left hand)". If counter-clockwise, simply calculate 360-[my expression]. 2. First and last vertex of a polygon are the same point and after the last vertex, there is no next vertex to draw a line - thus no angle.
    – Babel
    Sep 22, 2021 at 13:36
  • You're correct about #1; I missed the part in your answer about counter-clockwise (my apologies). But regarding #2, imagine an initial polygons that is roughly square; it therefore has four inner angles. What I'm seeing is that only the first three inner angles are symbolized.
    – Stu Smith
    Sep 23, 2021 at 14:52
  • OK, which one of the expressions did you use?
    – Babel
    Sep 23, 2021 at 14:54
  • When I use the first expression under "B. Get list of all the angles of a line or polygon", this is what I get: 4 angle values for a polygon with 4 vertices, see: i.sstatic.net/6lq79.png
    – Babel
    Sep 23, 2021 at 15:02
4

A solution for polygons

Let's assume there is a polygon layer 'poly_test' with a corresponding attribute table, see the image below

input

Step 1. Use the "Extract vertices" geoalogorithm step_1

Step 2. Procced with "Delete duplicate geometries" step_2

Step 3. In the Field Calculator create a new integer field

new_attribute

and fill it with the following expression (It may look creepy but it simply an implication of the 'Law of cosines' between current, previous and next vertex through a centroid of an intersection between original layer and a buffer around a vertex, see some explanations underneath).

Do not forget to (1) change the 'poly_test' into your actual polygon layer's name and (2) adjust 100 with an appropriate minimal size of a buffer around a vertex (It is a sort of an approximation factor, the lower the better but not always).

with_variable('buffer_size',100,
    with_variable('orignal_layer','poly_test',
        round(
            if("vertex_index" = minimum("vertex_index",group_by:="id"),
                -- first vertex
                degrees(
                    acos(
                        (
                        length(make_line($geometry,point_n(geometry(get_feature(@orignal_layer,'id',"id")),maximum("vertex_index",group_by:="id")+1)))^2
                        +
                        length(make_line($geometry,centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                        -
                        length(make_line(centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50))),point_n(geometry(get_feature(@orignal_layer,'id',"id")),maximum("vertex_index",group_by:="id")+1)))^2
                        )
                        /
                        (
                        2
                        *
                        length(make_line($geometry,point_n(geometry(get_feature(@orignal_layer,'id',"id")),maximum("vertex_index",group_by:="id")+1)))
                        *
                        length(make_line($geometry,centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))
                        )
                        )
                    )
                    +
                    degrees(
                        acos(
                            (
                            length(make_line($geometry,centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            +
                            length(make_line($geometry,point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index"+2)))^2
                            -
                            length(make_line(point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index"+2),centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            )
                            /
                            (
                            2
                            *
                            length(make_line($geometry,centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))
                            *
                            length(make_line($geometry,point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index"+2)))
                            )
                        )
                    ),
                -- last vertex
                if("vertex_index" = maximum("vertex_index",group_by:="id"),
                    degrees(
                        acos(
                            (
                            length(make_line($geometry,centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            +
                            length(make_line($geometry,point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index")))^2
                            -
                            length(make_line(point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index"),centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            )
                            /
                            (
                            2
                            *
                            length(make_line($geometry,centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))
                            *
                            length(make_line($geometry,point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index")))
                            )
                        )
                    )
                    +
                    degrees(
                        acos(
                            (
                            length(make_line($geometry,point_n(geometry(get_feature(@orignal_layer,'id',"id")),minimum("vertex_index",group_by:="id")+1)))^2
                            +
                            length(make_line($geometry,centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            -
                            length(make_line(centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50))),point_n(geometry(get_feature(@orignal_layer,'id',"id")),minimum("vertex_index",group_by:="id")+1)))^2
                            )
                            /
                            (
                            2
                            *
                            length(make_line($geometry,point_n(geometry(get_feature(@orignal_layer,'id',"id")),minimum("vertex_index",group_by:="id")+1)))
                            *
                            length(make_line($geometry,centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))
                            )
                        )
                    ),
                    -- everything in between
                    degrees(
                        acos(
                            (
                            length(make_line($geometry, point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index")))^2
                            +
                            length(make_line($geometry, centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            -
                            length(make_line(point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index"), centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            )
                            /
                            (
                            2
                            *
                            length(make_line($geometry, centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))
                            *
                            length(make_line($geometry, point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index")))
                            )
                        )
                    )
                    +
                    degrees(
                        acos(
                            (
                            length(make_line($geometry, centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            +
                            length(make_line($geometry, point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index"+2)))^2
                            -
                            length(make_line(point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index"+2), centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))^2
                            )
                            /
                            (
                            2
                            *
                            length(make_line($geometry, centroid(intersection(geometry(get_feature(@orignal_layer,'id',"id")), buffer($geometry,@buffer_size,50)))))
                            *
                            length(make_line($geometry, point_n(geometry(get_feature(@orignal_layer,'id',"id")),"vertex_index"+2)))
                            )
                        )
                    )
                )
            )
        )
    )
)

and get the final output

result


Explanations:

The value of interior angle achieved via a sum of two angles: ∠KAC + ∠KAB, where each is calculated using the 'Law of cosines' inside of corresponding triangle △AKB and △AKC.

explanation

Point (⋅)K is simply a centroid of an intersection between original layer and a buffer around corresponding vertex. More details about midpoint, center and centroid, please find here.

centroid(
    intersection(
        geometry(get_feature(@orignal_layer,'id',"id")),
        buffer($geometry,@buffer_size,50)
        )
    )
4
  • 1
    This works well and returns the same results as the other solution, so cross-verified that the basic considerations behind both solutions are correct. Still, if you don't mind, I consider the other approach to be less complex, more straightforward and universally valid for both polygons and lines.
    – Babel
    May 25, 2021 at 11:13
  • 1
    Thank you for confirming that my Frankenstein works! I do not mind because simple is better than complex.
    – Taras
    May 25, 2021 at 11:16
  • 1
    I must run your Frankenstein on my monster polygon: will see what results...?
    – Babel
    May 25, 2021 at 11:17
  • Dang! I followed your instructions to the letter - even digitized a polygon similar to yours - and all I get in the output is NULL. I tried buffer sizes from 0.1 to 1000 but no joy. Advice?
    – Stu Smith
    Sep 21, 2021 at 2:51
1

Currently working only for Single Part Polygons.

There is "LF Tools " Plugin, that possesses the "Calculate polygon angles" tool for that.

This algorithm calculates the inner and outer angles of the polygon vertices of a layer. The output layer corresponds to the points with the calculated angles stored in the respective attributes.

Input: input

Output: output

Additional application of the "Join attributes by location" geoalgorithm between points and original polygons can be useful.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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