Calculating interior angles of polygons or lines in QGIS

Question

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?

Answer 1

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:

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

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):

---

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

5. 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:

6. 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:

Answer 2

A solution for polygons

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

Step 1. Use the "Extract vertices" geoalogorithm

Step 2. Procced with "Delete duplicate geometries"

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

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

---

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.

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

Answer 3

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:

Output:

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