4

I'm trying to measure the maximum width of irregularly-shaped polygons. To me, the most logical way of doing this is to simply get the longest possible straight line through the polygon and then the longest possible line perpendicular to it (within the boundaries of the original polygon feature).

However, I can't figure out how to do either of these steps, and assumed the first would be pretty easy. The closest I've come so far is to get the length of the longest line using the minimum bounding geometry tool in Toolbox.

enter image description here

Also found this thread: The longest straight line (L) within a polygon and the longest perpendicular line (T) to L but the process is quite manual, and I need this calculation done on thousands of polygons.

Improve this question
3
  • I suggest you provide example data including the variety of shapes you need to calculate this for.
    – Erik
    Jul 22, 2022 at 9:20
  • Like the abstraction in the figure from the main post. Very irregular, many eggplanet-shaped. Assume there is virtually no symmetry.
    – hmnoidk
    Jul 22, 2022 at 9:24
  • see this post for finding the longest line (it can be simplified if you have no hole within the polygons)
    – JGH
    Jul 22, 2022 at 13:01

1 Answer 1

Reset to default
7

This python script will do what you want and you can play with the DENSIFY_EVERY_X_METER and EPSILON_ANGLE to either have a more accurate result or have your result faster

ex: enter image description here

import math
import numpy as np

DENSIFY_EVERY_X_METER = 0.05  # Distance max between two points of the polygon
# The process time scales heavily with this parameter so (the lower 
# the value the slower) so be careful or the process will takes hours
# you can start with a high value (if the perimeter of the polygon is
# 100m a good starting value would be 1m/0.5m) and try multiple 
# iterations with lower values

EPSILON_ANGLE = 5 # Angle in degrees maximum acceptable 
# between the longest line and the widest line

# Name of the layer containing the polygons you want to process
NAME_LAYER_IN_LEGEND = "polygon"

layer_polygon = QgsProject.instance().mapLayersByName(NAME_LAYER_IN_LEGEND)[0]


def azimuth(point1, point2):
    '''azimuth between 2 points (interval 0 - 360)'''
    angle = np.arctan2(point2.x() - point1.x(), point2.y() - point1.y())
    return np.degrees(angle) if angle >= 0 else np.degrees(angle) + 360

def angle_inter(tuple_pt1, tuple_pt2):
    '''Returns the angle between two lines
       tuple_pt1: first and last point of the first line
       tuple_pt2: first and last point of the second line 
    '''
    return abs(azimuth(*tuple_pt1) - azimuth(*tuple_pt2))

def find_lines(points):
    '''
    points : list of points that make up the densified polygon
    Returns two lines: 
        - the first line being the longest possible line 
            from one point to another in the 'points' list
        - the second line being the widest (longest line perpendicular 
            to the longest line) line between the all points in 
            'points' list
    '''
    geom_line_longest = None
    geom_line_width = None
    
    tuple_pt_longest = None
    distance_max = None
    for pt1 in points:
        # GeometryEngine is a optimized way to do spatial operations
        pt_eng1 = QgsGeometry.createGeometryEngine(QgsPoint(pt1))
        pt_eng1.prepareGeometry()
        
        for pt2 in points:
            # Calculate the distance between two points
            dist = pt_eng1.distance(QgsPoint(pt2))
            
            # If the two points are further appart than the previous
            # furthest points, the two points are saved 
            if distance_max is None or dist > distance_max:
                distance_max = dist
                tuple_pt_longest = (pt1, pt2)

    # Build the longest line from the two furthest points of the polygon
    geom_line_longest = QgsGeometry.fromPolylineXY(tuple_pt_longest)

    tuple_pt_width = None
    distance_max = None
    for pt1 in points:
        # GeometryEngine is a optimized way to do spatial operations
        pt_eng1 = QgsGeometry.createGeometryEngine(QgsPoint(pt1))
        pt_eng1.prepareGeometry()

        for pt2 in points:
            # Calculate the distance between two points
            dist = pt_eng1.distance(QgsPoint(pt2))

            # Calculate angle between the two current points and 
            # the longest line of the polygon
            tuple_pt_a_test = (pt1, pt2)
            angle = angle_inter(tuple_pt_longest, tuple_pt_a_test)

            # If the two points are further appart than the previous
            # furthest points, the two points are saved 
            # Only if the calculated angle between the two points is 
            # close to a 90° angle
            if (distance_max is None or dist > distance_max) \
                    and 90 - EPSILON_ANGLE <= angle <= 90 + EPSILON_ANGLE:
                distance_max = dist
                tuple_pt_width = (pt1, pt2)
    
    # Build the widest line from the two furthest points of the polygon
    # forming a line perpendicular to the longest line
    geom_line_width = QgsGeometry.fromPolylineXY(tuple_pt_width)
    return geom_line_longest, geom_line_width


layer_line = QgsVectorLayer("LineString?crs=EPSG:3857", "linesinpoly", "memory")

features = []
for polygon in layer_polygon.getFeatures():
    # Adds points to every segments of the polygon until there isn't
    # any two points more than DENSIFY_EVERY_X_METER meters away from each other
    geom_poly_densified = polygon.geometry().densifyByDistance(DENSIFY_EVERY_X_METER)

    # Extracts the points of the polygon into a set (list of unique items)
    points_poly = set(geom_poly_densified.asPolygon()[0])

    # Calculate the two lines from this set of unique points of the polygon
    geom_line_longest, geom_line_width = find_lines(points_poly)

    # Build features containing the two lines for the current polygon
    f_long = QgsFeature(layer_line.fields())
    f_long.setGeometry(geom_line_longest)
    f_width = QgsFeature(layer_line.fields())
    f_width.setGeometry(geom_line_width)
    # Adds the two features to the output list of line features
    features.extend([f_long, f_width])

# Adds the line features to the output layer
layer_line.dataProvider().addFeatures(features)

# Loads the output line layer to the legend
QgsProject.instance().addMapLayer(layer_line)

PS: If you don't want the lines to cross outside the polygon you can just add a contains verification between the poly and current line in both loops inside the if conditions like this:

enter image description here

import math
import numpy as np

DENSIFY_EVERY_X_METER = 0.05  # Distance max between two points of the polygon
# The process time scales heavily with this parameter so (the lower 
# the value the slower) so be careful or the process will takes hours
# you can start with a high value (if the perimeter of the polygon is
# 100m a good starting value would be 1m/0.5m) and try multiple 
# iterations with lower values

EPSILON_ANGLE = 5 # Angle in degrees maximum acceptable 
# between the longest line and the widest line

# Name of the layer containing the polygons you want to process
NAME_LAYER_IN_LEGEND = "polygon"

layer_polygon = QgsProject.instance().mapLayersByName(NAME_LAYER_IN_LEGEND)[0]


def azimuth(point1, point2):
    '''azimuth between 2 points (interval 0 - 360)'''
    angle = np.arctan2(point2.x() - point1.x(), point2.y() - point1.y())
    return np.degrees(angle) if angle >= 0 else np.degrees(angle) + 360

def angle_inter(tuple_pt1, tuple_pt2):
    '''Returns the angle between two lines
       tuple_pt1: first and last point of the first line
       tuple_pt2: first and last point of the second line 
    '''
    return abs(azimuth(*tuple_pt1) - azimuth(*tuple_pt2))

def find_lines(geom_poly_eng, points):
    '''
    geom_polygon: geometry of a polygon
    points : list of points that make up the densified polygon
    Returns two lines: 
        - the first line being the longest possible line 
            from one point to another in the 'points' list
        - the second line being the widest (longest line perpendicular 
            to the longest line) line between the all points in 
            'points' list
    '''
    geom_line_longest = None
    geom_line_width = None
    
    tuple_pt_longest = None
    distance_max = None
    for pt1 in points:
        # GeometryEngine is a optimized way to do spatial operations
        pt_eng1 = QgsGeometry.createGeometryEngine(QgsPoint(pt1))
        pt_eng1.prepareGeometry()
        
        for pt2 in points:
            # Calculate the distance between two points
            dist = pt_eng1.distance(QgsPoint(pt2))
            
            # If the two points are further appart than the previous
            # furthest points, the two points are saved 
            line_geom = QgsGeometry.fromPolylineXY([pt1, pt2])
            if (distance_max is None or dist > distance_max) \
                    and geom_poly_eng.contains(line_geom.constGet()):
                distance_max = dist
                tuple_pt_longest = (pt1, pt2)

    # Build the longest line from the two furthest points of the polygon
    geom_line_longest = QgsGeometry.fromPolylineXY(tuple_pt_longest)

    tuple_pt_width = None
    distance_max = None
    for pt1 in points:
        # GeometryEngine is a optimized way to do spatial operations
        pt_eng1 = QgsGeometry.createGeometryEngine(QgsPoint(pt1))
        pt_eng1.prepareGeometry()

        for pt2 in points:
            # Calculate the distance between two points
            dist = pt_eng1.distance(QgsPoint(pt2))

            # Calculate angle between the two current points and 
            # the longest line of the polygon
            tuple_pt_a_test = (pt1, pt2)
            angle = angle_inter(tuple_pt_longest, tuple_pt_a_test)

            # If the two points are further appart than the previous
            # furthest points, the two points are saved 
            # Only if the calculated angle between the two points is 
            # close to a 90° angle
            line_geom = QgsGeometry.fromPolylineXY([pt1, pt2])
            if (distance_max is None or dist > distance_max) \
                    and 90 - EPSILON_ANGLE <= angle <= 90 + EPSILON_ANGLE \
                    and geom_poly_eng.contains(line_geom.constGet()):
                distance_max = dist
                tuple_pt_width = (pt1, pt2)
    
    # Build the widest line from the two furthest points of the polygon
    # forming a line perpendicular to the longest line
    geom_line_width = QgsGeometry.fromPolylineXY(tuple_pt_width)
    return geom_line_longest, geom_line_width


layer_line = QgsVectorLayer("LineString?crs=EPSG:3857", "linesinpoly", "memory")

features = []
for polygon in layer_polygon.getFeatures():
    # Adds points to every segments of the polygon until there isn't
    # any two points more than DENSIFY_EVERY_X_METER meters away from each other
    geom_poly_densified = polygon.geometry().densifyByDistance(DENSIFY_EVERY_X_METER)

    # Extracts the points of the polygon into a set (list of unique items)
    points_poly = set(geom_poly_densified.asPolygon()[0])

    # Calculate the two lines from this set of unique points of the polygon
    geom_poly_eng = QgsGeometry.createGeometryEngine(geom_poly_densified.constGet())
    geom_poly_eng.prepareGeometry()
    geom_line_longest, geom_line_width = find_lines(geom_poly_eng, points_poly)

    # Build features containing the two lines for the current polygon
    f_long = QgsFeature(layer_line.fields())
    f_long.setGeometry(geom_line_longest)
    f_width = QgsFeature(layer_line.fields())
    f_width.setGeometry(geom_line_width)
    # Adds the two features to the output list of line features
    features.extend([f_long, f_width])

# Adds the line features to the output layer
layer_line.dataProvider().addFeatures(features)

# Loads the output line layer to the legend
QgsProject.instance().addMapLayer(layer_line)
Improve this answer
7
  • "[...] you can just add [...]" - you assume people know Python and/or easily understand your logic. While your answer may be correct and helpful, it currently not provides explanation nor furthers understanding. Could you explain a bit what your code does?
    – Erik
    Jul 22, 2022 at 10:29
  • @Erik Sure that's probably a good idea
    – Louis Cottereau
    Jul 22, 2022 at 11:30
  • And a time consuming one, so don't feel too obligated ;-)
    – Erik
    Jul 22, 2022 at 11:38
  • 1
    Updated, the best I could without spending two hours on it ^^
    – Louis Cottereau
    Jul 22, 2022 at 11:59
  • Thanks for this. Two questions: are you suggesting to start with DENSIFY_EVERY_X_METER = 2, for polygons around 100 m in perimeter? And EPSILON_ANGLE = 5 in your default variable. Shouldn't this be 90 degrees? Or is 5 degrees the tolerance (85 - 95 degrees are accepted)?
    – hmnoidk
    Aug 10, 2022 at 10:39

Your Answer

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

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