# How to get maximum possible straight + perpendicular line through polygon?

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

• 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. 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

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

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

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

# 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

# Loads the output line layer to the legend
``````

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

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

# 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